/*
   There are a total of 9 entries in the algebraic variable array.
   There are a total of 16 entries in each of the rate and state variable arrays.
   There are a total of 14 entries in the constant variable array.
 */
/*
 * VOI is time in component environment (day).
 * STATES[0] is B0 in component B0 (cells_per_GC).
 * CONSTANTS[0] is pr in component kinetic_parameters (dimensionless).
 * CONSTANTS[1] is mu in component kinetic_parameters (first_order_rate_constant).
 * CONSTANTS[2] is rho in component kinetic_parameters (first_order_rate_constant).
 * CONSTANTS[3] is delta_B in component kinetic_parameters (first_order_rate_constant).
 * ALGEBRAIC[4] is CT_star in component CT_star (cells_per_GC).
 * STATES[1] is B1 in component B1 (cells_per_GC).
 * ALGEBRAIC[5] is alpha_B in component alpha_B (dimensionless).
 * STATES[2] is B2 in component B2 (cells_per_GC).
 * STATES[3] is B3 in component B3 (cells_per_GC).
 * STATES[4] is B4 in component B4 (cells_per_GC).
 * STATES[5] is B5 in component B5 (cells_per_GC).
 * STATES[6] is B6 in component B6 (cells_per_GC).
 * STATES[7] is B7 in component B7 (cells_per_GC).
 * STATES[8] is B8 in component B8 (cells_per_GC).
 * STATES[9] is B9 in component B9 (cells_per_GC).
 * STATES[10] is B10 in component B10 (cells_per_GC).
 * ALGEBRAIC[0] is B_sum in component centroblasts_sum (cells_per_GC).
 * STATES[11] is C in component C (cells_per_GC).
 * CONSTANTS[4] is d in component C (first_order_rate_constant).
 * STATES[12] is C_star in component C_star (cells_per_GC).
 * ALGEBRAIC[3] is CA in component CA (cells_per_GC).
 * ALGEBRAIC[1] is C_starsum in component centrocytes_sum (cells_per_GC).
 * STATES[13] is M in component M (cells_per_GC).
 * STATES[14] is A in component A (cells_per_GC).
 * CONSTANTS[5] is z in component A (first_order_rate_constant).
 * CONSTANTS[6] is u in component A (dimensionless).
 * ALGEBRAIC[2] is log_A in component A (dimensionless).
 * STATES[15] is T in component T (cells_per_GC).
 * CONSTANTS[7] is p in component T (first_order_rate_constant).
 * CONSTANTS[8] is sigma in component T (first_order_rate_constant).
 * CONSTANTS[9] is delta_T in component T (first_order_rate_constant).
 * ALGEBRAIC[6] is alpha_T in component alpha_T (dimensionless).
 * CONSTANTS[10] is SA in component CA (dimensionless).
 * CONSTANTS[11] is ST in component CT_star (dimensionless).
 * CONSTANTS[12] is KB in component alpha_B (dimensionless).
 * CONSTANTS[13] is KT in component alpha_T (dimensionless).
 * ALGEBRAIC[7] is total in component total (cells_per_GC).
 * ALGEBRAIC[8] is log_total in component total (dimensionless).
 * RATES[0] is d/dt B0 in component B0 (cells_per_GC).
 * RATES[1] is d/dt B1 in component B1 (cells_per_GC).
 * RATES[2] is d/dt B2 in component B2 (cells_per_GC).
 * RATES[3] is d/dt B3 in component B3 (cells_per_GC).
 * RATES[4] is d/dt B4 in component B4 (cells_per_GC).
 * RATES[5] is d/dt B5 in component B5 (cells_per_GC).
 * RATES[6] is d/dt B6 in component B6 (cells_per_GC).
 * RATES[7] is d/dt B7 in component B7 (cells_per_GC).
 * RATES[8] is d/dt B8 in component B8 (cells_per_GC).
 * RATES[9] is d/dt B9 in component B9 (cells_per_GC).
 * RATES[10] is d/dt B10 in component B10 (cells_per_GC).
 * RATES[11] is d/dt C in component C (cells_per_GC).
 * RATES[12] is d/dt C_star in component C_star (cells_per_GC).
 * RATES[13] is d/dt M in component M (cells_per_GC).
 * RATES[14] is d/dt A in component A (cells_per_GC).
 * RATES[15] is d/dt T in component T (cells_per_GC).
 * There are a total of 0 condition variables.
 */
void
initConsts(double* CONSTANTS, double* RATES, double *STATES)
{
STATES[0] = 3;
CONSTANTS[0] = 0.15;
CONSTANTS[1] = 3;
CONSTANTS[2] = 4;
CONSTANTS[3] = 0.8;
STATES[1] = 0;
STATES[2] = 0;
STATES[3] = 0;
STATES[4] = 0;
STATES[5] = 0;
STATES[6] = 0;
STATES[7] = 0;
STATES[8] = 0;
STATES[9] = 0;
STATES[10] = 0;
STATES[11] = 0;
CONSTANTS[4] = 2;
STATES[12] = 0;
STATES[13] = 0;
STATES[14] = 500;
CONSTANTS[5] = 0.02;
CONSTANTS[6] = 0.15;
STATES[15] = 0;
CONSTANTS[7] = 2;
CONSTANTS[8] = 5;
CONSTANTS[9] = 0.8;
CONSTANTS[10] = 500;
CONSTANTS[11] = 50;
CONSTANTS[12] = 1e4;
CONSTANTS[13] = 100;
RATES[0] = 0.1001;
RATES[1] = 0.1001;
RATES[2] = 0.1001;
RATES[3] = 0.1001;
RATES[4] = 0.1001;
RATES[5] = 0.1001;
RATES[6] = 0.1001;
RATES[7] = 0.1001;
RATES[8] = 0.1001;
RATES[9] = 0.1001;
RATES[10] = 0.1001;
RATES[11] = 0.1001;
RATES[12] = 0.1001;
RATES[13] = 0.1001;
RATES[14] = 0.1001;
RATES[15] = 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]*CONSTANTS[1]*ALGEBRAIC[4] - ( CONSTANTS[2]*STATES[0]+ CONSTANTS[3]*STATES[0]);
resid[1] = RATES[1] -  CONSTANTS[2]*(1.00000+ALGEBRAIC[5])*STATES[0] - ( CONSTANTS[2]*STATES[1]+ CONSTANTS[3]*STATES[1]);
resid[2] = RATES[2] -  CONSTANTS[2]*(1.00000+ALGEBRAIC[5])*STATES[1] - ( CONSTANTS[2]*STATES[2]+ CONSTANTS[3]*STATES[2]);
resid[3] = RATES[3] -  CONSTANTS[2]*(1.00000+ALGEBRAIC[5])*STATES[2] - ( CONSTANTS[2]*STATES[3]+ CONSTANTS[3]*STATES[3]);
resid[4] = RATES[4] -  CONSTANTS[2]*(1.00000+ALGEBRAIC[5])*STATES[3] - ( CONSTANTS[2]*STATES[4]+ CONSTANTS[3]*STATES[4]);
resid[5] = RATES[5] -  CONSTANTS[2]*(1.00000+ALGEBRAIC[5])*STATES[4] - ( CONSTANTS[2]*STATES[5]+ CONSTANTS[3]*STATES[5]);
resid[6] = RATES[6] -  CONSTANTS[2]*(1.00000+ALGEBRAIC[5])*STATES[5] - ( CONSTANTS[2]*STATES[6]+ CONSTANTS[3]*STATES[6]);
resid[7] = RATES[7] -  CONSTANTS[2]*(1.00000+ALGEBRAIC[5])*STATES[6] - ( CONSTANTS[2]*STATES[7]+ CONSTANTS[3]*STATES[7]);
resid[8] = RATES[8] -  CONSTANTS[2]*(1.00000+ALGEBRAIC[5])*STATES[7] - ( CONSTANTS[2]*STATES[8]+ CONSTANTS[3]*STATES[8]);
resid[9] = RATES[9] -  CONSTANTS[2]*(1.00000+ALGEBRAIC[5])*STATES[8] - ( CONSTANTS[2]*STATES[9]+ CONSTANTS[3]*STATES[9]);
resid[10] = RATES[10] -  CONSTANTS[2]*(1.00000+ALGEBRAIC[5])*STATES[9] - ( CONSTANTS[2]*STATES[10]+ CONSTANTS[3]*STATES[10]);
resid[11] = RATES[11] -  CONSTANTS[4]*CONSTANTS[2]*STATES[10]*1.00000 -  CONSTANTS[1]*STATES[11];
resid[12] = RATES[12] -  CONSTANTS[1]*ALGEBRAIC[3] -  CONSTANTS[1]*STATES[12];
resid[13] = RATES[13] -  (1.00000 - CONSTANTS[0])*CONSTANTS[1]*ALGEBRAIC[4];
resid[14] = RATES[14] -  - CONSTANTS[5]*STATES[14] -  CONSTANTS[6]*ALGEBRAIC[3]*1.00000;
resid[15] = RATES[15] - ( CONSTANTS[8]*1.00000+ CONSTANTS[7]*ALGEBRAIC[6]*ALGEBRAIC[4]) -  CONSTANTS[9]*STATES[15];
}
void
computeVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[1] = STATES[11]+STATES[12];
ALGEBRAIC[2] = arbitrary_log(STATES[14]/1.00000, 10);
ALGEBRAIC[7] = ALGEBRAIC[0]+ALGEBRAIC[1];
ALGEBRAIC[8] = arbitrary_log(ALGEBRAIC[7]/1.00000+1.00000e-12, 10);
}
void
computeEssentialVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[3] = ( STATES[11]*STATES[14])/( CONSTANTS[10]*1.00000+STATES[14]);
ALGEBRAIC[4] = ( STATES[12]*STATES[15])/( CONSTANTS[11]*1.00000+STATES[12]);
ALGEBRAIC[0] = STATES[1]+STATES[1]+STATES[3]+STATES[4]+STATES[5]+STATES[6]+STATES[6]+STATES[7]+STATES[8]+STATES[9]+STATES[10];
ALGEBRAIC[5] = CONSTANTS[12]/(CONSTANTS[12]+ALGEBRAIC[0]/1.00000);
ALGEBRAIC[6] = CONSTANTS[13]/(CONSTANTS[13]+STATES[15]/1.00000);
}
void
getStateInformation(double* SI)
{
SI[0] = 1.0;
SI[1] = 1.0;
SI[2] = 1.0;
SI[3] = 1.0;
SI[4] = 1.0;
SI[5] = 1.0;
SI[6] = 1.0;
SI[7] = 1.0;
SI[8] = 1.0;
SI[9] = 1.0;
SI[10] = 1.0;
SI[11] = 1.0;
SI[12] = 1.0;
SI[13] = 1.0;
SI[14] = 1.0;
SI[15] = 1.0;
}
void
computeRoots(double VOI, double* CONSTANTS, double* RATES, double* OLDRATES, double* STATES,
             double* OLDSTATES, double* ALGEBRAIC, double* CONDVARS)
{
}