/* Modified extensively 1986-1989 Timothy L. Davis
 * $Source: /mit/hst/src/sim/RCS/eqns.c,v $
 * $Log:	eqns.c,v $
 * Revision 2.2  90/08/28  14:44:37  tldavis
 * Added function UpdateTiming() to calculate the portion of the cardiac
 * cycle and store it in an unsigned int.
 * 
 * Revision 1.1  90/08/19  16:25:52  tldavis
 * Initial revision
 * 
 * Revision 2.1  90/08/18  18:39:51  tldavis
 * New location of CVDefs in this directory.
 * 
 * Revision 2.0  90/08/06  18:02:38  tldavis
 * Final X11R3 version.
 * 
 * Revision 1.6  90/08/06  17:03:31  tldavis
 * Added a few new dependent variables: atrial pressures and instantaneous
 * ventricular capacitances.
 * 
 * Revision 1.5  90/05/31  16:25:16  tldavis
 * Pre-ACIS checkin.
 * 
 * Revision 1.4  90/02/20  00:01:55  tldavis
 * Fixed ventricular constant volume bug.
 * 
 * Revision 1.3  89/10/12  22:17:01  tldavis
 * Updated Fall, 1989.
 * 
 * Revision 1.2  89/08/21  15:23:03  tldavis
 * Works with baroreflex control in a single process.
 * 
 * Tweaked X11R3 version, for cpoon, with control system features.
 * 
 * Equations within the circulatory system model including:
 *    -	differential KCL equations at each node
 *    -	heart capacitance values as a function of time via a look-up table
 *       (without interpolation).  NOTE: this assumes the diastolic capacitance
 *	is constant, which is only valid at normal ventricular volumes.
 *
 */
#include "CVDefs.h"
#include <stdio.h>
#include <math.h>
#include "../event/Event.h" /* for UpdateParamWindow() */

extern double sigHR, sigR, sigV, sigCl, sigCr; /* from reflex.c */
extern void calcReflexes(); /* reflex.c */

/* 
 * 1/3 of systemic arterial capacitance is in the thorax: A_TH_FRAC
 */
#define A_TH_FRAC 0.333

/* nonstatic variables below are used by reflex.c exclusively */
static double Vl0;	/* left heart volume at 0 pressure */
static double Vr0;	/* right heart volume at 0 pressure */
static double Va0;	/* peripheral arterial volume at 0 pressure */
double Vv0;	/* peripheral venous volume at 0 pressure */
static double VPa0;	/* pulmonary arterial volume at 0 pressure */
static double VPv0;	/* pulmonary venous volume at 0 pressure */
static double bv;	/* total blood volume */
double Ra;	/* peripheral arterial resistance */
static double Rv;	/* peripheral venous resistance */
static double RPv;	/* pulmonary vascular resistance */
static double Rli;	/* left ventricular inflow resistance */
static double Rlo;	/* left ventricular outflow resistance */
static double Rro;	/* right ventricular outflow resistance */
double Ca;	/* peripheral arterial capacitance */
static double Cv;	/* peripheral venous capacitance */
static double CPa;	/* pulmonary arterial capacitance */
static double CPv;	/* pulmonary venous capacitance */
static double Cldias;	/* left heart diastolic capacitance */
double Clsys;	/* left heart systolic capacitance */
static double Crdias;	/* right heart diastolic capacitance */
double Crsys;	/* right heart systolic capacitance */
static double pth;	/* intrathoracic pressure */
double hr;	/* heart rate operating point in beats per minute */
static int DLI;		/* left heart inflow diode */
static int DLO;		/* left heart outflow diode */
static int DRI;		/* right heart inflow diode */
static int DRO;		/* right heart outflow diode */

/* The above variables are the signals used for calculations involving
 * hr, Ra, Vv0, Clsys, and Crsys.  The nominal values of hr, Ra, etc.
 * are used as setpoints for the linear control system if reflexes are
 * on.  If reflexes are off, hr etc. are copied directly into sigHR etc. 
 */

/*PUBLIC*/
void initeqnvars(p)
SIMULATION *p;
{
  Vl0 = p->Vl0;  Vr0 = p->Vr0;  Va0 = p->Va0;  Vv0 = p->Vv0;  
  VPa0 = p->VPa0;  VPv0 = p->VPv0;  bv = p->bv;
  Ra = p->Ra;Rv = p->Rv;RPv = p->RPv;Rli = p->Rli; Rlo = p->Rlo; Rro = p->Rro;
  Ca = p->Ca;  Cv = p->Cv; CPa = p->CPa; CPv = p->CPv;
  Cldias = p->Cldias; Clsys = p->Clsys; Crdias = p->Crdias; Crsys = p->Crsys;
  pth = p->pth;  hr = p->hr;
}

/*PUBLIC*/
SIMULATION *presentsim(p)
SIMULATION *p;
{
  p->Vl0 = Vl0;  p->Vr0 = Vr0;  p->Va0 = Va0;  p->Vv0 = Vv0;  p->VPa0 = VPa0;
  p->VPv0 = VPv0;  p->bv = bv;  p->Ra = Ra;  p->Rv = Rv;  p->RPv = RPv;
  p->Rli = Rli;  p->Rlo = Rlo;  p->Rro = Rro;
  p->Ca = Ca;  p->Cv = Cv;  p->CPa = CPa;  p->CPv = CPv;
  p->Cldias = Cldias;  p->Clsys = Clsys;  p->Crdias = Crdias; p->Crsys = Crsys;
  p->pth = pth;  p->hr = hr;
  return(p);
}

/*PUBLIC*/
void
eqns(dxdt, i, x) /* takes x[] pressures, and time index in systole
		  * then return changes in dxdt[]
		  */
double x[], dxdt[];
int i;
{
	double Cl(), Cr(), dCldt(), dCrdt(), Cinc();
	double q0, q1, q2, q3, q4, q5; /* These are FLOWs */

	q0 = DLI ? (xpPv - xpl) /Rli  : 0.;
	q1 = DLO ? (xpl - xpa)  /Rlo  : 0.;
	q2 =       (xpa - xpv)  /sigR;
	q3 = DRI ? (xpv - xpr)  /Rv   : 0.;
	q4 = DRO ? (xpr - xpPa) /Rro  : 0.;
	q5 =       (xpPa - xpPv)/RPv;
	dpldt = ((pth - xpl)*dCldt(i) + q0 - q1)/Cl(i);
	dpadt = (q1 - q2)/Ca;
	dpvdt = (q2 - q3)/Cv;
	dprdt = ((pth - xpr)*dCrdt(i) + q3 - q4)/Cr(i);
	dpPadt = (q4 - q5)/CPa;
	dpPvdt = (q5 - q0)/CPv;
}

static	struct	ctable         /* Table of heart contraction capacitance */
{
	double	fctr; /* Maximum capacitance for scaling purposes. Min == 0 */
	unsigned long val[TBLSEGSPLUS2];
} ctbl = {
#include "tblc"
};

static	struct	dctable        /* table of heart contraction cap. derivatives*/
{
	double	fctr; /* Not used. Calculated from ctbl.fctr * ds/dt */
	long val[TBLSEGSPLUS2];
} dctbl = {
#include "tbldcdt"
};

static double 
  k2oTsys,     /* time multiplier for Sys. table lookup*/
  kr1,         /* capacitance multiplier for RH systole*/         
  kl1,         /* capacitance multiplier for LH systole*/
  kr2,         /*     "          "        "  RH systole dC/dt */
  kl2;         /*     "          "        "  LH    "      "   */

/*PUBLIC*/
double Cinc(dt)
     double dt;
{
  return (dt * k2oTsys);
}

/*PUBLIC*/
double Cr(i)
     int i;
{
  if (i >= TBLSEGS)  return (Crdias);
  else  return (sigCr + ctbl.val[i] * kr1); /* ctbl.val is LONG! */
}

/*PUBLIC*/
double	Cl(i)
     int i;
{
  if (i >= TBLSEGS) return (Cldias);
  else return (sigCl + ctbl.val[i] * kl1);
}

/*PUBLIC*/
double	dCrdt(i)
     int i;
{
  if (i >= TBLSEGS) return (0.);
  else return (dctbl.val[i] * kr2);
}

/*PUBLIC*/
double dCldt(i)
     int i;
{
  if (i >= TBLSEGS) return (0.);
  else return (dctbl.val[i] * kl2);
}

/*PUBLIC*/
double Tsystole(hrate)
     double hrate;
{
  return(0.30*sqrt(60./hrate));
}

/* 
 * lookup tables for heart capacitance used from 0 <= t <= x1Tsys.
 * k2oTsys is the number of lookup table items per time in seconds.
 * kr1 is scaling factor for right heart capacitance.
 * kr2 is scaling factor for right heart d(capacitance)/dt in cap. per sec.
 * similar for kl1 and kl2.  
 */
/*PUBLIC*/
void btconst() 
{
  double Tsys = Tsystole(sigHR); /* time for contraction (seconds) */; 
  k2oTsys = SYSTOLESEGS/Tsys;  /* number of data file segs per time*/
  kl1 = (Cldias - sigCl)/ctbl.fctr; /* scaling factor for left heart */
  kl2 = kl1 * k2oTsys;  /* Scaling factor for dcdt.  Table is in
			 * dc/ds (per segment).  Convert dc/dt to cap 
			 * per time by chain rule 
			 */
  kr1 = (Crdias - sigCr)/ctbl.fctr; /* scaling factor for right heart */
  kr2 = kr1 * k2oTsys;  /* scaling factor for dcdt on right */
}


/*PUBLIC*/
void initvalves(x)
double x[];
{
	DLI = (xpPv >= xpl);
	DLO = (xpl >= xpa);
	DRI = (xpv >= xpr);
	DRO = (xpr >= xpPa);
}

/*PUBLIC*/
void fixvolume(x)
double x[];
{
      /* find deviation from total blood vol, 
       * correct at systemic venous reservoir 
       */
  double diff = (bv - (Vl0+Va0+sigV+Vr0+VPa0+VPv0) -
		 (xCl*(xpl-pth) + Ca*(xpa-A_TH_FRAC*pth) + Cv*xpv +
		  xCr*(xpr-pth) + CPa*(xpPa-pth) + CPv*(xpPv-pth)))/Cv;
#ifdef DEBUG
  if (fabs(diff) > 1.0)
    printf("Volume error:  %g ml   ",diff);
#endif
  xpv += diff;
}

/*PUBLIC*/
double VenVol(x)
     double x[];
{
  return(xpv * Cv);
}

/* Returns TRUE for valves state within given period if pd1 & pd2 are different
 * Returns TRUE at leading edge of pd1 if pd2 == 0.
 * Sets bit INPERIOD_POS on entrance, clears on exit.
 * FALSE otherwise.
 * BIT CODE:  0, 1: DLI, DLO beginning desired period, immediately following
	      2, 3: the previous DLI, DLO to bits 0, 1.
	      4, 5: DLI, DLO immediately following the end of the periods, with
	      6, 7: DLI, DLO the last valve state in the period.
	      8:    Flag set in this routine 1: in period, 0: not in.
	      9, 10: The current flag settings saved here from last time.
	      11:   1 if averages are being plotted instead of direct data
	      12:   1 at the start of each beat (beginning of diastole)
	      13:   1 during diastole
              14-15: reserved
 * This uses twice as many bits as necessary for the number of states,
 * but it is the best way to do it directly from pressures alone.
 */

/*PUBLIC*/
unsigned int UpdatePeriod(x, flagPtr)
double x[];
unsigned int *flagPtr;
{
#define INPERIOD_POS 0x8
#define PREV_POS 0x9
  unsigned int start = *flagPtr & 0x0F;
  unsigned int end = *flagPtr>>4 & 0x0F;
  unsigned int wasin = *flagPtr & 1<<INPERIOD_POS;
  unsigned int current = (*flagPtr & 3<<PREV_POS) >> (PREV_POS-2);
  
  initvalves(x);
  if (DLI && !(*flagPtr&0x2000)) *flagPtr |= 0x1000; /* Set start cycle flag */
  else *flagPtr &= 0xEFFF; /* Clear start cycle flag */
  if (DLI) *flagPtr |= 0x2000; /* Set (old ) diastole flag */
  else *flagPtr &= ~0x2000;

  if (!start) return TRUE; /* No starting point means all the time */
  current |= DLI | DLO<<1;
  *flagPtr = *flagPtr&0xF9FF | (current&0x3)<<PREV_POS; /* Update flags */
  if (current == start)
    *flagPtr |= 1<<INPERIOD_POS; /* Set the INPERIOD flag */
  else if (current == end)
    *flagPtr &= ~((unsigned int)(1<<INPERIOD_POS)); /* Clear INPERIOD flag */
  if (!end) /* Only return TRUE on up edge */
    return (!wasin && ((*flagPtr)>>INPERIOD_POS & 1)? TRUE : FALSE); 
  else
    return (1 & *flagPtr>>INPERIOD_POS ? TRUE : FALSE); 
}

/*PUBLIC*/
double getval(x, valplace, valtype) /*only returns error (huge number)
				      if totally unreasonable */
     double x[];
     int valplace, valtype;
{
  switch (valtype) {
  case (TIME): return (xt);
  case (PRESSURE): 
     switch (valplace) {
     case (SYSTEM): return(pth);
     case (PC): return (xpPv);
     case (SC): return (xpv);
     default:   return (x[valplace]);
     }
  case (VOLUME): switch (valplace) {
     case (LV): return (Vl0+(xpl-pth)*xCl);
     case (SA): return (Va0+(xpa-A_TH_FRAC*pth)*Ca);
     case (SV): return (sigV+xpv*Cv);
     case (RV): return (Vr0+(xpr-pth)*xCr);
     case (PA): return (VPa0+(xpPa-pth)*CPa);
     case (PV): return (VPv0+(xpPv-pth)*CPv);
     case (PC): return (0.0);
     case (SC): return (0.0);
     case (SYSTEM): return(bv);
     default: return(0.0);
     }
  case (FLOW): switch (valplace) {
     case (LV): 
     case (LVl):
     case (SA): initvalves(x); return (DLO*(xpl-xpa)/Rlo);
     case (SCl):
     case (SC): return ((xpa-xpv)/sigR);
     case (SVl):
     case (SV): initvalves(x); return (DRI*(xpv-xpr)/Rv);
     case (RVl):
     case (RV): 
     case (PA): initvalves(x); return (DRO*(xpr-xpPa)/Rro);
     case (PCl):
     case (PC): return ((xpPa-xpPv)/RPv);
     case (PVl):
     case (PV): initvalves(x); return (DLI*(xpPv-xpl)/Rli);
     case (SYSTEM): return((xpa-xpv)/sigR);
     default: break;
     }
  case (CAPACITANCE): switch(valplace) {
     case (LV): return (xCl);
     case (SA): return(Ca);
     case (SC): return(0.0);
     case (SV): return(Cv);
     case (RV): return (xCr);
     case (PA): return(CPa);
     case (PC): return(0.0);
     case (PV): return(CPv);
     case (SYSTEM): break;
     default: break;
     }   
  case (CAPACITANCE_SYS): switch(valplace) {
     case (LV): return(sigCl);
     case (RV): return(sigCr);
     default: return(getval(x, valplace, CAPACITANCE));
     }
  case (CAPACITANCE_DIAS): switch(valplace) {
     case (LV): return(Cldias);
     case (RV): return(Crdias);
     case (SV): initvalves(x); return(DRI ? xpr : xpv); /* ||| HORRIBLE HACK! Using this Capacitance_dias code for ATRIAL (Inflow ) PRESSURE! */
     case (PV): initvalves(x); return(DLI ? xpl : xpPv);
     default: return(getval(x, valplace, CAPACITANCE));
     }
  case (ZPVOLUME): switch(valplace) {
     case (LV): return(Vl0);
     case (SA): return(Va0);
     case (SC): return(0.0);
     case (SV): return(sigV);
     case (RV): return(Vr0);
     case (PA): return(VPa0);
     case (PC): return(0.0);
     case (PV): return(VPv0);
     case (SYSTEM): return(bv);
     default: return(0.0);
     }
  case (RESISTANCE): switch(valplace) {
     case (LV): return(Rlo);
     case (SA): return(0.0);
     case (SC): return(sigR);
     case (SV): return(Rv);
     case (RV): return(Rro);
     case (PA): return(0.0);
     case (PC): return(RPv);
     case (PV): return(Rli);
     case (SYSTEM): break;
     default: return(0.0);
     }
  case (HEARTRATE): return(sigHR);
  }
  return((double) 0x0ffffffff);
}

static double SetPlaceTypeValue(double x[], int valplace, int valtype, double value);

/*PUBLIC*/
double setval(x, valplace, valtype, value) /* returns value set or 0xffffffff */
     double x[];
     int valtype, valplace;
     double value;
{
  double SetPlaceTypeValue();
  double ret;
/*  if (!var_in_bounds(valtype, valplace, value)) doesn't work yet!
*    return((double)0x0ffffffff);
*/

  SimFlush(); /* in Event.c, to clear data queue of old simulated data. */
  ret = SetPlaceTypeValue(x, valplace, valtype, value);
  UpdateParamWindow();
  return (ret);
}

static double SetPlaceTypeValue(double x[], int valplace, int valtype, double value)
{
  double tmp;
  switch (valtype) {
  case (TIME):
  case (FLOW):
    break;
  case (PRESSURE): switch(valplace) {
    case (SYSTEM): 
      xpl += value - pth;
      xpr += value - pth;
      xpPa += value - pth;
      xpPv += value - pth;
      /* ||| Pa has 1/3 transthoracic pressure on it! */
      xpa += A_TH_FRAC*(value - pth);
      pth = value;
      return(pth);
    default:
      break;
    }
    break;
  case (VOLUME): switch (valplace) {
    case (SYSTEM): 
      bv = value;
      fixvolume(x);
      return(bv);
    default:
      break;
    }
    break;
  case (CAPACITANCE): switch(valplace) {
     case (SA): 
       xpa *= Ca/value;
       return(Ca = value);
     case (SV): 
       xpv *= Cv/value;
       return(Cv = value);
     case (PA): 
       xpPa *= CPa/value;
       return(CPa = value);
     case (PV): 
       xpPv *= CPv/value;
       return(CPv = value);
     case (LV):
     case (RV):
     case (SC):
     case (PC):
     case (SYSTEM): 
       break;
     }   
    break;
   case (CAPACITANCE_SYS): switch(valplace) {
     case (LV): 
       Clsys = value; calcReflexes(); btconst(); 
       tmp = Cl((int) x[9]);
       xpl *= xCl/tmp; xCl = tmp; return(Clsys);
     case (RV): 
       Crsys = value; calcReflexes(); btconst(); 
       tmp = Cr((int) x[9]);
       xpr *= xCr/tmp; xCr = tmp; return(Crsys);
     default: break;
     }
    break;
   case (CAPACITANCE_DIAS): switch(valplace) {
     case (LV): 
       Cldias = value; calcReflexes(); btconst(); 
       tmp = Cl((int) x[9]);
       xpl *= xCl/tmp; xCl = tmp; return(Cldias);
     case (RV): 
       Crdias = value; calcReflexes(); btconst(); 
       tmp = Cr((int) x[9]);
       xpr *= xCr/tmp; xCr = tmp; return(Crdias);
     default: break;
     }
    break;
   case (ZPVOLUME): switch(valplace) {
     case (LV): 
       xpl += (Vl0 - value)/xCl;
       return(Vl0 = value);
     case (SA): 
       xpa += (Va0 - value)/Ca;
       return(Va0 = value);
     case (SV): 
       tmp = sigV; Vv0 = value; calcReflexes();
       xpv += (tmp - sigV)/Cv;
       return(Vv0 = value);
     case (RV): 
       xpr += (Vr0 - value)/xCr;
       return(Vr0 = value);
     case (PA): 
       xpPa = (VPa0 - value)/CPa;
       return(VPa0= value);
     case (PV): 
       xpPv = (VPv0 - value)/CPv;
       return(VPv0 = value);
     case (SC): 
     case (PC):
     case (SYSTEM):
       break;
     }
    break;
   case (RESISTANCE): switch(valplace) {
     case (LV): return(Rlo = value);
     case (SC): 
       Ra = value; 
       calcReflexes();
       return(Ra);
     case (SV): return(Rv = value);
     case (RV): return(Rro = value);
     case (PC): return(RPv = value);
     case (PV): return(Rli = value);
     case (SA):
     case (PA):
     case (SYSTEM): 
       break;
     }
    break;
   case (HEARTRATE): hr = value; return(hr); 
  }
  return((double) 0x0ffffffff); /* error: bad value or location */
  }