//	This program demonstrates a Markov model.  A simple two-state model of
//	a ligand-gated channel is created.  The ligand-gating is imposed simply
//	by setting the forward rate to a positive value initially, and setting
//	it to 0.0 after 100 ms.  This simulates a 100ms pulse of ligand.
//	The result should be a PSP-like "hump" of channels in the Open state.

// include the necessary header files
#include <iostream.h>			// for output
#include "Markov.h"				// Markov model class declaration

#ifdef macintosh
#include <console.h>			// for ccommand() function
#endif

// define global parameters
const real TMAX = 1.0;			// simulation end time (s)
const real DT = 0.005;			// simulation step size (s)
const int N = 2;				// define number of states
const int C = 0;				// define state names (for our convenience)
const int P = 1;				// (use "P" for oPen; "O" is confusing)
	
Markov *MakeMarkov( void )
{
	// create the demo Markov model
	
	Markov *m = new Markov( N );	// create an empty model
	if (m) {
		m->SetRate( C,P, 16.0 );		// set rate constants
		m->SetRate( P,C,  4.7 );
	
		m->SetValue( C, 1 );			// start with 100% in the Closed state
	}
	
	return m;
}

main( int argc, char **argv )
{
	// give Mac Symantec C++ users a chance to redirect input & output
	#ifdef macintosh
		ccommand( &argv );
	#endif
	
	// create the model
	Markov *model = MakeMarkov();		// model

	// check for failed creation; if good, do simulation
	if (model && model->GetQtyStates()) {

		// create a "ligand concentration" variable, and bind forward rate to it
		real Ligand = 1.0;					// "ligand" concentration
		model->SetCoeff( C,P, &Ligand );	// coefficient for Closed-->Open rate

		// set up the output stream
		cout.precision(4);
		cout.width(12);
		cout.setf(ios::showpoint);
		
		cout << "Time\tClosed\tOpen\n";

		// run the simulation
		for (real t=0; t<TMAX; t += DT) {
			cout << t << '\t' << model->GetValue(C) <<
						 '\t' << model->GetValue(P) << '\n';
 
			model->Step(DT);		// could have used gStepmaster.StepAll(DT)

			if (t>0.099 && t < 0.110) {
				// at this point, the ligand disappears...
				Ligand = 0.0;
			}
		}
	}

	return 1;
}