/*
   There are a total of 8 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 10 entries in the constant variable array.
 */
/*
 * VOI is time in component environment (minute).
 * CONSTANTS[0] is BFM in component muscle_O2_delivery (L_per_minute).
 * CONSTANTS[1] is OVA in component muscle_O2_delivery (mL_per_L).
 * CONSTANTS[2] is HM in component muscle_O2_delivery (dimensionless).
 * CONSTANTS[3] is AOM in component muscle_O2_delivery (dimensionless).
 * CONSTANTS[9] is O2ARTM in component M_O2_blood_supply (mL_per_minute).
 * ALGEBRAIC[4] is RMO in component delivery_of_O2_to_M_tissues (mL_per_minute).
 * ALGEBRAIC[5] is PVO in component M_venous_O2_content (mmHg).
 * CONSTANTS[4] is EXC in component parameter_values (dimensionless).
 * CONSTANTS[5] is EXCXP2 in component parameter_values (dimensionless).
 * ALGEBRAIC[6] is OVS in component M_venous_O2_content (dimensionless).
 * ALGEBRAIC[1] is PMO in component pressure_of_O2_in_M_tissue_cells (mmHg).
 * ALGEBRAIC[3] is MMO in component metabolic_O2_consumption_by_M_tissue (mL_per_minute).
 * ALGEBRAIC[2] is P2O in component metabolic_O2_consumption_by_M_tissue (mmHg).
 * CONSTANTS[6] is OMM in component parameter_values (mL_per_minute).
 * CONSTANTS[7] is PM5 in component parameter_values (per_mmHg).
 * ALGEBRAIC[0] is QOM in component volume_of_O2_in_M_tissue (mL).
 * ALGEBRAIC[7] is DO2M in component volume_of_O2_in_M_tissue (mL_per_minute).
 * STATES[0] is QOM1 in component volume_of_O2_in_M_tissue (mL).
 * CONSTANTS[8] is PK2 in component parameter_values (mmHg_per_mL).
 * RATES[0] is d/dt QOM1 in component volume_of_O2_in_M_tissue (mL).
 */
void
initConsts(double* CONSTANTS, double* RATES, double *STATES)
{
CONSTANTS[0] = 0.989949;
CONSTANTS[1] = 204.497;
CONSTANTS[2] = 40.0381;
CONSTANTS[3] = 1.00002;
CONSTANTS[4] = 1;
CONSTANTS[5] = 0.17;
CONSTANTS[6] = 57.1;
CONSTANTS[7] = 30;
STATES[0] = 48.0839;
CONSTANTS[8] = 0.79167;
CONSTANTS[9] =  CONSTANTS[1]*CONSTANTS[0];
}
void
computeRates(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[0] = (STATES[0]<0.000100000 ? 0.000100000 : STATES[0]);
ALGEBRAIC[1] =  CONSTANTS[8]*ALGEBRAIC[0];
rootfind_0(VOI, CONSTANTS, RATES, STATES, ALGEBRAIC, pret);
ALGEBRAIC[2] = (ALGEBRAIC[1]>38.0000 ? 38.0000 : ALGEBRAIC[1]);
ALGEBRAIC[3] =  CONSTANTS[3]*CONSTANTS[6]*CONSTANTS[4]*(1.00000 - pow(38.0001 - ALGEBRAIC[2], 3.00000)/54872.0);
ALGEBRAIC[7] = ALGEBRAIC[4] - ALGEBRAIC[3];
RATES[0] = ALGEBRAIC[7];
}
void
computeVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[0] = (STATES[0]<0.000100000 ? 0.000100000 : STATES[0]);
ALGEBRAIC[1] =  CONSTANTS[8]*ALGEBRAIC[0];
ALGEBRAIC[2] = (ALGEBRAIC[1]>38.0000 ? 38.0000 : ALGEBRAIC[1]);
ALGEBRAIC[3] =  CONSTANTS[3]*CONSTANTS[6]*CONSTANTS[4]*(1.00000 - pow(38.0001 - ALGEBRAIC[2], 3.00000)/54872.0);
ALGEBRAIC[7] = ALGEBRAIC[4] - ALGEBRAIC[3];
}

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[9] - ALGEBRAIC[4])/( CONSTANTS[2]*5.25000*CONSTANTS[0]);
  hx[1] = ALGEBRAIC[5] -  57.1400*ALGEBRAIC[6]*pow(CONSTANTS[4], CONSTANTS[5]);
  hx[2] = ALGEBRAIC[4] -  (ALGEBRAIC[5] - ALGEBRAIC[1])*CONSTANTS[7]*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] = {57.064,39.9701,0.699673};
  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];
}