//----------------------------------------------------------------------------
//	Copyright (C) 2002, 2003, 2004 Humboldt-Universitaet zu Berlin
//
//	This library is free software; you can redistribute it and/or
//	modify it under the terms of the GNU Lesser General Public
//	License as published by the Free Software Foundation; either
//	version 2.1 of the License, or (at your option) any later version.
//
//	This library is distributed in the hope that it will be useful,
//	but WITHOUT ANY WARRANTY; without even the implied warranty of
//	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//	Lesser General Public License for more details.
//
//	You should have received a copy of the GNU Lesser General Public
//	License along with this library; if not, write to the Free Software
//	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//
//----------------------------------------------------------------------------
/**	\file Example_Continuous.cpp

	\author Ralf Gerstenberger

	\date created at 2003/07/02

	\brief Example for continuous simulations

	This example demonstrates the use of Continuous and	ContuTrace.

	\since 1.0
*/
/** \example Example_Continuous.cpp

	The basic continuous simulation techniques of ODEMx are introduced.

	Apart from discrete processes a simulation can also contain so called
	continuous processes. A discrete process does change the system state
	only at discrete moments in simulation time. It uses the
	actions defined in its 'main()' function. The progress of time for a
	process is realised with functions like 'holdFor', 'activateAt' and so on.
	See the basic simulation example for details.
	A continuous process instead changes the state of a system continuously.
	Such a process could be for example the melting of a block of metal
	inside an oven.
	At every time, during the melting process, the temperature of the
	metal is changed. That means \b every time you measure the temperature you
	will get a new value.
	Of course real continuity is quite impossible in the discrete world of
	our computers. So ODEMx has to approximate continuous state changes with
	a step by step computation.
*/

#include <odemx/base/Continuous.h>
#include <odemx/base/Simulation.h>
#include <odemx/base/DefaultSimulation.h>
using namespace odemx::base;

#ifdef ODEMX_USE_CONTINUOUS

#include <cmath>

//
// The first Continuous process is very simple. It only defines
// a sinus/co-sinus oscillator.
//
class Oscillator : public Continuous {
public:
	//
	// The construction of a Continuous object is a little different to
	// that of a Process object. You also have to provide the number
	// of state-variables used by your process. In case of Oscillator
	// we need two variables.
	//
	Oscillator() : Continuous( odemx::getDefaultSimulation(), "Oscillator", 2 ) {};

protected:
	//
	// The main-function of a continuous process is quite similar to that of
	// a discrete process. You can do everything that is possible in Process.
	// That's why a continuous can behave just like a discrete process. But
	// it can also go through phases of continuous state changes.
	//
	virtual int main() {
		//
		// Before you can start the solver, which is computing the continuous
		// state changes, you will have to initialise the state variables.
		//
		state[0]=1;
		state[1]=0;

		//
		// Than you should set some parameters for the internal solver. These
		// parameters include error sensitivity and the step length.
		//
		// The step length is set with 'setStepLength()'. The first parameter
		// sets the minimal step length, while the second defines the maximal
		// step length. The internal solver will compute new states in steps.
		// Between every step the process holds. The actual length of the step
		// taken will depend on numerical errors, peer processes, state events,
		// and the parameters provided with 'setStepLength()'. The actual step
		// length will not exceed your provided maximum. If the step length has
		// to be reduced because of numerical errors or state events it will not
		// be reduced below your provided minimum.
		//
		setStepLength(0.01, 0.1);

		//
		// The error sensitivity is set with 'setErrorlimit()'. The first parameter
		// defines whether the errors should be considered relative to the value of
		// the state variables (0) or absolute (1). The second parameter
		// sets the maximum error acceptable (relative or absolute). If the actual
		// error exceeds this value the solver will try to reduce the step length.
		// If this fails, because the step length is already to small, you will get
		// a simulation error.
		//
		setErrorlimit(0, 0.1);

        	//
		// Finally, the continuous phase is started with 'integrate()'. It is stopped
		// either by a time event, a state event or an interrupt from another process.
		//
		// The time event is set by the first parameter. If it is 0 the solver will
		// run for ever. Otherwise it will run to the given absolute time. If the
		// provided time has already passed it will return at once. The return value of
		// integrate will be 0 if the time event was hit. If the process is interrupted
		// 'integrate()' returns 2.
		//
		integrate(20.0, 0);

		return 0;
	}

	//
	// Every Continuous process has to provide its specific 'derivatives()' function.
	// In this function you define how the state changes during the time. You do this
	// by setting the rate in which a state variable is changed. Although you don't have
	// to, you can include the time provided by t in your computation. But never use
	// 'getCurrentTime()' to get the time.
	//
    virtual void derivatives (double t) {
		rate[0]=-state[1];
		rate[1]=state[0];
	}
};

//
// The FreeFall continuous process needs only one state variable. It demonstrates
// the use of parameter t in 'derivatives()'. The result is an idealistic free fall.
//
class FreeFall : public Continuous {
public:
	FreeFall() : Continuous( odemx::getDefaultSimulation(), "FreeFall", 1) {};

protected:
	virtual int main() {
		state[0]=0.0;

		setStepLength(.1,1);
		integrate(20.0, 0);

		return 0;
	}

	//
	// Again, never use getCurrentTime() to include the current time in your
	// computation. The reason for this is, that 'derivatives()' is called multiple
	// times for each step and these calls are not synchronised to the
	// 'official' time in the simulation.
	//
    	virtual void derivatives (double t) {
		double g=9.81;

		rate[0]=t*g;
	}
};

//
// RealFall simulates a 'real fall' which is slowed down by friction.
//
class RealFall : public Continuous {
public:
	RealFall() : Continuous( odemx::getDefaultSimulation(), "RealFall", 2 ) {};

protected:
	//
	// As in FreeFall we don't set the error limits. We can do so because
	// ODEMx uses default settings for error limits and step length. The
	// default for the error limits is a relative (0) error limit of 0.1 .
	//
	virtual int main() {
		state[0]=2.0;
		state[1]=0.0;

		setStepLength(.1, 1);
		integrate(20.0, 0);

		return 0;
	}

    virtual void derivatives (double t) {
		double k=0.5;
		double g=-9.81;

		rate[0]=state[1];
		rate[1]=g - k*state[1];
	}
};

//
// RealBounce finally demonstrates the use of state events.
//
class RealBounce : public Continuous {
public:
	RealBounce() : Continuous( odemx::getDefaultSimulation(), "RealBounce", 3 ) {};

	//
	// 'hitGround()' is used to check a state event. The signature
	// of functions that can be used as state-event-functions is:
	//     bool(Process::*)()
	// ODEMx uses pointer to member functions if it needs a call-back.
	// The advantage is, you can use member functions with full access
	// to internal data of your classes to check state events. The
	// costs of this design decision is that you will always have to
	// cast the address of your state-event-function to the type 'Condition'.
	//
	// A state function has to return true if a state event has occurred.
	// Remember, the computation of state changes is done in steps. Because
	// of that it is unlikely to hit a state event exactly. This has to
	// be considered when programming a state-event-function. If a state
	// event has occurred (or passed) the internal solver starts a binary
	// search to get closer to the exact event time. Finally, if it gets close
	// enough (minimum step length) the computation is stopped and 'integrate()'
	// returns 1.
	//
	bool hitGround() {
		return state[0]<=0.0;
	}

protected:
	virtual int main() {
		double g=-9.81;

		state[0]=2.0;
		state[1]=0.0;
		state[2]=g;

		setStepLength(0.01, 0.1);
		//
		// RealBounce uses the return value of 'integrate()' to control
		// the computation. Remember 'integrate()' returns 0 if a
		// time event occurred, 2 if the process was interrupted and 1
		// if a state event was detected.
		// 'integrate()' is called with the time event 20.0 and the
		// state-event-function 'hitGround()'. If the state event
		// is hit we reflect the movement 'state[1]=-(0.8*state[1])'
		// and continue until the simulation time exceeds 20.0.
		//
		while (integrate(20.0, (odemx::base::Condition)&RealBounce::hitGround)==1)
		{
			//
			// The ball hit the ground and is reflected.
			//
			if (fabs(state[1])<0.01 && state[0]<0.01) {
				//
				// We must stop the ball if it has lost to much energy.
				// Otherwise we would produce an annoying loop. The
				// state event would be hit in every step.
				//
				state[1]=0.0;
				state[2]=0.0;
			}
			state[1]=-(0.8*state[1]);
		}

		return 0;
	}

    virtual void derivatives (double t) {
		double k=0.5;

		rate[0]=state[1];
		rate[1]=state[2] - k*fabs(state[1]);
		rate[2]=0.0;
	}
};

int main(int argc, char* argv[]) {

	Oscillator* osci = new Oscillator();
	FreeFall* free = new FreeFall();
	RealFall* real = new RealFall();
	RealBounce* bounce = new RealBounce();

	//
	// To follow the state changes we use the class ContuTrace.
	// ContuTrace observes a provided continuous process and logs all
	// state changes into a text file. The file is managed by ContuTrace.
	// The name is either set in the constructor or build from the
	// name of the observed continuous process. If you run the simulation
	// you will find the 4 text files:
	//     Oscillator_trace.txt
	//     FreeFall_trace.txt
	//     RealFall_trace.txt
	//     RealBounce_trace.txt
	//

	ContuTrace osciTrace(osci), freeTrace(free), realTrace(real), bounceTrace(bounce);

	//
	// Continuous processes are activated just like discrete processes.
	//
	osci->activate();
	free->activate();
	real->activate();
	bounce->activate();

	//
	// We can run the simulation without a time limit because all
	// our continuous processes end at 20.0 .
	//
	odemx::getDefaultSimulation().run();

	delete bounce;
	delete real;
	delete free;
	delete osci;

	return 0;
}

#else // ODEMX_USE_CONTINUOUS not defined

#include <iostream>

int main( int argc, char* argv[] )
{
	std::cerr << "ODEMx compiled without support for continuous processes."
			  << std::endl;
	return 0;
}

#endif /* ODEMX_USE_CONTINUOUS */
