/*
   There are a total of 15 entries in the algebraic variable array.
   There are a total of 3 entries in each of the rate and state variable arrays.
   There are a total of 21 entries in the constant variable array.
 */
/*
 * VOI is t in component environment (second).
 * ALGEBRAIC[7] is Pi in component TempCDa (UnitP).
 * ALGEBRAIC[11] is Pi in component TempR (UnitP).
 * ALGEBRAIC[10] is Qo in component TempRC (UnitQ).
 * ALGEBRAIC[12] is Qo in component TempCDv (UnitQ).
 * ALGEBRAIC[3] is Pi in component TempCDv (UnitP).
 * ALGEBRAIC[9] is Qo in component TempCDa (UnitQ).
 * CONSTANTS[0] is CVao in component ParaLeftHeart (UnitCV).
 * ALGEBRAIC[2] is E in component EVentricle (UnitE).
 * STATES[0] is V in component TempCDv (UnitV).
 * CONSTANTS[1] is PlvIni in component ParaLeftHeart (UnitP).
 * CONSTANTS[2] is VlvIni in component ParaLeftHeart (UnitV).
 * ALGEBRAIC[13] is Tao in component TempCDv (dimensionless).
 * CONSTANTS[3] is Vlv0 in component ParaLeftHeart (UnitV).
 * CONSTANTS[4] is CVmi in component ParaLeftHeart (UnitCV).
 * ALGEBRAIC[6] is E in component EAtrium (UnitE).
 * STATES[1] is V in component TempCDa (UnitV).
 * CONSTANTS[5] is PlaIni in component ParaLeftHeart (UnitP).
 * CONSTANTS[6] is VlaIni in component ParaLeftHeart (UnitV).
 * ALGEBRAIC[8] is Tao in component TempCDa (dimensionless).
 * CONSTANTS[7] is Vla0 in component ParaLeftHeart (UnitV).
 * CONSTANTS[8] is ElvMax in component ParaLeftHeart (UnitE).
 * CONSTANTS[9] is ElvMin in component ParaLeftHeart (UnitE).
 * CONSTANTS[10] is T in component ParaLeftHeart (second).
 * CONSTANTS[11] is Ts1 in component ParaLeftHeart (dimensionless).
 * CONSTANTS[12] is Ts2 in component ParaLeftHeart (dimensionless).
 * ALGEBRAIC[0] is mt in component EVentricle (second).
 * ALGEBRAIC[1] is et in component EVentricle (dimensionless).
 * CONSTANTS[13] is ElaMax in component ParaLeftHeart (UnitE).
 * CONSTANTS[14] is ElaMin in component ParaLeftHeart (UnitE).
 * CONSTANTS[15] is Tpwb in component ParaLeftHeart (dimensionless).
 * CONSTANTS[16] is Tpww in component ParaLeftHeart (dimensionless).
 * ALGEBRAIC[4] is mt in component EAtrium (second).
 * ALGEBRAIC[5] is et in component EAtrium (dimensionless).
 * STATES[2] is Pi in component TempRC (UnitP).
 * ALGEBRAIC[14] is Qo in component TempR (UnitQ).
 * CONSTANTS[17] is Rc in component ParaWestkessel (UnitR).
 * CONSTANTS[18] is Rv in component ParaWestkessel (UnitR).
 * CONSTANTS[19] is Cv in component ParaWestkessel (UnitC).
 * CONSTANTS[20] is P0v in component ParaWestkessel (UnitP).
 * RATES[0] is d/dt V in component TempCDv (UnitV).
 * RATES[1] is d/dt V in component TempCDa (UnitV).
 * RATES[2] is d/dt Pi in component TempRC (UnitP).
 */
void
initConsts(double* CONSTANTS, double* RATES, double *STATES)
{
CONSTANTS[0] = 350.;
CONSTANTS[1] = 1.0;
CONSTANTS[2] = 5.0;
CONSTANTS[3] = 500;
CONSTANTS[4] = 400.;
CONSTANTS[5] = 1.0;
CONSTANTS[6] = 4.0;
CONSTANTS[7] = 20;
CONSTANTS[8] = 2.5;
CONSTANTS[9] = 0.1;
CONSTANTS[10] = 1.0;
CONSTANTS[11] = 0.3;
CONSTANTS[12] = 0.45;
CONSTANTS[13] = 0.25;
CONSTANTS[14] = 0.15;
CONSTANTS[15] = 0.92;
CONSTANTS[16] = 0.09;
CONSTANTS[17] = 0.005;
CONSTANTS[18] = 0.65;
CONSTANTS[19] = 2.6;
CONSTANTS[20] = 0.;
STATES[0] = CONSTANTS[3];
STATES[1] = CONSTANTS[7];
STATES[2] = CONSTANTS[20];
}
void
computeRates(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[4] = VOI -  CONSTANTS[10]*floor(VOI/CONSTANTS[10]);
ALGEBRAIC[5] = (ALGEBRAIC[4]>=0.00000&&ALGEBRAIC[4]<= ((CONSTANTS[15]+CONSTANTS[16]) - 1.00000)*CONSTANTS[10] ? 1.00000 - cos(( 2.00000*3.14159*(ALGEBRAIC[4] -  (CONSTANTS[15] - 1.00000)*CONSTANTS[10]))/( CONSTANTS[16]*CONSTANTS[10])) : ALGEBRAIC[4]> ((CONSTANTS[15]+CONSTANTS[16]) - 1.00000)*CONSTANTS[10]&&ALGEBRAIC[4]<= CONSTANTS[15]*CONSTANTS[10] ? 0.00000 : ALGEBRAIC[4]> CONSTANTS[15]*CONSTANTS[10]&&ALGEBRAIC[4]<=CONSTANTS[10] ? 1.00000 - cos(( 2.00000*3.14159*(ALGEBRAIC[4] -  CONSTANTS[15]*CONSTANTS[10]))/( CONSTANTS[16]*CONSTANTS[10])) : 0.0/0.0);
ALGEBRAIC[6] = CONSTANTS[14]+( ALGEBRAIC[5]*(CONSTANTS[13] - CONSTANTS[14]))/2.00000;
ALGEBRAIC[7] = CONSTANTS[5]+ ALGEBRAIC[6]*(STATES[1] - CONSTANTS[6]);
ALGEBRAIC[10] = (STATES[2] - ALGEBRAIC[7])/CONSTANTS[18];
ALGEBRAIC[0] = VOI -  CONSTANTS[10]*floor(VOI/CONSTANTS[10]);
ALGEBRAIC[1] = (ALGEBRAIC[0]>=0.00000&&ALGEBRAIC[0]<= CONSTANTS[11]*CONSTANTS[10] ? 1.00000 - cos(( 3.14159*ALGEBRAIC[0])/( CONSTANTS[11]*CONSTANTS[10])) : ALGEBRAIC[0]> CONSTANTS[11]*CONSTANTS[10]&&ALGEBRAIC[0]<= CONSTANTS[12]*CONSTANTS[10] ? 1.00000+cos(( 3.14159*(ALGEBRAIC[0] -  CONSTANTS[11]*CONSTANTS[10]))/( (CONSTANTS[12] - CONSTANTS[11])*CONSTANTS[10])) : ALGEBRAIC[0]> CONSTANTS[12]*CONSTANTS[10]&&ALGEBRAIC[0]<CONSTANTS[10] ? 0.00000 : 0.0/0.0);
ALGEBRAIC[2] = CONSTANTS[9]+( ALGEBRAIC[1]*(CONSTANTS[8] - CONSTANTS[9]))/2.00000;
ALGEBRAIC[3] = CONSTANTS[1]+ ALGEBRAIC[2]*(STATES[0] - CONSTANTS[2]);
ALGEBRAIC[8] = (ALGEBRAIC[7]>=ALGEBRAIC[3] ? 1.00000 : ALGEBRAIC[7]<ALGEBRAIC[3] ? 0.00000 : 0.0/0.0);
ALGEBRAIC[9] = (ALGEBRAIC[7]>=ALGEBRAIC[3] ?  CONSTANTS[4]*ALGEBRAIC[8]*pow(fabs(ALGEBRAIC[7] - ALGEBRAIC[3]), 0.500000) : ALGEBRAIC[7]<ALGEBRAIC[3] ?  -1.00000*CONSTANTS[4]*ALGEBRAIC[8]*pow(fabs(ALGEBRAIC[3] - ALGEBRAIC[7]), 0.500000) : 0.0/0.0);
RATES[1] = ALGEBRAIC[10] - ALGEBRAIC[9];
rootfind_0(VOI, CONSTANTS, RATES, STATES, ALGEBRAIC, pret);
RATES[0] = ALGEBRAIC[9] - ALGEBRAIC[12];
ALGEBRAIC[14] = ALGEBRAIC[12];
RATES[2] = (ALGEBRAIC[14] - ALGEBRAIC[10])/CONSTANTS[19];
}
void
computeVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[4] = VOI -  CONSTANTS[10]*floor(VOI/CONSTANTS[10]);
ALGEBRAIC[5] = (ALGEBRAIC[4]>=0.00000&&ALGEBRAIC[4]<= ((CONSTANTS[15]+CONSTANTS[16]) - 1.00000)*CONSTANTS[10] ? 1.00000 - cos(( 2.00000*3.14159*(ALGEBRAIC[4] -  (CONSTANTS[15] - 1.00000)*CONSTANTS[10]))/( CONSTANTS[16]*CONSTANTS[10])) : ALGEBRAIC[4]> ((CONSTANTS[15]+CONSTANTS[16]) - 1.00000)*CONSTANTS[10]&&ALGEBRAIC[4]<= CONSTANTS[15]*CONSTANTS[10] ? 0.00000 : ALGEBRAIC[4]> CONSTANTS[15]*CONSTANTS[10]&&ALGEBRAIC[4]<=CONSTANTS[10] ? 1.00000 - cos(( 2.00000*3.14159*(ALGEBRAIC[4] -  CONSTANTS[15]*CONSTANTS[10]))/( CONSTANTS[16]*CONSTANTS[10])) : 0.0/0.0);
ALGEBRAIC[6] = CONSTANTS[14]+( ALGEBRAIC[5]*(CONSTANTS[13] - CONSTANTS[14]))/2.00000;
ALGEBRAIC[7] = CONSTANTS[5]+ ALGEBRAIC[6]*(STATES[1] - CONSTANTS[6]);
ALGEBRAIC[10] = (STATES[2] - ALGEBRAIC[7])/CONSTANTS[18];
ALGEBRAIC[0] = VOI -  CONSTANTS[10]*floor(VOI/CONSTANTS[10]);
ALGEBRAIC[1] = (ALGEBRAIC[0]>=0.00000&&ALGEBRAIC[0]<= CONSTANTS[11]*CONSTANTS[10] ? 1.00000 - cos(( 3.14159*ALGEBRAIC[0])/( CONSTANTS[11]*CONSTANTS[10])) : ALGEBRAIC[0]> CONSTANTS[11]*CONSTANTS[10]&&ALGEBRAIC[0]<= CONSTANTS[12]*CONSTANTS[10] ? 1.00000+cos(( 3.14159*(ALGEBRAIC[0] -  CONSTANTS[11]*CONSTANTS[10]))/( (CONSTANTS[12] - CONSTANTS[11])*CONSTANTS[10])) : ALGEBRAIC[0]> CONSTANTS[12]*CONSTANTS[10]&&ALGEBRAIC[0]<CONSTANTS[10] ? 0.00000 : 0.0/0.0);
ALGEBRAIC[2] = CONSTANTS[9]+( ALGEBRAIC[1]*(CONSTANTS[8] - CONSTANTS[9]))/2.00000;
ALGEBRAIC[3] = CONSTANTS[1]+ ALGEBRAIC[2]*(STATES[0] - CONSTANTS[2]);
ALGEBRAIC[8] = (ALGEBRAIC[7]>=ALGEBRAIC[3] ? 1.00000 : ALGEBRAIC[7]<ALGEBRAIC[3] ? 0.00000 : 0.0/0.0);
ALGEBRAIC[9] = (ALGEBRAIC[7]>=ALGEBRAIC[3] ?  CONSTANTS[4]*ALGEBRAIC[8]*pow(fabs(ALGEBRAIC[7] - ALGEBRAIC[3]), 0.500000) : ALGEBRAIC[7]<ALGEBRAIC[3] ?  -1.00000*CONSTANTS[4]*ALGEBRAIC[8]*pow(fabs(ALGEBRAIC[3] - ALGEBRAIC[7]), 0.500000) : 0.0/0.0);
ALGEBRAIC[14] = ALGEBRAIC[12];
}

void objfunc_0(double *p, double *hx, int m, int n, void *adata)
{
  struct rootfind_info* rfi = (struct rootfind_info*)adata;
#define VOI rfi->aVOI
#define CONSTANTS rfi->aCONSTANTS
#define RATES rfi->aRATES
#define STATES rfi->aSTATES
#define ALGEBRAIC rfi->aALGEBRAIC
#define pret rfi->aPRET
  ALGEBRAIC[11] = p[0];
  ALGEBRAIC[12] = p[1];
  ALGEBRAIC[13] = p[2];
  hx[0] = ALGEBRAIC[13] - (ALGEBRAIC[3]>=ALGEBRAIC[11] ? 1.00000 : ALGEBRAIC[3]<ALGEBRAIC[11] ? 0.00000 : 0.0/0.0);
  hx[1] = ALGEBRAIC[12] - (ALGEBRAIC[3]>=ALGEBRAIC[11] ?  CONSTANTS[0]*ALGEBRAIC[13]*pow(fabs(ALGEBRAIC[3] - ALGEBRAIC[11]), 0.500000) : ALGEBRAIC[3]<ALGEBRAIC[11] ?  -1.00000*CONSTANTS[0]*ALGEBRAIC[13]*pow(fabs(ALGEBRAIC[11] - ALGEBRAIC[3]), 0.500000) : 0.0/0.0);
  hx[2] = ALGEBRAIC[11] - (STATES[2]+ CONSTANTS[17]*ALGEBRAIC[12]);
#undef VOI
#undef CONSTANTS
#undef RATES
#undef STATES
#undef ALGEBRAIC
#undef pret
}
void rootfind_0(double VOI, double* CONSTANTS, double* RATES,
double* STATES, double* ALGEBRAIC, int* pret)
{
  static double p[3] = {0.1,0.1,0.1};
  double bp[3], work[LM_DIF_WORKSZ(3, 3)];
  struct rootfind_info rfi;
  rfi.aVOI = VOI;
  rfi.aCONSTANTS = CONSTANTS;
  rfi.aRATES = RATES;
  rfi.aSTATES = STATES;
  rfi.aALGEBRAIC = ALGEBRAIC;
  rfi.aPRET = pret;
  do_levmar(objfunc_0, p, bp, work, pret, 3, &rfi);
  ALGEBRAIC[11] = p[0];
  ALGEBRAIC[12] = p[1];
  ALGEBRAIC[13] = p[2];
}