/*
   There are a total of 11 entries in the algebraic variable array.
   There are a total of 4 entries in each of the rate and state variable arrays.
   There are a total of 21 entries in the constant variable array.
 */
/*
 * VOI is time in component environment (hour).
 * ALGEBRAIC[0] is Metin in component environment (flux).
 * STATES[0] is Met in component Met (micromolar).
 * ALGEBRAIC[7] is V_MS in component V_MS (flux).
 * ALGEBRAIC[10] is V_BHMT in component V_BHMT (flux).
 * ALGEBRAIC[1] is V_MATI in component V_MATI (flux).
 * ALGEBRAIC[3] is V_MATIII in component V_MATIII (flux).
 * STATES[1] is AdoMet in component AdoMet (micromolar).
 * ALGEBRAIC[6] is V_METH in component V_METH (flux).
 * ALGEBRAIC[4] is V_GNMT in component V_GNMT (flux).
 * STATES[2] is AdoHcy in component AdoHcy (micromolar).
 * ALGEBRAIC[8] is V_AH in component V_AH (flux).
 * STATES[3] is Hcy in component Hcy (micromolar).
 * ALGEBRAIC[9] is V_CBS in component V_CBS (flux).
 * CONSTANTS[0] is V_MATImax in component V_MATI (flux).
 * CONSTANTS[1] is Km_MATI in component V_MATI (micromolar).
 * CONSTANTS[2] is Ki_MATI in component V_MATI (micromolar).
 * CONSTANTS[3] is V_MATIIImax in component V_MATIII (flux).
 * ALGEBRAIC[2] is Km1_MATIII in component V_MATIII (micromolar).
 * CONSTANTS[4] is Km2_MATIII in component V_MATIII (micromolar).
 * CONSTANTS[5] is V_GNMTmax in component V_GNMT (flux).
 * CONSTANTS[6] is Km_GNMT in component V_GNMT (micromolar).
 * CONSTANTS[7] is Ki_GNMT in component V_GNMT (micromolar).
 * CONSTANTS[8] is V_METHmax in component V_METH (flux).
 * ALGEBRAIC[5] is Km1_METH in component V_METH (micromolar).
 * CONSTANTS[9] is Km2_METH_A in component V_METH (dimensionless).
 * CONSTANTS[10] is five_mTHF in component V_MS (micromolar).
 * CONSTANTS[11] is V_MSmax in component V_MS (flux).
 * CONSTANTS[12] is Kd_MS in component V_MS (micromolar).
 * CONSTANTS[13] is Km_Hcy_MS in component V_MS (micromolar).
 * CONSTANTS[14] is Km_five_mTHF_MS in component V_MS (micromolar).
 * CONSTANTS[15] is alpha1 in component V_AH (first_order_rate_constant).
 * CONSTANTS[16] is alpha2 in component V_AH (dimensionless).
 * CONSTANTS[17] is beta1 in component V_CBS (second_order_rate_constant).
 * CONSTANTS[18] is beta2 in component V_CBS (first_order_rate_constant).
 * CONSTANTS[19] is V_BHMTmax in component V_BHMT (flux).
 * CONSTANTS[20] is Km_BHMT in component V_BHMT (micromolar).
 * RATES[0] is d/dt Met in component Met (micromolar).
 * RATES[1] is d/dt AdoMet in component AdoMet (micromolar).
 * RATES[2] is d/dt AdoHcy in component AdoHcy (micromolar).
 * RATES[3] is d/dt Hcy in component Hcy (micromolar).
 */
void
initConsts(double* CONSTANTS, double* RATES, double *STATES)
{
STATES[0] = 53.5;
STATES[1] = 137.6;
STATES[2] = 13.2;
STATES[3] = 0.88;
CONSTANTS[0] = 561;
CONSTANTS[1] = 41;
CONSTANTS[2] = 50;
CONSTANTS[3] = 22870;
CONSTANTS[4] = 21.1;
CONSTANTS[5] = 10600;
CONSTANTS[6] = 4500;
CONSTANTS[7] = 20;
CONSTANTS[8] = 4521;
CONSTANTS[9] = 10;
CONSTANTS[10] = 5.2;
CONSTANTS[11] = 500;
CONSTANTS[12] = 1;
CONSTANTS[13] = 0.1;
CONSTANTS[14] = 25;
CONSTANTS[15] = 100;
CONSTANTS[16] = 10;
CONSTANTS[17] = 1.7;
CONSTANTS[18] = 30;
CONSTANTS[19] = 2500;
CONSTANTS[20] = 12;
}
void
computeRates(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[1] = CONSTANTS[0]/(1.00000+ (CONSTANTS[1]/STATES[0])*(1.00000+STATES[1]/CONSTANTS[2]));
ALGEBRAIC[2] = 20000.0/(1.00000+ 5.70000*pow(STATES[1]/(STATES[1]+600.000), 2.00000));
ALGEBRAIC[3] = CONSTANTS[3]/(1.00000+( ALGEBRAIC[2]*CONSTANTS[4])/(pow(STATES[0], 2.00000)+ STATES[0]*CONSTANTS[4]));
ALGEBRAIC[5] =  1.00000*(1.00000+STATES[2]/4.00000);
ALGEBRAIC[6] = CONSTANTS[8]/(1.00000+ALGEBRAIC[5]/STATES[1]+CONSTANTS[9]+( CONSTANTS[9]*ALGEBRAIC[5])/STATES[1]);
ALGEBRAIC[4] = ( (CONSTANTS[5]/(1.00000+pow(CONSTANTS[6]/STATES[1], 2.30000)))*1.00000)/(1.00000+STATES[2]/CONSTANTS[7]);
RATES[1] = (ALGEBRAIC[1]+ALGEBRAIC[3]) - (ALGEBRAIC[6]+ALGEBRAIC[4]);
ALGEBRAIC[8] =  CONSTANTS[15]*(STATES[2] -  CONSTANTS[16]*STATES[3]);
RATES[2] = (ALGEBRAIC[6]+ALGEBRAIC[4]) - ALGEBRAIC[8];
ALGEBRAIC[0] = (VOI<2.00000||VOI>=8.00000 ? 200.000 : VOI>=2.00000&&VOI<5.00000 ? 300.000 : VOI>=5.00000&&VOI<8.00000 ? 100.000 : 200.000);
ALGEBRAIC[7] = ( CONSTANTS[11]*CONSTANTS[10]*STATES[3])/( CONSTANTS[12]*CONSTANTS[13]+ CONSTANTS[13]*CONSTANTS[10]+ CONSTANTS[14]*STATES[3]+ CONSTANTS[10]*STATES[3]);
ALGEBRAIC[10] = ( (0.700000 -  0.0250000*((STATES[1]+STATES[2]) - 150.000))*CONSTANTS[19]*STATES[3])/(CONSTANTS[20]+STATES[3]);
RATES[0] = (ALGEBRAIC[7]+ALGEBRAIC[10]+ALGEBRAIC[0]) - (ALGEBRAIC[1]+ALGEBRAIC[3]);
ALGEBRAIC[9] =  ( CONSTANTS[17]*(STATES[1]+STATES[2]) - CONSTANTS[18])*STATES[3];
RATES[3] = ALGEBRAIC[8] - (ALGEBRAIC[9]+ALGEBRAIC[7]+ALGEBRAIC[10]);
}
void
computeVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[1] = CONSTANTS[0]/(1.00000+ (CONSTANTS[1]/STATES[0])*(1.00000+STATES[1]/CONSTANTS[2]));
ALGEBRAIC[2] = 20000.0/(1.00000+ 5.70000*pow(STATES[1]/(STATES[1]+600.000), 2.00000));
ALGEBRAIC[3] = CONSTANTS[3]/(1.00000+( ALGEBRAIC[2]*CONSTANTS[4])/(pow(STATES[0], 2.00000)+ STATES[0]*CONSTANTS[4]));
ALGEBRAIC[5] =  1.00000*(1.00000+STATES[2]/4.00000);
ALGEBRAIC[6] = CONSTANTS[8]/(1.00000+ALGEBRAIC[5]/STATES[1]+CONSTANTS[9]+( CONSTANTS[9]*ALGEBRAIC[5])/STATES[1]);
ALGEBRAIC[4] = ( (CONSTANTS[5]/(1.00000+pow(CONSTANTS[6]/STATES[1], 2.30000)))*1.00000)/(1.00000+STATES[2]/CONSTANTS[7]);
ALGEBRAIC[8] =  CONSTANTS[15]*(STATES[2] -  CONSTANTS[16]*STATES[3]);
ALGEBRAIC[0] = (VOI<2.00000||VOI>=8.00000 ? 200.000 : VOI>=2.00000&&VOI<5.00000 ? 300.000 : VOI>=5.00000&&VOI<8.00000 ? 100.000 : 200.000);
ALGEBRAIC[7] = ( CONSTANTS[11]*CONSTANTS[10]*STATES[3])/( CONSTANTS[12]*CONSTANTS[13]+ CONSTANTS[13]*CONSTANTS[10]+ CONSTANTS[14]*STATES[3]+ CONSTANTS[10]*STATES[3]);
ALGEBRAIC[10] = ( (0.700000 -  0.0250000*((STATES[1]+STATES[2]) - 150.000))*CONSTANTS[19]*STATES[3])/(CONSTANTS[20]+STATES[3]);
ALGEBRAIC[9] =  ( CONSTANTS[17]*(STATES[1]+STATES[2]) - CONSTANTS[18])*STATES[3];
}