/*
   There are a total of 13 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 14 entries in the constant variable array.
 */
/*
 * CONSTANTS[0] is stimPeriod in component stimulus_protocol_params (ms).
 * CONSTANTS[1] is stimDuration in component stimulus_protocol_params (ms).
 * CONSTANTS[2] is stimCurrent in component stimulus_protocol_params (uA_per_cmsq).
 * VOI is time in component time (ms).
 * ALGEBRAIC[11] is INa in component INa (uA_per_cmsq).
 * ALGEBRAIC[12] is IK in component IK (uA_per_cmsq).
 * ALGEBRAIC[4] is Ileak in component Ileak (uA_per_cmsq).
 * STATES[0] is V in component action_potential (mV).
 * ALGEBRAIC[0] is minus_V in component action_potential (mV).
 * STATES[1] is X in component hh_gating_variable (dimensionless).
 * STATES[2] is X in component hh_gating_variable (dimensionless).
 * STATES[3] is X in component hh_gating_variable (dimensionless).
 * CONSTANTS[3] is V_initial in component initial_conditions (mV).
 * CONSTANTS[4] is n_initial in component initial_conditions (dimensionless).
 * CONSTANTS[5] is m_initial in component initial_conditions (dimensionless).
 * CONSTANTS[6] is h_initial in component initial_conditions (dimensionless).
 * ALGEBRAIC[1] is IStim in component stimulus_protocol (uA_per_cmsq).
 * CONSTANTS[7] is VNa in component parameters (mV).
 * CONSTANTS[8] is VK in component parameters (mV).
 * CONSTANTS[9] is Vleak in component parameters (mV).
 * CONSTANTS[10] is gNa_max in component parameters (mS_per_cmsq).
 * CONSTANTS[11] is gK_max in component parameters (mS_per_cmsq).
 * CONSTANTS[12] is gleak_max in component parameters (mS_per_cmsq).
 * CONSTANTS[13] is Cm in component parameters (uF_per_cmsq).
 * ALGEBRAIC[2] is gNa in component INa (mS_per_cmsq).
 * ALGEBRAIC[3] is gK in component IK (mS_per_cmsq).
 * ALGEBRAIC[5] is alpha in component alpha_m (per_ms).
 * ALGEBRAIC[6] is beta in component beta_m (per_ms).
 * ALGEBRAIC[7] is alpha in component alpha_h (per_ms).
 * ALGEBRAIC[8] is beta in component beta_h (per_ms).
 * ALGEBRAIC[9] is alpha in component alpha_n (per_ms).
 * ALGEBRAIC[10] is beta in component beta_n (per_ms).
 * RATES[0] is d/dt V in component action_potential (mV).
 * RATES[2] is d/dt X in component hh_gating_variable (dimensionless).
 * RATES[1] is d/dt X in component hh_gating_variable (dimensionless).
 * RATES[3] is d/dt X in component hh_gating_variable (dimensionless).
 * There are a total of 1 condition variables.
 */
void
initConsts(double* CONSTANTS, double* RATES, double *STATES)
{
CONSTANTS[0] = 50;
CONSTANTS[1] = 0.5;
CONSTANTS[2] = -10.0;
CONSTANTS[3] = 0.0;
CONSTANTS[4] = 0.315;
CONSTANTS[5] = 0.042;
CONSTANTS[6] = 0.608;
CONSTANTS[7] = -115.0;
CONSTANTS[8] = 12.0;
CONSTANTS[9] = -10.613;
CONSTANTS[10] = 120.0;
CONSTANTS[11] = 36.0;
CONSTANTS[12] = 0.3;
CONSTANTS[13] = 1.0;
STATES[0] = CONSTANTS[3];
STATES[1] = CONSTANTS[6];
STATES[2] = CONSTANTS[5];
STATES[3] = CONSTANTS[4];
RATES[0] = 0.1001;
RATES[2] = 0.1001;
RATES[1] = 0.1001;
RATES[3] = 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[1] - (ALGEBRAIC[11]+ALGEBRAIC[12]+ALGEBRAIC[4]))/CONSTANTS[13];
resid[1] = RATES[2] -  ALGEBRAIC[5]*(1.00000 - STATES[2]) -  ALGEBRAIC[6]*STATES[2];
resid[2] = RATES[1] -  ALGEBRAIC[7]*(1.00000 - STATES[1]) -  ALGEBRAIC[8]*STATES[1];
resid[3] = RATES[3] -  ALGEBRAIC[9]*(1.00000 - STATES[3]) -  ALGEBRAIC[10]*STATES[3];
}
void
computeVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[0] = - STATES[0];
}
void
computeEssentialVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[1] = (CONDVAR[0]<0.00000 ? CONSTANTS[2] : 0.00000);
ALGEBRAIC[4] =  CONSTANTS[12]*(STATES[0] - CONSTANTS[9]);
ALGEBRAIC[5] = ( 0.100000*(STATES[0]+25.0000))/(exp((STATES[0]+25.0000)/10.0000) - 1.00000);
ALGEBRAIC[6] =  4.00000*exp(STATES[0]/18.0000);
ALGEBRAIC[7] =  0.0700000*exp(STATES[0]/20.0000);
ALGEBRAIC[8] = 1.00000/(exp((STATES[0]+30.0000)/10.0000)+1.00000);
ALGEBRAIC[9] = ( 0.0100000*(STATES[0]+10.0000))/(exp((STATES[0]+10.0000)/10.0000) - 1.00000);
ALGEBRAIC[10] =  0.125000*exp(STATES[0]/80.0000);
ALGEBRAIC[2] =  CONSTANTS[10]*pow(STATES[2], 3.00000)*STATES[1];
ALGEBRAIC[11] =  ALGEBRAIC[2]*(STATES[0] - CONSTANTS[7]);
ALGEBRAIC[3] =  CONSTANTS[11]*pow(STATES[3], 4.00000);
ALGEBRAIC[12] =  ALGEBRAIC[3]*(STATES[0] - CONSTANTS[8]);
}
void
getStateInformation(double* SI)
{
SI[0] = 1.0;
SI[1] = 1.0;
SI[2] = 1.0;
SI[3] = 1.0;
}
void
computeRoots(double VOI, double* CONSTANTS, double* RATES, double* OLDRATES, double* STATES,
             double* OLDSTATES, double* ALGEBRAIC, double* CONDVARS)
{
CONDVAR[0] =  (int)(VOI) % (int)(CONSTANTS[0]) - CONSTANTS[1];
}