/*
   There are a total of 9 entries in the algebraic variable array.
   There are a total of 5 entries in each of the rate and state variable arrays.
   There are a total of 20 entries in the constant variable array.
 */
/*
 * VOI is time in component environment (millisecond).
 * CONSTANTS[0] is tau_c in component nucleotides (second).
 * CONSTANTS[1] is eta in component nucleotides (dimensionless).
 * CONSTANTS[2] is v in component nucleotides (dimensionless).
 * CONSTANTS[3] is k in component nucleotides (dimensionless).
 * ALGEBRAIC[0] is phi in component nucleotides (dimensionless).
 * STATES[0] is ADP in component nucleotides (dimensionless).
 * STATES[1] is ATP in component nucleotides (dimensionless).
 * CONSTANTS[4] is C_m in component membrane (femtofarad).
 * ALGEBRAIC[3] is I_Ca in component Ca_current (femtoampere).
 * ALGEBRAIC[4] is I_K in component K_current (femtoampere).
 * ALGEBRAIC[7] is I_KCa in component Ca_activated_K_current (femtoampere).
 * ALGEBRAIC[8] is I_KATP in component ATP_sensitive_K_current (femtoampere).
 * STATES[2] is V in component membrane (millivolt).
 * CONSTANTS[5] is g_Ca_ in component Ca_current (picosiemens).
 * CONSTANTS[6] is V_Ca in component Ca_current (millivolt).
 * CONSTANTS[7] is v_m in component Ca_current (millivolt).
 * CONSTANTS[8] is s_m in component Ca_current (millivolt).
 * ALGEBRAIC[1] is m_infinity in component Ca_current (dimensionless).
 * CONSTANTS[9] is g_K_ in component K_current (picosiemens).
 * CONSTANTS[10] is V_K in component K_current (millivolt).
 * STATES[3] is n in component K_channel_activation (dimensionless).
 * CONSTANTS[11] is g_KCa_ in component Ca_activated_K_current (picosiemens).
 * CONSTANTS[12] is k_D in component Ca_activated_K_current (micromolar).
 * STATES[4] is c in component cytosolic_Ca (micromolar).
 * ALGEBRAIC[6] is omega in component Ca_activated_K_current (dimensionless).
 * CONSTANTS[13] is g_KATP_ in component ATP_sensitive_K_current (picosiemens).
 * CONSTANTS[14] is tau_n in component K_channel_activation (millisecond).
 * CONSTANTS[15] is v_n in component K_channel_activation (millivolt).
 * CONSTANTS[16] is s_n in component K_channel_activation (millivolt).
 * ALGEBRAIC[2] is n_infinity in component K_channel_activation (dimensionless).
 * ALGEBRAIC[5] is J_mem in component Ca_influx (micromolar_per_ms).
 * CONSTANTS[17] is f in component Ca_influx (dimensionless).
 * CONSTANTS[18] is alpha in component Ca_influx (micromolar_per_fA_ms).
 * CONSTANTS[19] is k_c in component Ca_influx (per_millisecond).
 * RATES[1] is d/dt ATP in component nucleotides (dimensionless).
 * RATES[0] is d/dt ADP in component nucleotides (dimensionless).
 * RATES[2] is d/dt V in component membrane (millivolt).
 * RATES[3] is d/dt n in component K_channel_activation (dimensionless).
 * RATES[4] is d/dt c in component cytosolic_Ca (micromolar).
 */
void
initConsts(double* CONSTANTS, double* RATES, double *STATES)
{
CONSTANTS[0] = 1200;
CONSTANTS[1] = 185;
CONSTANTS[2] = 10;
CONSTANTS[3] = 20;
STATES[0] = 0.085817;
STATES[1] = 2.1047;
CONSTANTS[4] = 5300;
STATES[2] = -67.018;
CONSTANTS[5] = 1200;
CONSTANTS[6] = 25;
CONSTANTS[7] = -20;
CONSTANTS[8] = 12;
CONSTANTS[9] = 3000;
CONSTANTS[10] = -75;
STATES[3] = 0.00011;
CONSTANTS[11] = 300;
CONSTANTS[12] = 0.3;
STATES[4] = 0.15666;
CONSTANTS[13] = 350;
CONSTANTS[14] = 16;
CONSTANTS[15] = -16;
CONSTANTS[16] = 5.6;
CONSTANTS[17] = 0.001;
CONSTANTS[18] = 0.00000225;
CONSTANTS[19] = 0.1;
}
void
computeRates(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[0] =  STATES[1]*pow(1.00000+ CONSTANTS[3]*STATES[0], 2.00000);
RATES[1] = (CONSTANTS[2] - ALGEBRAIC[0])/( 1000.00*CONSTANTS[0]);
RATES[0] = (ALGEBRAIC[0] -  CONSTANTS[1]*STATES[0])/( 1000.00*CONSTANTS[0]);
ALGEBRAIC[2] = 1.00000/(1.00000+exp((CONSTANTS[15] - STATES[2])/CONSTANTS[16]));
RATES[3] = (ALGEBRAIC[2] - STATES[3])/CONSTANTS[14];
ALGEBRAIC[1] = 1.00000/(1.00000+exp((CONSTANTS[7] - STATES[2])/CONSTANTS[8]));
ALGEBRAIC[3] =  CONSTANTS[5]*ALGEBRAIC[1]*(STATES[2] - CONSTANTS[6]);
ALGEBRAIC[5] =  - CONSTANTS[17]*( CONSTANTS[18]*ALGEBRAIC[3]+ CONSTANTS[19]*STATES[4]);
RATES[4] = ALGEBRAIC[5];
ALGEBRAIC[4] =  CONSTANTS[9]*STATES[3]*(STATES[2] - CONSTANTS[10]);
ALGEBRAIC[6] = 1.00000/(1.00000+CONSTANTS[12]/STATES[4]);
ALGEBRAIC[7] =  CONSTANTS[11]*ALGEBRAIC[6]*(STATES[2] - CONSTANTS[10]);
ALGEBRAIC[8] = ( (STATES[2] - CONSTANTS[10])*CONSTANTS[13])/STATES[1];
RATES[2] = - (ALGEBRAIC[3]+ALGEBRAIC[4]+ALGEBRAIC[7]+ALGEBRAIC[8])/CONSTANTS[4];
}
void
computeVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[0] =  STATES[1]*pow(1.00000+ CONSTANTS[3]*STATES[0], 2.00000);
ALGEBRAIC[2] = 1.00000/(1.00000+exp((CONSTANTS[15] - STATES[2])/CONSTANTS[16]));
ALGEBRAIC[1] = 1.00000/(1.00000+exp((CONSTANTS[7] - STATES[2])/CONSTANTS[8]));
ALGEBRAIC[3] =  CONSTANTS[5]*ALGEBRAIC[1]*(STATES[2] - CONSTANTS[6]);
ALGEBRAIC[5] =  - CONSTANTS[17]*( CONSTANTS[18]*ALGEBRAIC[3]+ CONSTANTS[19]*STATES[4]);
ALGEBRAIC[4] =  CONSTANTS[9]*STATES[3]*(STATES[2] - CONSTANTS[10]);
ALGEBRAIC[6] = 1.00000/(1.00000+CONSTANTS[12]/STATES[4]);
ALGEBRAIC[7] =  CONSTANTS[11]*ALGEBRAIC[6]*(STATES[2] - CONSTANTS[10]);
ALGEBRAIC[8] = ( (STATES[2] - CONSTANTS[10])*CONSTANTS[13])/STATES[1];
}