/*
   There are a total of 3 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 13 entries in the constant variable array.
 */
/*
 * VOI is time in component enviroment (second).
 * STATES[0] is P in component P (micromolar).
 * CONSTANTS[0] is V_p in component P (per_second).
 * CONSTANTS[1] is k_p in component P (micromolar).
 * CONSTANTS[2] is IPR_3_flux in component P (flux).
 * STATES[1] is c in component c (micromolar).
 * ALGEBRAIC[2] is J_flux in component J_flux (flux).
 * ALGEBRAIC[0] is J_pump in component J_pump (flux).
 * CONSTANTS[12] is J_leak in component J_leak (flux).
 * CONSTANTS[3] is k_flux in component J_flux (micromolar_per_second).
 * ALGEBRAIC[1] is mu in component mu (dimensionless).
 * STATES[2] is h in component h (dimensionless).
 * CONSTANTS[4] is b in component J_flux (dimensionless).
 * CONSTANTS[5] is k_1 in component J_flux (micromolar).
 * CONSTANTS[6] is gamma in component J_pump (micromolar_per_second).
 * CONSTANTS[7] is k_gamma in component J_pump (micromolar).
 * CONSTANTS[8] is beta in component J_leak (flux).
 * CONSTANTS[9] is k_mu in component mu (micromolar).
 * CONSTANTS[10] is k_2 in component h (micromolar).
 * CONSTANTS[11] is tau_h in component h (second).
 * RATES[0] is d/dt P in component P (micromolar).
 * RATES[1] is d/dt c in component c (micromolar).
 * RATES[2] is d/dt h in component h (dimensionless).
 */
void
initConsts(double* CONSTANTS, double* RATES, double *STATES)
{
STATES[0] = 0;
CONSTANTS[0] = 0.08;
CONSTANTS[1] = 1;
CONSTANTS[2] = 0.72;
STATES[1] = 0.3;
CONSTANTS[3] = 3;
STATES[2] = 1;
CONSTANTS[4] = 0.11;
CONSTANTS[5] = 0.7;
CONSTANTS[6] = 1;
CONSTANTS[7] = 0.27;
CONSTANTS[8] = 0.15;
CONSTANTS[9] = 0.01;
CONSTANTS[10] = 0.7;
CONSTANTS[11] = 0.2;
CONSTANTS[12] = CONSTANTS[8];
}
void
computeRates(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
RATES[0] = (VOI<=15.0000 ? CONSTANTS[2] - ( CONSTANTS[0]*STATES[0]*CONSTANTS[1])/(CONSTANTS[1]+STATES[0]) : ( - CONSTANTS[0]*STATES[0]*CONSTANTS[1])/(CONSTANTS[1]+STATES[0]));
RATES[2] = (pow(CONSTANTS[10], 2.00000)/(pow(CONSTANTS[10], 2.00000)+pow(STATES[1], 2.00000)) - STATES[2])/CONSTANTS[11];
ALGEBRAIC[1] = pow(STATES[0], 3.00000)/(pow(CONSTANTS[9], 3.00000)+pow(STATES[0], 3.00000));
ALGEBRAIC[2] =  CONSTANTS[3]*ALGEBRAIC[1]*STATES[2]*(CONSTANTS[4]+( (1.00000 - CONSTANTS[4])*STATES[1])/(CONSTANTS[5]+STATES[1]));
ALGEBRAIC[0] = ( CONSTANTS[6]*pow(STATES[1], 2.00000))/(pow(CONSTANTS[7], 2.00000)+pow(STATES[1], 2.00000));
RATES[1] = (ALGEBRAIC[2] - ALGEBRAIC[0])+CONSTANTS[12];
}
void
computeVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[1] = pow(STATES[0], 3.00000)/(pow(CONSTANTS[9], 3.00000)+pow(STATES[0], 3.00000));
ALGEBRAIC[2] =  CONSTANTS[3]*ALGEBRAIC[1]*STATES[2]*(CONSTANTS[4]+( (1.00000 - CONSTANTS[4])*STATES[1])/(CONSTANTS[5]+STATES[1]));
ALGEBRAIC[0] = ( CONSTANTS[6]*pow(STATES[1], 2.00000))/(pow(CONSTANTS[7], 2.00000)+pow(STATES[1], 2.00000));
}