/*
   There are a total of 9 entries in the algebraic variable array.
   There are a total of 1 entries in each of the rate and state variable arrays.
   There are a total of 6 entries in the constant variable array.
 */
/*
 * VOI is time in component environment (minute).
 * CONSTANTS[0] is BFN in component non_muscle_O2_delivery (L_per_minute).
 * CONSTANTS[1] is OVA in component non_muscle_O2_delivery (mL_per_L).
 * CONSTANTS[2] is HM in component non_muscle_O2_delivery (dimensionless).
 * CONSTANTS[3] is AOM in component non_muscle_O2_delivery (dimensionless).
 * CONSTANTS[5] is O2ARTN in component NM_O2_blood_supply (mL_per_minute).
 * ALGEBRAIC[4] is DOB in component delivery_of_O2_to_NM_tissues (mL_per_minute).
 * ALGEBRAIC[5] is POV in component NM_venous_O2_content (mmHg).
 * ALGEBRAIC[6] is OSV in component NM_venous_O2_content (dimensionless).
 * ALGEBRAIC[1] is POT in component pressure_of_O2_in_NM_tissue_cells (mmHg).
 * ALGEBRAIC[3] is MO2 in component O2_consumption_by_NM_tissue (mL_per_minute).
 * CONSTANTS[4] is O2M in component parameter_values (mL_per_minute).
 * ALGEBRAIC[2] is P1O in component O2_consumption_by_NM_tissue (mmHg).
 * ALGEBRAIC[0] is QO2 in component volume_of_O2_in_NM_tissue (mL).
 * ALGEBRAIC[8] is DO2N in component volume_of_O2_in_NM_tissue (mL_per_minute).
 * ALGEBRAIC[7] is DO2N1 in component volume_of_O2_in_NM_tissue (mL_per_minute).
 * STATES[0] is QO2T in component volume_of_O2_in_NM_tissue (mL).
 * RATES[0] is d/dt QO2T in component volume_of_O2_in_NM_tissue (mL).
 */
void
initConsts(double* CONSTANTS, double* RATES, double *STATES)
{
CONSTANTS[0] = 2.79521;
CONSTANTS[1] = 204.497;
CONSTANTS[2] = 40.0381;
CONSTANTS[3] = 1.00002;
CONSTANTS[4] = 164;
STATES[0] = 72.2362;
CONSTANTS[5] =  CONSTANTS[1]*CONSTANTS[0];
}
void
computeRates(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[0] = (STATES[0]<0.00000 ? 0.00000 : STATES[0]);
ALGEBRAIC[1] =  ALGEBRAIC[0]*0.486110;
rootfind_0(VOI, CONSTANTS, RATES, STATES, ALGEBRAIC, pret);
ALGEBRAIC[2] = (ALGEBRAIC[1]>35.0000 ? 35.0000 : ALGEBRAIC[1]);
ALGEBRAIC[3] =  CONSTANTS[3]*CONSTANTS[4]*(1.00000 - pow(35.0001 - ALGEBRAIC[2], 3.00000)/42875.0);
ALGEBRAIC[7] = ALGEBRAIC[4] - ALGEBRAIC[3];
ALGEBRAIC[8] = (ALGEBRAIC[0]<6.00000&&ALGEBRAIC[7]<0.00000 ?  ALGEBRAIC[7]*0.100000 : ALGEBRAIC[7]);
RATES[0] = ALGEBRAIC[8];
}
void
computeVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[0] = (STATES[0]<0.00000 ? 0.00000 : STATES[0]);
ALGEBRAIC[1] =  ALGEBRAIC[0]*0.486110;
ALGEBRAIC[2] = (ALGEBRAIC[1]>35.0000 ? 35.0000 : ALGEBRAIC[1]);
ALGEBRAIC[3] =  CONSTANTS[3]*CONSTANTS[4]*(1.00000 - pow(35.0001 - ALGEBRAIC[2], 3.00000)/42875.0);
ALGEBRAIC[7] = ALGEBRAIC[4] - ALGEBRAIC[3];
ALGEBRAIC[8] = (ALGEBRAIC[0]<6.00000&&ALGEBRAIC[7]<0.00000 ?  ALGEBRAIC[7]*0.100000 : ALGEBRAIC[7]);
}

void objfunc_0(double *p, double *hx, int m, int n, void *adata)
{
  struct rootfind_info* rfi = (struct rootfind_info*)adata;
#define VOI rfi->aVOI
#define CONSTANTS rfi->aCONSTANTS
#define RATES rfi->aRATES
#define STATES rfi->aSTATES
#define ALGEBRAIC rfi->aALGEBRAIC
#define pret rfi->aPRET
  ALGEBRAIC[4] = p[0];
  ALGEBRAIC[5] = p[1];
  ALGEBRAIC[6] = p[2];
  hx[0] = ALGEBRAIC[6] - (CONSTANTS[5] - ALGEBRAIC[4])/( CONSTANTS[2]*5.25000*CONSTANTS[0]);
  hx[1] = ALGEBRAIC[5] -  ALGEBRAIC[6]*57.1400;
  hx[2] = ALGEBRAIC[4] -  (ALGEBRAIC[5] - ALGEBRAIC[1])*12.8570*CONSTANTS[0];
#undef VOI
#undef CONSTANTS
#undef RATES
#undef STATES
#undef ALGEBRAIC
#undef pret
}
void rootfind_0(double VOI, double* CONSTANTS, double* RATES,
double* STATES, double* ALGEBRAIC, int* pret)
{
  static double p[3] = {0.1,0.1,0.1};
  double bp[3], work[LM_DIF_WORKSZ(3, 3)];
  struct rootfind_info rfi;
  rfi.aVOI = VOI;
  rfi.aCONSTANTS = CONSTANTS;
  rfi.aRATES = RATES;
  rfi.aSTATES = STATES;
  rfi.aALGEBRAIC = ALGEBRAIC;
  rfi.aPRET = pret;
  do_levmar(objfunc_0, p, bp, work, pret, 3, &rfi);
  ALGEBRAIC[4] = p[0];
  ALGEBRAIC[5] = p[1];
  ALGEBRAIC[6] = p[2];
}