//	This program simulates a passive dendritic cable.

// include the necessary header files
#include <iostream.h>			// for output
#include "Cylinder.h"			// Compartment class declaration
#include "Link.h"				// Link class declaration
#include "Injector.h"			// Injector class declaration

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

// define global parameters
const real TMAX = 0.1;				// simulation end time (s)
const real DT = 0.005;				// simulation step size (s)
const real Ri = 0.01;				// specific internal (axial) resistance (ohm m)
const real Len = 0.00001;			// cable length (m)
const int  NComp = 10;				// number of compartments
const real EREST = -0.07;			// resting membrane potential (volts)
const real ELEAK = EREST + 0.0106;	// "leak" potential (volts)


void Connect( Compartment& pComp1, Compartment& pComp2)
{
	//	connect the two given compartments via a pair of Links
	
	// center-to-center conductance is mean of axial conductance of each compartment
	real conductance = (pComp1.Gi + pComp2.Gi) / 2.0;
	
	Link *temp = new Link( &pComp1, &pComp2, conductance );
	temp = new Link( &pComp2, &pComp1, conductance );
}
	
main( int argc, char **argv )
{
	// give Mac Symantec C++ users a chance to redirect input & output
	#ifdef macintosh
		ccommand( &argv );
	#endif
	
	// set up the output stream
	cout.precision(4);
	cout.width(7);
	cout.setf(ios::showpoint);
	
	// create the cable
	
	int cmpnum;							// compartment counter
	Cylinder *compArray;				// array of cylindrical compartments
	
	compArray = new Cylinder[NComp];				// create the cylinders
	for (cmpnum=0; cmpnum<NComp; cmpnum++) {		// for each one...
		compArray[cmpnum].SetLength( Len / (real)NComp );		// set the length
		compArray[cmpnum].SetV( ELEAK );			// and the initial voltage
		if (cmpnum > 0)								// then connect to neighbors
			Connect( compArray[cmpnum-1], compArray[cmpnum]);	
	}

	// provide a current injection of 0.3 uA in the first compartment
	Injector inject( &compArray[0], 0.3E-6 );      

	// now the cable's created; run the simulation and see what happens!
	
	for (real t=0; t<TMAX; t += DT ) {
		// update the simulation
		gStepmaster.StepAll(DT);
		
		// output the voltage of each compartment
		for (cmpnum=0; cmpnum<NComp; cmpnum++) {
			cout << compArray[cmpnum].GetV();
			cout << (cmpnum+1 < NComp ? '\t' : '\n');
		}
	}
	
	return 1;
}