/*
   There are a total of 12 entries in the algebraic variable array.
   There are a total of 15 entries in each of the rate and state variable arrays.
   There are a total of 49 entries in the constant variable array.
 */
/*
 * VOI is time in component Environment (second).
 * CONSTANTS[0] is Eta_b in component Parameters (mmHg_sec).
 * CONSTANTS[1] is Rho_b in component Parameters (gram_per_cubic_meter).
 * CONSTANTS[2] is E_p in component Parameters (mmHg).
 * CONSTANTS[3] is r_p in component Parameters (meter).
 * CONSTANTS[4] is l_p in component Parameters (meter).
 * CONSTANTS[5] is h_p in component Parameters (meter).
 * CONSTANTS[6] is r_a0 in component Parameters (meter).
 * CONSTANTS[7] is h_a0 in component Parameters (meter).
 * CONSTANTS[8] is sigma_ae0 in component Parameters (mmHg).
 * CONSTANTS[9] is K_a_sigma in component Parameters (dimensionless).
 * CONSTANTS[10] is sigma_ac in component Parameters (mmHg).
 * CONSTANTS[11] is T_a_max0 in component Parameters (mmHg_m).
 * CONSTANTS[12] is r_am in component Parameters (meter).
 * CONSTANTS[13] is r_at in component Parameters (meter).
 * CONSTANTS[14] is n_am in component Parameters (dimensionless).
 * CONSTANTS[15] is Eta_a in component Parameters (mmHg_sec).
 * CONSTANTS[16] is Q_bl in component Parameters (cubic_m_per_sec).
 * CONSTANTS[17] is P_ic in component Parameters (mmHg).
 * CONSTANTS[18] is R_v in component Parameters (mmHg_sec_per_cubic_m).
 * CONSTANTS[19] is C_v in component Parameters (cubic_m_per_mmHg).
 * CONSTANTS[20] is P_ai in component Parameters (mmHg).
 * CONSTANTS[21] is T_myo0 in component Parameters (mmHg_m).
 * CONSTANTS[22] is T_myo_s in component Parameters (mmHg_m).
 * CONSTANTS[23] is Tau_myo in component Parameters (second).
 * CONSTANTS[24] is Tau_shear in component Parameters (second).
 * CONSTANTS[25] is G_myo in component Parameters (dimensionless).
 * CONSTANTS[26] is G_shear in component Parameters (dimensionless).
 * CONSTANTS[27] is G_neuro in component Parameters (dimensionless).
 * CONSTANTS[28] is G_meta in component Parameters (dimensionless).
 * CONSTANTS[29] is x_ini in component Parameters (dimensionless).
 * CONSTANTS[30] is Tau_neuro in component Parameters (second).
 * CONSTANTS[31] is Tau_meta in component Parameters (second).
 * CONSTANTS[32] is conc_CO2_a in component Parameters (mol_per_m3).
 * CONSTANTS[33] is M_CO2_0 in component Parameters (mol_per_sec).
 * CONSTANTS[34] is alpha_tv in component Parameters (dimensionless).
 * CONSTANTS[35] is beta_tv in component Parameters (mol_per_m3).
 * CONSTANTS[36] is vol_PVC in component Parameters (cubic_m).
 * CONSTANTS[37] is K_shear in component Parameters (second).
 * CONSTANTS[38] is K_Ra in component Parameters (mmHg_sec_m).
 * CONSTANTS[39] is conc_CO2_t_sp in component Parameters (mol_per_m3).
 * CONSTANTS[40] is f_Q in component Parameters (dimensionless).
 * CONSTANTS[41] is n_QM in component Parameters (dimensionless).
 * CONSTANTS[42] is K_Va in component Parameters (meter).
 * CONSTANTS[43] is P_vo in component Parameters (mmHg).
 * STATES[0] is Q_Lp in component PCA (cubic_m_per_sec).
 * STATES[1] is P_p in component PCA (mmHg).
 * CONSTANTS[46] is C_p in component PCA_Components (cubic_m_per_mmHg).
 * CONSTANTS[47] is L_p in component PCA_Components (mmHg_sec2_per_m3).
 * CONSTANTS[48] is R_p in component PCA_Components (mmHg_sec_per_cubic_m).
 * ALGEBRAIC[0] is R_a in component AC_Components (mmHg_sec_per_cubic_m).
 * STATES[9] is P_a in component Arteriolar_Circulation (mmHg).
 * ALGEBRAIC[1] is C_a in component AC_Components (cubic_m_per_mmHg).
 * STATES[10] is Q_a in component AC_Flow (cubic_m_per_sec).
 * STATES[2] is P_v in component Venous_Circulation (mmHg).
 * STATES[3] is r_a in component Arteriolar_Circulation (meter).
 * ALGEBRAIC[8] is T_ae in component AC_Tension (mmHg_m).
 * ALGEBRAIC[11] is T_am in component AC_Tension (mmHg_m).
 * STATES[11] is T_av in component AC_Tension (mmHg_m).
 * STATES[12] is T_a in component AC_Tension (mmHg_m).
 * ALGEBRAIC[2] is h_a in component AC_Tension_Variables (meter).
 * ALGEBRAIC[10] is T_a_max in component AC_Tension_Variables (mmHg_m).
 * ALGEBRAIC[9] is M_s in component CBF_Regulation (dimensionless).
 * ALGEBRAIC[3] is M_s1 in component CBF_Regulation (dimensionless).
 * STATES[4] is x_myo in component Myo_Regulation (dimensionless).
 * STATES[5] is x_shear in component Shear_Regulation (dimensionless).
 * STATES[6] is x_neuro in component Neuro_Regulation (dimensionless).
 * STATES[7] is x_meta in component Meta_Regulation (dimensionless).
 * STATES[13] is A_myo in component Myo_Regulation (dimensionless).
 * STATES[14] is A_shear in component Shear_Regulation (dimensionless).
 * ALGEBRAIC[4] is A_neuro in component Neuro_Regulation (dimensionless).
 * ALGEBRAIC[5] is A_Meta in component Meta_Regulation (dimensionless).
 * CONSTANTS[44] is G_a_Meta in component Meta_Regulation (m3_per_mol).
 * STATES[8] is conc_CO2_t in component Meta_CO2_Kinetics (mol_per_m3).
 * ALGEBRAIC[6] is conc_CO2_v in component Meta_CO2_Kinetics (mol_per_m3).
 * ALGEBRAIC[7] is M_CO2 in component VC_Meta_Stim (mol_per_sec).
 * CONSTANTS[45] is unity in component VC_Meta_Stim (dimensionless).
 * RATES[1] is d/dt P_p in component PCA (mmHg).
 * RATES[0] is d/dt Q_Lp in component PCA (cubic_m_per_sec).
 * RATES[3] is d/dt r_a in component Arteriolar_Circulation (meter).
 * RATES[2] is d/dt P_v in component Venous_Circulation (mmHg).
 * RATES[4] is d/dt x_myo in component Myo_Regulation (dimensionless).
 * RATES[5] is d/dt x_shear in component Shear_Regulation (dimensionless).
 * RATES[6] is d/dt x_neuro in component Neuro_Regulation (dimensionless).
 * RATES[7] is d/dt x_meta in component Meta_Regulation (dimensionless).
 * RATES[8] is d/dt conc_CO2_t in component Meta_CO2_Kinetics (mol_per_m3).
 * There are a total of 1 condition variables.
 */
void
initConsts(double* CONSTANTS, double* RATES, double *STATES)
{
CONSTANTS[0] = 3E-5;
CONSTANTS[1] = 1.05E6;
CONSTANTS[2] = 12E3;
CONSTANTS[3] = 1.05E-3;
CONSTANTS[4] = 8.6E-2;
CONSTANTS[5] = 2.6E-4;
CONSTANTS[6] = 7.5E-5;
CONSTANTS[7] = 2.5E-5;
CONSTANTS[8] = 11.19;
CONSTANTS[9] = 4.5;
CONSTANTS[10] = 41.32;
CONSTANTS[11] = 1.50E-2;
CONSTANTS[12] = 1.28E-4;
CONSTANTS[13] = 1.74E-4;
CONSTANTS[14] = 1.75;
CONSTANTS[15] = 47.8;
CONSTANTS[16] = 8.8E-7;
CONSTANTS[17] = 10;
CONSTANTS[18] = 1.6E7;
CONSTANTS[19] = 2.5E-8;
CONSTANTS[20] = 96.0;
CONSTANTS[21] = 3.6E-3;
CONSTANTS[22] = 3E-3;
CONSTANTS[23] = 7;
CONSTANTS[24] = 60;
CONSTANTS[25] = 4;
CONSTANTS[26] = 0;
CONSTANTS[27] = -0.7;
CONSTANTS[28] = 0;
CONSTANTS[29] = 0;
CONSTANTS[30] = 5;
CONSTANTS[31] = 15;
CONSTANTS[32] = 20.65;
CONSTANTS[33] = 1.35E-7;
CONSTANTS[34] = 0.96;
CONSTANTS[35] = 8.9;
CONSTANTS[36] = 5.99E-6;
CONSTANTS[37] = 5.2E-7;
CONSTANTS[38] = 2.5E-9;
CONSTANTS[39] = 12.41;
CONSTANTS[40] = 1.25;
CONSTANTS[41] = 2.2;
CONSTANTS[42] = 120;
CONSTANTS[43] = 14;
STATES[0] = 8.8E-7;
STATES[1] = 91;
STATES[2] = 21;
STATES[3] = 7.7E-5;
STATES[4] = 0;
STATES[5] = 0;
STATES[6] = 0;
STATES[7] = 0;
CONSTANTS[44] = 0.59;
STATES[8] = 12.41;
CONSTANTS[45] = 1;
CONSTANTS[46] =  (( 3.00000* 3.14159265358979*pow(CONSTANTS[3], 2.00000)*pow(CONSTANTS[3]/CONSTANTS[5]+1.00000, 2.00000))/( CONSTANTS[2]*(( 2.00000*CONSTANTS[3])/CONSTANTS[5]+1.00000)))*CONSTANTS[4];
CONSTANTS[47] =  (CONSTANTS[1]/(  3.14159265358979*pow(CONSTANTS[3], 2.00000)))*CONSTANTS[4]*7.50000e-06;
CONSTANTS[48] =  (( 8.00000*CONSTANTS[0])/(  3.14159265358979*pow(CONSTANTS[3], 4.00000)))*CONSTANTS[4];
STATES[9] = 0.1001;
STATES[10] = 0.1001;
STATES[11] = 0.1001;
STATES[12] = 0.1001;
STATES[13] = 0.1001;
STATES[14] = 0.1001;
RATES[1] = 0.1001;
RATES[0] = 0.1001;
RATES[3] = 0.1001;
RATES[2] = 0.1001;
RATES[4] = 0.1001;
RATES[5] = 0.1001;
RATES[6] = 0.1001;
RATES[7] = 0.1001;
RATES[8] = 0.1001;
}
void
computeResiduals(double VOI, double* CONSTANTS, double* RATES, double* OLDRATES, double* STATES,
                 double* OLDSTATES, double* ALGEBRAIC, double* CONDVARS)
{
resid[0] = RATES[1] - ( 2.00000*(( STATES[0]*ALGEBRAIC[0] -  2.00000*STATES[1])+ 2.00000*STATES[9]))/( ALGEBRAIC[0]*CONSTANTS[46]);
resid[1] = RATES[0] - ((CONSTANTS[20] - STATES[1]) -  STATES[0]*CONSTANTS[48])/CONSTANTS[47];
resid[2] = RATES[3] - ( CONSTANTS[6]*(( STATES[1]*(ALGEBRAIC[0]+CONSTANTS[18])*STATES[3]+ STATES[2]*ALGEBRAIC[0]*STATES[3]) -  (ALGEBRAIC[8]+ALGEBRAIC[11]+ CONSTANTS[17]*(STATES[3]+ALGEBRAIC[2]))*( 2.00000*ALGEBRAIC[0]+CONSTANTS[18])))/( CONSTANTS[6]*pow(STATES[3], 2.00000)*CONSTANTS[42]*ALGEBRAIC[0]*(ALGEBRAIC[0]+CONSTANTS[18])+ CONSTANTS[15]*ALGEBRAIC[2]*( 2.00000*ALGEBRAIC[0]+CONSTANTS[18]));
resid[3] = STATES[9] - (STATES[12]+ CONSTANTS[17]*(STATES[3]+ALGEBRAIC[2]))/STATES[3];
resid[4] = STATES[10] - ( 2.00000*(STATES[9] - STATES[2]))/(ALGEBRAIC[0]+CONSTANTS[18]);
resid[5] = RATES[2] -  (2.00000/CONSTANTS[19])*((STATES[9] - STATES[2])/(ALGEBRAIC[0]+CONSTANTS[18]) - (STATES[2] - CONSTANTS[43])/CONSTANTS[18]);
resid[6] = STATES[12] - ALGEBRAIC[8]+ALGEBRAIC[11]+STATES[11];
resid[7] = STATES[11] -  (CONSTANTS[15]/CONSTANTS[6])*RATES[3]*ALGEBRAIC[2];
resid[8] = STATES[13] - (STATES[12] - CONSTANTS[21])/CONSTANTS[22];
resid[9] = RATES[4] - (STATES[13] - STATES[4])/CONSTANTS[23];
resid[10] = STATES[14] - ( CONSTANTS[37]*STATES[10])/pow(STATES[3], 3.00000) - 1.00000;
resid[11] = RATES[5] - (STATES[14] - STATES[5])/CONSTANTS[24];
resid[12] = RATES[6] - (ALGEBRAIC[4] - STATES[6])/CONSTANTS[30];
resid[13] = RATES[7] - (ALGEBRAIC[5] - STATES[7])/CONSTANTS[31];
resid[14] = RATES[8] -  (1.00000/CONSTANTS[36])*(ALGEBRAIC[7] -  STATES[10]*(ALGEBRAIC[6] - CONSTANTS[32]));
}
void
computeVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[1] = 1.03120/ALGEBRAIC[0];
}
void
computeEssentialVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[0] = CONSTANTS[38]/pow(STATES[3], 4.00000);
ALGEBRAIC[2] =  pow((pow(STATES[3], 2.00000)+ 2.00000*CONSTANTS[6]*CONSTANTS[7]+pow(CONSTANTS[7], 2.00000)), 1.0 / 2) - STATES[3];
ALGEBRAIC[4] = (CONDVAR[0]>=0.00000 ? 1.00000 : 0.00000);
ALGEBRAIC[5] =  CONSTANTS[44]*(STATES[8] - CONSTANTS[39]);
ALGEBRAIC[6] =  CONSTANTS[34]*STATES[8]+CONSTANTS[35];
ALGEBRAIC[7] =  CONSTANTS[33]*(CONSTANTS[45]+ ALGEBRAIC[4]*((CONSTANTS[40] - CONSTANTS[45])/CONSTANTS[41]));
ALGEBRAIC[8] =  ALGEBRAIC[2]*( CONSTANTS[8]*(exp(( CONSTANTS[9]*(STATES[3] - CONSTANTS[6]))/CONSTANTS[6]) - 1.00000) - CONSTANTS[10]);
ALGEBRAIC[3] =  CONSTANTS[25]*STATES[4]+ CONSTANTS[26]*STATES[5]+ CONSTANTS[27]*STATES[6]+ CONSTANTS[28]*STATES[7]+CONSTANTS[29];
ALGEBRAIC[9] = (exp( 2.00000*ALGEBRAIC[3]) - 1.00000)/(exp( 2.00000*ALGEBRAIC[3])+1.00000);
ALGEBRAIC[10] =  CONSTANTS[11]*(1.00000+ALGEBRAIC[9]);
ALGEBRAIC[11] =  ALGEBRAIC[10]*exp(- pow(fabs((STATES[3] - CONSTANTS[12])/(CONSTANTS[13] - CONSTANTS[12])), CONSTANTS[14]));
}
void
getStateInformation(double* SI)
{
SI[0] = 1.0;
SI[1] = 1.0;
SI[9] = 0.0;
SI[10] = 0.0;
SI[2] = 1.0;
SI[3] = 1.0;
SI[11] = 0.0;
SI[12] = 0.0;
SI[4] = 1.0;
SI[5] = 1.0;
SI[6] = 1.0;
SI[7] = 1.0;
SI[13] = 0.0;
SI[14] = 0.0;
SI[8] = 1.0;
}
void
computeRoots(double VOI, double* CONSTANTS, double* RATES, double* OLDRATES, double* STATES,
             double* OLDSTATES, double* ALGEBRAIC, double* CONDVARS)
{
CONDVAR[0] = VOI - 1.00000;
}