• Backlinks
• Edit
• History
• Print

TstepEM.c

 /* Tstep Module for the Explicit Monodomain case */

 /* Copyright John B. Pormann, 21 July 2000, all rights reserved */


 /* you can link directly to the documentation with: */
 /*			Dev.TstepEM				*/
 /*			User.TstepEM			*/


 #include "CardioWave.h"

 /* global variable */
 sim_t SimType = Monodomain;

 static char* RCSID = "$Id: TstepEM.c,v 1.5 2004/08/24 19:52:57 jpormann Exp $";

 /* local variables */
 static real   TT = 0.0;
 static real   Tf = 5.000;
 static real   dT = 0.001;
 static vector Fv = _INIT_VECTOR;
 static vector Fq = _INIT_VECTOR;

 static rword resources[] = {
 	{ "deltat",	2001 },
 	{ "dt",		2001 },
 	{ "t0",		2003 },
 	{ "tf",		2002 },
 	{ "tfinal",	2002 },
 	{ "tinit",	2003 },
 	{ "tinitial",	2003 },
 	{ NULL, 0 }
 };

 int InitTimeStepper( char** res, real* t ) {
 	int   i,cmd;

 	i = 0;
 	while( res[i] != NULL ) {
 		cmd = FindCommand( resources, res[i] );
 		switch( cmd ) {
 			case 2001:
 				dT = GetRealValue( res[i] );
 				break;
 			case 2002:
 				Tf = GetRealValue( res[i] );
 				break;
 			case 2003:
 				TT = GetRealValue( res[i] );
 				break;
 		}
 		i++;
 	}

 	CreateMatrix();

 	if( vecalloc(&Fv,Tissue) < 0 ) {
 		return( -1 );
 	}
 	if( vecalloc(&Fq,State) < 0 ) {
 		return( -2 );
 	}

 	*t = TT;

 	if( (DebugLevel>0) and (SelfPE==0) ) {
 		printf("TstepEM?: dt = le\n",dT,Tf);
 		if( DebugLevel > 1 ) {
 			printf("TstepEM?: RCSID: %s\n",RCSID);
 		}
 	}

 	return( 0 );
 }

 /* what to do at exit of program */
 void ExitTimeStepper( void ) {
 	if( SaveState and (SelfPE==0) ) {
 		fprintf( FpResources, "t0=%le\n", TT );
 	}

 	return;
 }

 int CreateMatrix( void ) {
 	/* : for Explicit-Monodomain, no work needs be done   */
 	/*   on any matrices, but we do need to compute 1/Vol */
 	if( InvVolume.size == 0 ) {
 		/* : since we don't really need Volume directly we'll */
 		/*   just swap instead of allocating and copying      */
 		vswap( &Volume, &InvVolume );
 		vrecip( 1.0, InvVolume );
 	}

 	return( 0 );
 }

 int Update( real* t, vector Vm, vector Vx, vector Q, vector Av ) {
 	/* update time count */
 	TT += dT;
 	*t = TT;

 	/* compute the diffusion currents */
 	sprvec( 1.0/Beta, MatrixINT, Vm, 0.0, Fv );

 	/* : adjust for voronoi volumes/areas */
 	vvscale( 1.0, InvVolume, Fv );

 	/* : now, add the ionic currents into Fv & Fq */
 	GetF( Vm, Q, Fv, Fq, Av );

 	/* : this applies 1.0*Istim */
 	ApplyStimulus( Intracellular, TT, 1.0, Fv, Fq );

 	/* update Vm & Q */
 	vvadd( dT/Cm, Fv, Vm );
 	vvadd( dT, Fq, Q );

 	if( TT >= Tf ) {
 		return( -1 );
 	}

 	return( 0 );
 }