/*
   There are a total of 6 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 20 entries in the constant variable array.
 */
/*
 * VOI is time in component environment (ms).
 * ALGEBRAIC[1] is Ca_i in component intracellular_ion_concentrations (uM).
 * ALGEBRAIC[0] is mtime in component intracellular_ion_concentrations (dimensionless).
 * STATES[0] is xb in component crossbridges (dimensionless).
 * STATES[1] is TRPN in component troponin (dimensionless).
 * CONSTANTS[0] is k_xb in component crossbridges (per_ms).
 * CONSTANTS[1] is nperm in component crossbridges (dimensionless).
 * CONSTANTS[2] is perm50 in component crossbridges (dimensionless).
 * ALGEBRAIC[2] is permtot in component crossbridges (dimensionless).
 * CONSTANTS[3] is Ca_50ref in component troponin (uM).
 * CONSTANTS[4] is beta_1 in component troponin (dimensionless).
 * CONSTANTS[5] is k_off in component troponin (per_ms).
 * CONSTANTS[6] is n_TRPN in component troponin (dimensionless).
 * CONSTANTS[16] is lambda_m in component filament_overlap (dimensionless).
 * CONSTANTS[17] is Ca_50 in component troponin (uM).
 * CONSTANTS[14] is lambda in component Myofilaments (dimensionless).
 * CONSTANTS[15] is dlambdadt in component Myofilaments (per_ms).
 * CONSTANTS[19] is overlap in component filament_overlap (dimensionless).
 * CONSTANTS[7] is beta_0 in component filament_overlap (dimensionless).
 * CONSTANTS[18] is lambda_s in component filament_overlap (dimensionless).
 * CONSTANTS[8] is T_ref in component isometric_tension (kPa).
 * ALGEBRAIC[3] is T_0 in component isometric_tension (kPa).
 * ALGEBRAIC[4] is Q in component dynamic_stiffness (dimensionless).
 * CONSTANTS[9] is a in component dynamic_stiffness (dimensionless).
 * STATES[2] is Q_1 in component dynamic_stiffness (dimensionless).
 * STATES[3] is Q_2 in component dynamic_stiffness (dimensionless).
 * CONSTANTS[10] is A_1 in component dynamic_stiffness (dimensionless).
 * CONSTANTS[11] is A_2 in component dynamic_stiffness (dimensionless).
 * CONSTANTS[12] is alpha_1 in component dynamic_stiffness (per_ms).
 * CONSTANTS[13] is alpha_2 in component dynamic_stiffness (per_ms).
 * ALGEBRAIC[5] is Tension in component dynamic_stiffness (kPa).
 * RATES[0] is d/dt xb in component crossbridges (dimensionless).
 * RATES[1] is d/dt TRPN in component troponin (dimensionless).
 * RATES[2] is d/dt Q_1 in component dynamic_stiffness (dimensionless).
 * RATES[3] is d/dt Q_2 in component dynamic_stiffness (dimensionless).
 * There are a total of 4 condition variables.
 */
void
initConsts(double* CONSTANTS, double* RATES, double *STATES)
{
STATES[0] = 0.00046;
STATES[1] = 0.0752;
CONSTANTS[0] = 0.1;
CONSTANTS[1] = 5;
CONSTANTS[2] = 0.35;
CONSTANTS[3] = 0.8;
CONSTANTS[4] = -1.5;
CONSTANTS[5] = 0.1;
CONSTANTS[6] = 2;
CONSTANTS[7] = 1.65;
CONSTANTS[8] = 120;
CONSTANTS[9] = 0.35;
STATES[2] = 0;
STATES[3] = 0;
CONSTANTS[10] = -29;
CONSTANTS[11] = 116;
CONSTANTS[12] = 0.1;
CONSTANTS[13] = 0.5;
CONSTANTS[14] = 1.00000;
CONSTANTS[15] = 0.00000;
CONSTANTS[16] = (CONSTANTS[14]>1.20000 ? 1.20000 : CONSTANTS[14]);
CONSTANTS[17] =  CONSTANTS[3]*(1.00000+ CONSTANTS[4]*(CONSTANTS[16] - 1.00000));
CONSTANTS[18] = (CONSTANTS[16]>=0.870000 ? 0.870000 : CONSTANTS[16]);
CONSTANTS[19] = 1.00000+ CONSTANTS[7]*((CONSTANTS[16]+CONSTANTS[18]) - 1.87000);
RATES[0] = 0.1001;
RATES[1] = 0.1001;
RATES[2] = 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] -  CONSTANTS[0]*( ALGEBRAIC[2]*(1.00000 - STATES[0]) -  (1.00000/ALGEBRAIC[2])*STATES[0]);
resid[1] = RATES[1] -  CONSTANTS[5]*( pow(ALGEBRAIC[1]/CONSTANTS[17], CONSTANTS[6])*(1.00000 - STATES[1]) - STATES[1]);
resid[2] = RATES[2] -  CONSTANTS[10]*CONSTANTS[15] -  CONSTANTS[12]*STATES[2];
resid[3] = RATES[3] -  CONSTANTS[11]*CONSTANTS[15] -  CONSTANTS[13]*STATES[3];
}
void
computeVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[3] =  CONSTANTS[8]*STATES[0]*CONSTANTS[19];
ALGEBRAIC[4] = STATES[2]+STATES[3];
ALGEBRAIC[5] = (CONDVAR[3]<0.00000 ? ( ALGEBRAIC[3]*( CONSTANTS[9]*ALGEBRAIC[4]+1.00000))/(1.00000 - ALGEBRAIC[4]) : ( ALGEBRAIC[3]*(1.00000+ (CONSTANTS[9]+2.00000)*ALGEBRAIC[4]))/(1.00000+ALGEBRAIC[4]));
}
void
computeEssentialVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[0] = (VOI -  167.000*floor(VOI/167.000))/1.00000;
ALGEBRAIC[1] = (CONDVAR[0]>=0.00000&&CONDVAR[1]<0.00000 ?  1.00000*1.85358e-05*pow(ALGEBRAIC[0], 3.00000)+ - 0.00159034*pow(ALGEBRAIC[0], 2.00000)+ 0.0436459*pow(ALGEBRAIC[0], 1.00000)+0.167079 : CONDVAR[2]>=0.00000 ? (( 1.00000*- 5.74585e-08*pow(ALGEBRAIC[0], 3.00000)+ 3.11222e-05*pow(ALGEBRAIC[0], 2.00000)) -  0.00661849*pow(ALGEBRAIC[0], 1.00000))+0.720442 : 0.216000);
ALGEBRAIC[2] =  pow(pow(STATES[1]/CONSTANTS[2], CONSTANTS[1]), 1.0 / 2);
}
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] = ALGEBRAIC[0] - 1.17000;
CONDVAR[1] = ALGEBRAIC[0] - 30.8400;
CONDVAR[2] = ALGEBRAIC[0] - 30.8400;
CONDVAR[3] = ALGEBRAIC[4] - 0.00000;
}