//	This program demonstrates Cylinder, Injector, and StdChan
//	by creating a compartment with standard Hodgkin-Huxley
//	channels and injecting a current, resulting in a train
//	of action potentials.  Reproduces Genesis "tutorial3.g".

// include the necessary header files
#include <iostream.h>		// for output
#include "Cylinder.h"		// for cylindrical compartment
#include "Injector.h"		// for current injection
#include "ChanStd.h"		// for Hodgkin-Huxley channels

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

// define some global paramaters -- all in SI units

const real RM = 0.33333;		// specific membrane resistance (ohm m^2)
const real CM = 0.01;			// specific membrane capacitance (farad/m^2)
const real RA = 0.3;			// specific axial resistance (ohms m)
const real EREST = -0.07;		// resting membrane potential (volts)
const real ELEAK = EREST + 0.0106;	// "leak" potential (volts)
const real ENA = 0.045;			// sodium equilibrium potential
const real EK = -0.082;			// potassium equilibrium potential

const real SOMA_L = 30E-6;		// cylinder length (m)
const real SOMA_D = 30E-6;		// cylinder diameter (m)

const real TMAX = 0.100;		// simulation time (sec)
const real DT =   0.0001;		// simulation time step (sec)

const char DELIM = '\t';		// field delimiter for output

void SetUpNaChannel( ChanStd &pChan )
{
	// set the equilibrium potential
	pChan.SetE( ENA );

	// set up "M" (activation) gating variable
	pChan.Mexp = 3;
	pChan.SetFunc( kStdAlphaM, kLinForm, -0.1E6, -0.010, EREST + 0.025 );
	pChan.SetFunc( kStdBetaM, kExpForm, 4.0E3, -18.0E-3, EREST );
	
	// set up "H" (inactivation) gating variable
	pChan.Hexp = 1;
	pChan.SetFunc( kStdAlphaH, kExpForm, 70.0, -20.0E-3, EREST );
	pChan.SetFunc( kStdBetaH, kSigForm, 1.0E3, -10.0E-3, EREST + 30.0E-3 );
}

void SetUpKChannel( ChanStd &pChan )
{
	// set the equilibrium potential
	pChan.SetE( EK );

	// set up "M" (activation) gating variable
	pChan.Mexp = 4;
	pChan.SetFunc( kStdAlphaM, kLinForm, -10.0E3, -10.0E-3, EREST + 10.0E-3 );
	pChan.SetFunc( kStdBetaM, kExpForm, 125.0, -80.0E-3, EREST );
	
	// set up "H" (inactivation) gating variable
	pChan.Hexp = 0;
	// (since exponent is 0, no need to set the functions)
}

Cylinder *MakeDemoSoma( void )
{
	// create cylinder by radius and length
	Cylinder *soma = new Cylinder(SOMA_D/2, SOMA_L);

	// copy membrane parameters from global variables
	soma->SetE( ELEAK );
	soma->SetRm( RM );
	soma->SetCm( CM );
	soma->SetRa( RA );
	
	soma->SetV( ELEAK );		// start at rest

	// add Hodgkin-Huxley Na and K channels
	ChanStd *NaChan = new ChanStd( soma, 1200*soma->GetArea() );
	SetUpNaChannel( *NaChan );
	
	ChanStd *KChan = new ChanStd( soma, 360*soma->GetArea() );
	SetUpKChannel( *KChan );
	
	return soma;
}

main( int argc, char **argv )
{
	// give Mac Symantec C++ users a chance to redirect input & output
	#ifdef macintosh
		ccommand( &argv );
	#endif
	
	// create the soma
	Cylinder *soma = MakeDemoSoma();
	
	// provide a current injection of 0.3 nA
	Injector inject( soma, 0.3E-9 );
	
	// set up the output stream
	cout.precision(4);
	cout.width(7);
	
	// run the simulation
	cout << "msec" << DELIM << "mV" << endl;
	for (real t=0; t<TMAX; t += DT) {
		gStepmaster.StepAll(DT);
		cout << t*1000 << DELIM << soma->GetV()*1000 << endl;
	}
	
	return 1;
}