//	This is similar to the Hodgkin-Huxley demo, but it uses
//	two compartments, joined by a Delay link.

// 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
#include "Delay.h"		// for Delay link

#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 COMP_L = 30E-6;		// cylinder length (m)
const real COMP_D = 30E-6;		// cylinder diameter (m)

const real TMAX = 0.050;		// simulation time (sec)
const real DT =   0.0001;		// simulation time step (sec)
const real DELAY= 0.015;		// delay time (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 *MakeCompartment( void )
{
	// create cylinder by radius and length
	Cylinder *comp = new Cylinder(COMP_D/2, COMP_L);

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

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

	return comp;
}

main( int argc, char **argv )
{
	// give Mac Symantec C++ users a chance to redirect input & output
	#ifdef macintosh
		ccommand( &argv );
	#endif
	
	// create the compartments and join them with a one-way delay
	Cylinder *comp1 = MakeCompartment();
	Cylinder *comp2 = MakeCompartment();
	Delay delayLink( comp1, comp2, comp1->Gi, (int)(DELAY/DT) );

	// attach an injector to comp1, but don't apply current yet
	Injector inject( comp1, 0.0 );
	
	// set up the output stream
	cout.precision(4);
	cout.width(7);
	
	// pre-run the simulation to let initial conditions settle down
	for (real pret=0; pret<0.050; pret += DT)
		gStepmaster.StepAll(DT);

	// run the simulation and generate output
	cout << "msec" << DELIM << "mV1" << DELIM << "mV2" << endl;
	for (real t=0; t<TMAX; t += DT) {
		cout << t*1000 << DELIM << comp1->GetV()*1000 
			<< DELIM << comp2->GetV()*1000 << endl;
		gStepmaster.StepAll(DT);

		// apply the current pulse between t=(0.010,0.015)
		if (t>=0.010 && t<=0.015) inject.SetI( 0.3E-9 );
		else inject.SetI( 0.0 );
	}
	
	return 1;
}