/*
   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 29 entries in the constant variable array.
 */
/*
 * VOI is time in component environment (second).
 * ALGEBRAIC[8] is F_CE in component F_CE (newton).
 * ALGEBRAIC[0] is f_L_CE in component f_L_CE (newton).
 * ALGEBRAIC[5] is g_V_CE in component g_V_CE (dimensionless).
 * CONSTANTS[0] is a in component user_defined_constants (dimensionless).
 * CONSTANTS[1] is F_min in component f_L_CE (newton).
 * CONSTANTS[2] is F_max in component user_defined_constants (newton).
 * STATES[0] is L_CE in component L_CE (metre).
 * CONSTANTS[3] is L_CE_opt in component user_defined_constants (metre).
 * CONSTANTS[4] is W in component f_L_CE (dimensionless).
 * CONSTANTS[28] is lambda_a in component lambda_a (second_per_metre).
 * CONSTANTS[5] is V_max in component g_V_CE (metre_per_second).
 * ALGEBRAIC[6] is V_CE in component V_CE (dimensionless).
 * CONSTANTS[6] is A in component g_V_CE (dimensionless).
 * CONSTANTS[7] is g_max in component g_V_CE (dimensionless).
 * CONSTANTS[24] is d1 in component d1 (dimensionless).
 * CONSTANTS[26] is d2 in component d2 (dimensionless).
 * CONSTANTS[27] is d3 in component d3 (dimensionless).
 * CONSTANTS[8] is gamma in component g_V_CE (dimensionless).
 * CONSTANTS[9] is V_max in component d1 (metre_per_second).
 * CONSTANTS[10] is A in component d1 (dimensionless).
 * CONSTANTS[11] is g_max in component d1 (dimensionless).
 * CONSTANTS[12] is S in component d1 (metre_per_second).
 * CONSTANTS[13] is S in component d2 (metre_per_second).
 * CONSTANTS[14] is A in component d2 (dimensionless).
 * CONSTANTS[15] is V_max in component d2 (metre_per_second).
 * CONSTANTS[16] is gamma in component d2 (dimensionless).
 * CONSTANTS[17] is g_max in component d3 (dimensionless).
 * CONSTANTS[18] is gamma in component d3 (dimensionless).
 * ALGEBRAIC[4] is F_SEE in component F_SEE (newton).
 * CONSTANTS[19] is k_SEE in component F_SEE (newton_per_metre2).
 * ALGEBRAIC[3] is L_SEE in component L_SEE (metre).
 * CONSTANTS[20] is L_slack in component F_SEE (metre).
 * ALGEBRAIC[1] is F_PEE in component F_PEE (newton).
 * CONSTANTS[25] is k_PEE in component k_PEE (newton_per_metre2).
 * CONSTANTS[21] is L_slack in component F_PEE (metre).
 * CONSTANTS[22] is W in component k_PEE (dimensionless).
 * CONSTANTS[23] is L_CE_opt in component k_PEE (metre).
 * ALGEBRAIC[2] is L_m in component L_m (metre).
 * ALGEBRAIC[7] is F_m in component F_m (newton).
 * RATES[0] is d/dt L_CE in component L_CE (metre).
 */
void
initConsts(double* CONSTANTS, double* RATES, double *STATES)
{
CONSTANTS[0] = 0.8;
CONSTANTS[1] = 10;
CONSTANTS[2] = 7000;
STATES[0] = 0.038;
CONSTANTS[3] = 0.093;
CONSTANTS[4] = 0.63;
CONSTANTS[5] = 0.93;
CONSTANTS[6] = 0.25;
CONSTANTS[7] = 1.5;
CONSTANTS[8] = 5.67;
CONSTANTS[9] = 0.93;
CONSTANTS[10] = 0.25;
CONSTANTS[11] = 1.5;
CONSTANTS[12] = 2;
CONSTANTS[13] = 2;
CONSTANTS[14] = 0.25;
CONSTANTS[15] = 0.93;
CONSTANTS[16] = 5.67;
CONSTANTS[17] = 1.5;
CONSTANTS[18] = 5.67;
CONSTANTS[19] = 1000000;
CONSTANTS[20] = 0.0025;
CONSTANTS[21] = 0.0025;
CONSTANTS[22] = 0.63;
CONSTANTS[23] = 0.01;
CONSTANTS[24] = ( CONSTANTS[9]*CONSTANTS[10]*(CONSTANTS[11] - 1.00000))/( CONSTANTS[12]*(CONSTANTS[10]+1.00000));
CONSTANTS[25] = CONSTANTS[2]/pow( CONSTANTS[22]*CONSTANTS[23], 2.00000);
CONSTANTS[26] = ( CONSTANTS[13]*(CONSTANTS[14]+1.00000))/( CONSTANTS[15]*CONSTANTS[14]*pow(CONSTANTS[16]+1.00000, 2.00000));
CONSTANTS[27] = ( (CONSTANTS[17] - 1.00000)*pow(CONSTANTS[18], 2.00000))/pow(CONSTANTS[18]+1.00000, 2.00000)+1.00000;
CONSTANTS[28] =  1.00000*((1.00000 - exp( - 3.82000*CONSTANTS[0]))+ CONSTANTS[0]*exp(- 3.82000));
}
void
computeRates(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[0] = ( CONSTANTS[2]*( 1.00000*(1.00000 - STATES[0]) - pow(CONSTANTS[3], 2.00000)))/( pow(CONSTANTS[4], 2.00000)*pow(CONSTANTS[3], 2.00000));
ALGEBRAIC[2] = (VOI<=1.00000 ? 0.0380000 : VOI>1.00000&&VOI<2.00000 ? 0.0380000+ 0.00200000*(VOI - 1.00000) : 0.0400000);
ALGEBRAIC[3] = ALGEBRAIC[2] - STATES[0];
ALGEBRAIC[4] = (ALGEBRAIC[3]<=CONSTANTS[20] ? 0.00000 :  CONSTANTS[19]*pow(ALGEBRAIC[3] - CONSTANTS[20], 2.00000));
ALGEBRAIC[1] = (STATES[0]<=CONSTANTS[21] ? 0.00000 :  CONSTANTS[25]*pow(STATES[0] - CONSTANTS[21], 2.00000));
rootfind_0(VOI, CONSTANTS, RATES, STATES, ALGEBRAIC, pret);
RATES[0] =  1.00000*ALGEBRAIC[6];
}
void
computeVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[0] = ( CONSTANTS[2]*( 1.00000*(1.00000 - STATES[0]) - pow(CONSTANTS[3], 2.00000)))/( pow(CONSTANTS[4], 2.00000)*pow(CONSTANTS[3], 2.00000));
ALGEBRAIC[2] = (VOI<=1.00000 ? 0.0380000 : VOI>1.00000&&VOI<2.00000 ? 0.0380000+ 0.00200000*(VOI - 1.00000) : 0.0400000);
ALGEBRAIC[3] = ALGEBRAIC[2] - STATES[0];
ALGEBRAIC[4] = (ALGEBRAIC[3]<=CONSTANTS[20] ? 0.00000 :  CONSTANTS[19]*pow(ALGEBRAIC[3] - CONSTANTS[20], 2.00000));
ALGEBRAIC[1] = (STATES[0]<=CONSTANTS[21] ? 0.00000 :  CONSTANTS[25]*pow(STATES[0] - CONSTANTS[21], 2.00000));
ALGEBRAIC[7] = ALGEBRAIC[4];
ALGEBRAIC[8] =  ALGEBRAIC[0]*ALGEBRAIC[5]*CONSTANTS[0];
}

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[5] = p[0];
  ALGEBRAIC[6] = p[1];
  hx[0] = ALGEBRAIC[5] - (ALGEBRAIC[6]<=0.00000 ? ( CONSTANTS[28]*CONSTANTS[5]+ALGEBRAIC[6])/( CONSTANTS[28]*CONSTANTS[5] - ALGEBRAIC[6]/CONSTANTS[6]) : 0.00000<ALGEBRAIC[6]&&ALGEBRAIC[6]<= CONSTANTS[8]*CONSTANTS[24] ? ( CONSTANTS[7]*ALGEBRAIC[6]+CONSTANTS[24])/(ALGEBRAIC[6]+CONSTANTS[24]) : ALGEBRAIC[6]> CONSTANTS[8]*CONSTANTS[24] ? CONSTANTS[27]+ CONSTANTS[26]*ALGEBRAIC[6] : 0.0/0.0);
  hx[1] = ALGEBRAIC[6] -  1.00000*(( (1.00000/ALGEBRAIC[5])*( ALGEBRAIC[4]*(ALGEBRAIC[2] - STATES[0]) -  ALGEBRAIC[1]*STATES[0]))/( CONSTANTS[0]*ALGEBRAIC[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[2] = {0.1,0.1};
  double bp[2], work[LM_DIF_WORKSZ(2, 2)];
  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, 2, &rfi);
  ALGEBRAIC[5] = p[0];
  ALGEBRAIC[6] = p[1];
}