/*
   There are a total of 12 entries in the algebraic variable array.
   There are a total of 3 entries in each of the rate and state variable arrays.
   There are a total of 19 entries in the constant variable array.
 */
/*
 * VOI is time in component environment (dimensionless).
 * STATES[0] is C in component C (dimensionless).
 * CONSTANTS[0] is gamma in component C (dimensionless).
 * ALGEBRAIC[5] is Vs in component Vs (dimensionless).
 * ALGEBRAIC[0] is V1 in component V1 (dimensionless).
 * ALGEBRAIC[10] is V2 in component V2 (dimensionless).
 * ALGEBRAIC[6] is Vd in component Vd (dimensionless).
 * STATES[1] is K in component K (dimensionless).
 * STATES[2] is RP in component RP (dimensionless).
 * ALGEBRAIC[11] is V3 in component V3 (dimensionless).
 * ALGEBRAIC[1] is V4 in component V4 (dimensionless).
 * ALGEBRAIC[2] is E in component E (dimensionless).
 * CONSTANTS[1] is sigma in component E (dimensionless).
 * ALGEBRAIC[3] is RE in component RE (dimensionless).
 * ALGEBRAIC[7] is KP in component KP (dimensionless).
 * ALGEBRAIC[8] is KPI in component KPI (dimensionless).
 * CONSTANTS[2] is thetaI in component KPI (dimensionless).
 * ALGEBRAIC[9] is I in component I (dimensionless).
 * CONSTANTS[3] is lambda in component I (dimensionless).
 * CONSTANTS[4] is thetaE in component RE (dimensionless).
 * ALGEBRAIC[4] is R in component R (dimensionless).
 * CONSTANTS[5] is VCs in component Vs (dimensionless).
 * CONSTANTS[6] is Vsm in component Vs (dimensionless).
 * CONSTANTS[7] is KsE in component Vs (dimensionless).
 * CONSTANTS[8] is V1m in component V1 (dimensionless).
 * CONSTANTS[9] is K1C in component V1 (dimensionless).
 * CONSTANTS[10] is K1 in component V1 (dimensionless).
 * CONSTANTS[11] is V2m in component V2 (dimensionless).
 * CONSTANTS[12] is K2 in component V2 (dimensionless).
 * CONSTANTS[13] is V3m in component V3 (dimensionless).
 * CONSTANTS[14] is K3 in component V3 (dimensionless).
 * CONSTANTS[15] is V4m in component V4 (dimensionless).
 * CONSTANTS[16] is K4 in component V4 (dimensionless).
 * CONSTANTS[17] is VdEm in component Vd (dimensionless).
 * CONSTANTS[18] is KdC in component Vd (dimensionless).
 * RATES[0] is d/dt C in component C (dimensionless).
 * RATES[1] is d/dt K in component K (dimensionless).
 * RATES[2] is d/dt RP in component RP (dimensionless).
 */
void
initConsts(double* CONSTANTS, double* RATES, double *STATES)
{
STATES[0] = 0.38;
CONSTANTS[0] = 1.0;
STATES[1] = 0.1;
STATES[2] = 1.0;
CONSTANTS[1] = 10;
CONSTANTS[2] = 1.0;
CONSTANTS[3] = 1.0;
CONSTANTS[4] = 0.01;
CONSTANTS[5] = 1.0;
CONSTANTS[6] = 1.0;
CONSTANTS[7] = 0.1;
CONSTANTS[8] = 50.0;
CONSTANTS[9] = 0.1;
CONSTANTS[10] = 0.0001;
CONSTANTS[11] = 40;
CONSTANTS[12] = 0.0001;
CONSTANTS[13] = 3000;
CONSTANTS[14] = 0.0001;
CONSTANTS[15] = 3.0;
CONSTANTS[16] = 0.0001;
CONSTANTS[17] = 1000.0;
CONSTANTS[18] = 0.005;
}
void
computeRates(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
rootfind_0(VOI, CONSTANTS, RATES, STATES, ALGEBRAIC, pret);
ALGEBRAIC[5] = CONSTANTS[5]+ CONSTANTS[6]*(ALGEBRAIC[2]/(CONSTANTS[7]+ALGEBRAIC[2]));
ALGEBRAIC[0] =  CONSTANTS[8]*(STATES[0]/(CONSTANTS[9]+STATES[0]))*(STATES[1]/(CONSTANTS[10]+STATES[1]));
rootfind_1(VOI, CONSTANTS, RATES, STATES, ALGEBRAIC, pret);
ALGEBRAIC[10] =  CONSTANTS[11]*(ALGEBRAIC[7]/(CONSTANTS[12]+ALGEBRAIC[7]));
ALGEBRAIC[6] = STATES[0]+ CONSTANTS[17]*ALGEBRAIC[2]*(STATES[0]/(CONSTANTS[18]+STATES[0]));
RATES[0] = (ALGEBRAIC[5]+ CONSTANTS[0]*ALGEBRAIC[10]) - ( CONSTANTS[0]*ALGEBRAIC[0]+ALGEBRAIC[6]);
RATES[1] = ALGEBRAIC[10] - ALGEBRAIC[0];
ALGEBRAIC[11] =  CONSTANTS[13]*ALGEBRAIC[7]*(ALGEBRAIC[3]/(CONSTANTS[14]+ALGEBRAIC[3]));
ALGEBRAIC[1] =  CONSTANTS[15]*(STATES[2]/(CONSTANTS[16]+STATES[2]));
RATES[2] = ALGEBRAIC[11] - ALGEBRAIC[1];
}
void
computeVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[5] = CONSTANTS[5]+ CONSTANTS[6]*(ALGEBRAIC[2]/(CONSTANTS[7]+ALGEBRAIC[2]));
ALGEBRAIC[0] =  CONSTANTS[8]*(STATES[0]/(CONSTANTS[9]+STATES[0]))*(STATES[1]/(CONSTANTS[10]+STATES[1]));
ALGEBRAIC[10] =  CONSTANTS[11]*(ALGEBRAIC[7]/(CONSTANTS[12]+ALGEBRAIC[7]));
ALGEBRAIC[6] = STATES[0]+ CONSTANTS[17]*ALGEBRAIC[2]*(STATES[0]/(CONSTANTS[18]+STATES[0]));
ALGEBRAIC[11] =  CONSTANTS[13]*ALGEBRAIC[7]*(ALGEBRAIC[3]/(CONSTANTS[14]+ALGEBRAIC[3]));
ALGEBRAIC[1] =  CONSTANTS[15]*(STATES[2]/(CONSTANTS[16]+STATES[2]));
}

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[2] = p[0];
  ALGEBRAIC[3] = p[1];
  ALGEBRAIC[4] = p[2];
  hx[0] = ALGEBRAIC[2] - (1.00000 -  CONSTANTS[1]*ALGEBRAIC[3]);
  hx[1] = ALGEBRAIC[3] -  CONSTANTS[4]*ALGEBRAIC[4]*ALGEBRAIC[2];
  hx[2] = ALGEBRAIC[4] - (1.00000 - (STATES[2]+ALGEBRAIC[3]));
#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[2] = p[0];
  ALGEBRAIC[3] = p[1];
  ALGEBRAIC[4] = p[2];
}

void objfunc_1(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[7] = p[0];
  ALGEBRAIC[8] = p[1];
  ALGEBRAIC[9] = p[2];
  hx[0] = ALGEBRAIC[7] - (1.00000 - (ALGEBRAIC[8]+STATES[1]));
  hx[1] = ALGEBRAIC[8] -  CONSTANTS[2]*ALGEBRAIC[7]*ALGEBRAIC[9];
  hx[2] = ALGEBRAIC[9] - (1.00000 -  CONSTANTS[3]*ALGEBRAIC[8]);
#undef VOI
#undef CONSTANTS
#undef RATES
#undef STATES
#undef ALGEBRAIC
#undef pret
}
void rootfind_1(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_1, p, bp, work, pret, 3, &rfi);
  ALGEBRAIC[7] = p[0];
  ALGEBRAIC[8] = p[1];
  ALGEBRAIC[9] = p[2];
}