/*
   There are a total of 26 entries in the algebraic variable array.
   There are a total of 8 entries in each of the rate and state variable arrays.
   There are a total of 46 entries in the constant variable array.
 */
/*
 * VOI is time in component environment (second).
 * CONSTANTS[0] is volume_i0 in component volume_i (volume).
 * CONSTANTS[43] is volume_i in component volume_i (volume).
 * CONSTANTS[1] is gamma in component volume_er (dimensionless).
 * CONSTANTS[44] is volume_er in component volume_er (volume).
 * STATES[0] is volumeCa_i in component volumeCa_i (microMVolume).
 * ALGEBRAIC[10] is Jipt in component Jipt (microMvolumepersec).
 * ALGEBRAIC[13] is Jryr in component Jryr (microMvolumepersec).
 * ALGEBRAIC[14] is Jer in component Jer (microMvolumepersec).
 * ALGEBRAIC[15] is Jserca in component Jserca (microMvolumepersec).
 * ALGEBRAIC[21] is Jin in component Jin (microMvolumepersec).
 * ALGEBRAIC[22] is Jpm in component Jpm (microMvolumepersec).
 * ALGEBRAIC[0] is Ca_i in component Ca_i (microM).
 * STATES[1] is volumeCa_er in component volumeCa_er (microMVolume).
 * ALGEBRAIC[1] is Ca_er in component Ca_er (microM).
 * STATES[2] is volume_iIP3 in component volume_iIP3 (microMVolume).
 * ALGEBRAIC[6] is J_ip3P in component J_ip3P (microMpersec).
 * ALGEBRAIC[8] is J_ip3D in component J_ip3D (microMpersec).
 * ALGEBRAIC[2] is IP3 in component IP3 (microM).
 * STATES[3] is PIP2 in component PIP2 (dimensionless).
 * CONSTANTS[2] is PIP2_Total in component PIP2 (microM).
 * CONSTANTS[45] is J_ip3R in component J_ip3R (per_sec).
 * ALGEBRAIC[3] is IPX in component IPX (dimensionless).
 * CONSTANTS[3] is PIP2_Total in component IPX (microM).
 * CONSTANTS[4] is J_ip3R0 in component J_ip3R (per_sec).
 * STATES[4] is ATP_e in component receptor (microM).
 * ALGEBRAIC[5] is V_IP3 in component receptor (microMpersec).
 * CONSTANTS[5] is MechanicalStimulation in component receptor (dimensionless).
 * ALGEBRAIC[4] is V_mech in component receptor (dimensionless).
 * CONSTANTS[6] is V_mech0 in component receptor (dimensionless).
 * CONSTANTS[7] is C_1 in component receptor (microMpersec).
 * CONSTANTS[8] is C_2 in component receptor (microM).
 * CONSTANTS[9] is C_3 in component receptor (microMpersec).
 * CONSTANTS[10] is C_4 in component receptor (microM).
 * CONSTANTS[11] is V_ATP in component receptor (microMpersec).
 * CONSTANTS[12] is K_ATP in component receptor (microM).
 * CONSTANTS[13] is Beta_1 in component J_ip3D (microMpersec).
 * CONSTANTS[14] is Beta_2 in component J_ip3D (microMpersec).
 * ALGEBRAIC[7] is F_Ca in component J_ip3D (dimensionless).
 * CONSTANTS[15] is K_rc in component J_ip3D (microM).
 * CONSTANTS[16] is K_IP3 in component J_ip3D (microM).
 * ALGEBRAIC[9] is Pipt in component Jipt (dimensionless).
 * STATES[5] is O in component O (dimensionless).
 * CONSTANTS[17] is Jipt0 in component Jipt (per_sec).
 * CONSTANTS[18] is Ka in component Jryr (microM).
 * CONSTANTS[19] is Kb in component Jryr (microM).
 * CONSTANTS[20] is Kc in component Jryr (dimensionless).
 * CONSTANTS[21] is Jryr0 in component Jryr (per_sec).
 * ALGEBRAIC[12] is Pryr in component Jryr (dimensionless).
 * ALGEBRAIC[11] is W_inf in component Jryr (dimensionless).
 * CONSTANTS[22] is Jer0 in component Jer (per_sec).
 * CONSTANTS[23] is Vserca in component Jserca (microM2persec).
 * CONSTANTS[24] is Kserca in component Jserca (microM).
 * CONSTANTS[25] is Vm in component Jin (milliV).
 * ALGEBRAIC[16] is V_Ca in component Jin (milliV).
 * CONSTANTS[26] is R in component Jin (millijoule_per_mole_kelvin).
 * CONSTANTS[27] is T in component Jin (kelvin).
 * CONSTANTS[28] is F in component Jin (coulomb_per_mole).
 * CONSTANTS[29] is z_Ca in component Jin (dimensionless).
 * CONSTANTS[30] is Ca_ext in component Jin (microM).
 * CONSTANTS[31] is P_pm_1 in component Jin (conductancepervolume).
 * CONSTANTS[32] is P_pm_2 in component Jin (conductancepervolume).
 * ALGEBRAIC[17] is Iin_1 in component Jin (current).
 * ALGEBRAIC[18] is Iin_2 in component Jin (current).
 * ALGEBRAIC[19] is Jin_1 in component Jin (microMvolumepersec).
 * ALGEBRAIC[20] is Jin_2 in component Jin (microMvolumepersec).
 * CONSTANTS[33] is Vpm in component Jpm (microMpersec).
 * CONSTANTS[34] is Kpm in component Jpm (microM).
 * ALGEBRAIC[23] is S in component O (dimensionless).
 * STATES[6] is I1 in component I1 (dimensionless).
 * STATES[7] is I2 in component I2 (dimensionless).
 * ALGEBRAIC[24] is k1 in component constants (per_secmicroM).
 * CONSTANTS[35] is k1_ in component constants (first_order_rate_constant).
 * CONSTANTS[36] is k2 in component constants (first_order_rate_constant).
 * CONSTANTS[37] is k3 in component constants (first_order_rate_constant).
 * ALGEBRAIC[25] is k4 in component constants (first_order_rate_constant).
 * CONSTANTS[38] is k5 in component constants (first_order_rate_constant).
 * CONSTANTS[39] is alpha_1 in component constants (per_secmicroM).
 * CONSTANTS[40] is beta_1 in component constants (microM).
 * CONSTANTS[41] is alpha_4 in component constants (first_order_rate_constant).
 * CONSTANTS[42] is beta_4 in component constants (microM).
 * RATES[0] is d/dt volumeCa_i in component volumeCa_i (microMVolume).
 * RATES[1] is d/dt volumeCa_er in component volumeCa_er (microMVolume).
 * RATES[2] is d/dt volume_iIP3 in component volume_iIP3 (microMVolume).
 * RATES[3] is d/dt PIP2 in component PIP2 (dimensionless).
 * RATES[4] is d/dt ATP_e in component receptor (microM).
 * RATES[5] is d/dt O in component O (dimensionless).
 * RATES[6] is d/dt I1 in component I1 (dimensionless).
 * RATES[7] is d/dt I2 in component I2 (dimensionless).
 * There are a total of 3 condition variables.
 */
void
initConsts(double* CONSTANTS, double* RATES, double *STATES)
{
CONSTANTS[0] = 1450;
CONSTANTS[1] = 4.9;
STATES[0] = 75.8335890358239;
STATES[1] = 1199.73747802144;
STATES[2] = -3.28287434380887e-93;
STATES[3] = 1;
CONSTANTS[2] = 5;
CONSTANTS[3] = 25;
CONSTANTS[4] = 10;
STATES[4] = 0;
CONSTANTS[5] = 0;
CONSTANTS[6] = 1;
CONSTANTS[7] = 2.6;
CONSTANTS[8] = 0.02;
CONSTANTS[9] = 2.3;
CONSTANTS[10] = 10;
CONSTANTS[11] = 0.001;
CONSTANTS[12] = 0.0001;
CONSTANTS[13] = 0;
CONSTANTS[14] = 36;
CONSTANTS[15] = 0.3;
CONSTANTS[16] = 0.01;
STATES[5] = -2.42013500345263e-89;
CONSTANTS[17] = 28;
CONSTANTS[18] = 0.3722419;
CONSTANTS[19] = 0.636005;
CONSTANTS[20] = 0.0571428;
CONSTANTS[21] = 0;
CONSTANTS[22] = 0.0035;
CONSTANTS[23] = 0.09;
CONSTANTS[24] = 0.04;
CONSTANTS[25] = -39;
CONSTANTS[26] = 8314;
CONSTANTS[27] = 310;
CONSTANTS[28] = 96485;
CONSTANTS[29] = 2;
CONSTANTS[30] = 1300;
CONSTANTS[31] = 0.6;
CONSTANTS[32] = 28;
CONSTANTS[33] = 0.072;
CONSTANTS[34] = 0.6;
STATES[6] = 2.24343309680416e-81;
STATES[7] = 1.01276768045278e-40;
CONSTANTS[35] = 0.88;
CONSTANTS[36] = 0.5;
CONSTANTS[37] = 0.5;
CONSTANTS[38] = 0.02;
CONSTANTS[39] = 40;
CONSTANTS[40] = 0.8;
CONSTANTS[41] = 0.06;
CONSTANTS[42] = 0.01;
CONSTANTS[43] = CONSTANTS[0];
CONSTANTS[44] = CONSTANTS[43]/CONSTANTS[1];
CONSTANTS[45] = CONSTANTS[4];
RATES[0] = 0.1001;
RATES[1] = 0.1001;
RATES[2] = 0.1001;
RATES[3] = 0.1001;
RATES[4] = 0.1001;
RATES[5] = 0.1001;
RATES[6] = 0.1001;
RATES[7] = 0.1001;
}
void
computeResiduals(double VOI, double* CONSTANTS, double* RATES, double* OLDRATES, double* STATES,
                 double* OLDSTATES, double* ALGEBRAIC, double* CONDVARS)
{
resid[0] = RATES[0] - ((ALGEBRAIC[10]+ALGEBRAIC[13]+ALGEBRAIC[14]) - ALGEBRAIC[15])+(ALGEBRAIC[21] - ALGEBRAIC[22]);
resid[1] = RATES[1] - - ((ALGEBRAIC[10]+ALGEBRAIC[13]+ALGEBRAIC[14]) - ALGEBRAIC[15]);
resid[2] = RATES[2] -  (ALGEBRAIC[6] - ALGEBRAIC[8])*CONSTANTS[43];
resid[3] = RATES[3] -  CONSTANTS[45]*ALGEBRAIC[3] - ALGEBRAIC[6]/CONSTANTS[2];
resid[4] = RATES[4] - ( - CONSTANTS[11]*STATES[4])/(CONSTANTS[12]+STATES[4]);
resid[5] = RATES[5] -  ALGEBRAIC[24]*ALGEBRAIC[2]*ALGEBRAIC[23] - ( CONSTANTS[35]*STATES[5]+ CONSTANTS[36]*STATES[5]);
resid[6] = RATES[6] -  CONSTANTS[36]*STATES[5] -  (CONSTANTS[37]+ALGEBRAIC[25])*STATES[6];
resid[7] = RATES[7] -  ALGEBRAIC[25]*STATES[6] -  CONSTANTS[38]*STATES[7];
}
void
computeVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[11] = (1.00000+pow(CONSTANTS[18]/ALGEBRAIC[0], 4.00000)+pow(ALGEBRAIC[0]/CONSTANTS[19], 3.00000))/(1.00000+1.00000/CONSTANTS[20]+pow(ALGEBRAIC[0]/CONSTANTS[19], 3.00000)+pow(CONSTANTS[18]/ALGEBRAIC[0], 4.00000));
ALGEBRAIC[12] = ( ALGEBRAIC[11]*(1.00000+pow(CONSTANTS[18]/ALGEBRAIC[0], 4.00000)+pow(ALGEBRAIC[0]/CONSTANTS[19], 3.00000)))/(1.00000+1.00000/CONSTANTS[20]+pow(ALGEBRAIC[0]/CONSTANTS[19], 3.00000)+pow(CONSTANTS[18]/ALGEBRAIC[0], 4.00000));
}
void
computeEssentialVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[2] = STATES[2]/CONSTANTS[43];
ALGEBRAIC[3] = (1.00000 - ALGEBRAIC[2]/CONSTANTS[3]) - STATES[3];
ALGEBRAIC[5] = ( CONSTANTS[7]*STATES[4])/(CONSTANTS[8]+STATES[4])+( CONSTANTS[9]*pow(STATES[4], 2.00000))/(pow(CONSTANTS[10], 2.00000)+pow(STATES[4], 2.00000));
ALGEBRAIC[6] =  ALGEBRAIC[5]*STATES[3];
ALGEBRAIC[0] = STATES[0]/CONSTANTS[43];
ALGEBRAIC[7] = ALGEBRAIC[0]/(CONSTANTS[15]+ALGEBRAIC[0]);
ALGEBRAIC[8] = ( (CONSTANTS[13]+ CONSTANTS[14]*pow(ALGEBRAIC[7], 2.00000))*ALGEBRAIC[2])/(CONSTANTS[16]+ALGEBRAIC[2]);
ALGEBRAIC[1] = STATES[1]/CONSTANTS[44];
ALGEBRAIC[9] = pow(STATES[5], 4.00000);
ALGEBRAIC[10] =  CONSTANTS[17]*CONSTANTS[0]*ALGEBRAIC[9]*(ALGEBRAIC[1] - ALGEBRAIC[0]);
ALGEBRAIC[13] =  CONSTANTS[21]*CONSTANTS[0]*(ALGEBRAIC[1] - ALGEBRAIC[0]);
ALGEBRAIC[14] =  CONSTANTS[22]*CONSTANTS[0]*(ALGEBRAIC[1] - ALGEBRAIC[0]);
ALGEBRAIC[15] = ( (( CONSTANTS[0]*CONSTANTS[23]*pow(ALGEBRAIC[0], 2.00000))/(pow(CONSTANTS[24], 2.00000)+pow(ALGEBRAIC[0], 2.00000)))*1.00000)/ALGEBRAIC[1];
ALGEBRAIC[16] =  (( CONSTANTS[26]*CONSTANTS[27])/( CONSTANTS[29]*CONSTANTS[28]))*log(CONSTANTS[30]/ALGEBRAIC[0]);
ALGEBRAIC[17] =  CONSTANTS[0]*CONSTANTS[31]*(ALGEBRAIC[16] - CONSTANTS[25]);
ALGEBRAIC[19] = ALGEBRAIC[17]/( CONSTANTS[29]*CONSTANTS[28]);
ALGEBRAIC[4] = (CONDVAR[0]>0.00000&&CONDVAR[1]<0.00000&&CONDVAR[2]>0.00000 ? CONSTANTS[6] : 0.00000);
ALGEBRAIC[18] =  CONSTANTS[0]*CONSTANTS[32]*ALGEBRAIC[4]*(ALGEBRAIC[16] - CONSTANTS[25]);
ALGEBRAIC[20] = ALGEBRAIC[18]/( CONSTANTS[29]*CONSTANTS[28]);
ALGEBRAIC[21] = ALGEBRAIC[19]+ALGEBRAIC[20];
ALGEBRAIC[22] = ( CONSTANTS[33]*CONSTANTS[0]*pow(ALGEBRAIC[0], 2.00000))/(pow(CONSTANTS[34], 2.00000)+pow(ALGEBRAIC[0], 2.00000));
ALGEBRAIC[23] = 1.00000 - (STATES[5]+STATES[6]+STATES[7]);
ALGEBRAIC[24] = ( CONSTANTS[39]*pow(ALGEBRAIC[0], 3.00000))/(pow(CONSTANTS[40], 3.00000)+pow(ALGEBRAIC[0], 3.00000));
ALGEBRAIC[25] = ( CONSTANTS[41]*ALGEBRAIC[2])/(CONSTANTS[42]+ALGEBRAIC[2]);
}
void
getStateInformation(double* SI)
{
SI[0] = 1.0;
SI[1] = 1.0;
SI[2] = 1.0;
SI[3] = 1.0;
SI[4] = 1.0;
SI[5] = 1.0;
SI[6] = 1.0;
SI[7] = 1.0;
}
void
computeRoots(double VOI, double* CONSTANTS, double* RATES, double* OLDRATES, double* STATES,
             double* OLDSTATES, double* ALGEBRAIC, double* CONDVARS)
{
CONDVAR[0] = VOI - 10.0000;
CONDVAR[1] = VOI - 25.0000;
CONDVAR[2] = CONSTANTS[5] - 0.00000;
}