/*
   There are a total of 7 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 19 entries in the constant variable array.
 */
/*
 * VOI is time in component environment (second).
 * CONSTANTS[0] is L_1 in component contraction (mm).
 * CONSTANTS[1] is L_2 in component contraction (mm).
 * CONSTANTS[2] is L_3 in component contraction (mm).
 * CONSTANTS[3] is L_4 in component contraction (mm).
 * CONSTANTS[4] is f_c in component contraction (newton).
 * CONSTANTS[5] is v_max in component contraction (mm_per_second).
 * CONSTANTS[6] is curv in component contraction (dimensionless).
 * CONSTANTS[7] is k_1 in component contraction (newton).
 * CONSTANTS[8] is k_2 in component contraction (per_mm).
 * CONSTANTS[9] is F_1 in component contraction (newton).
 * CONSTANTS[10] is d_LSEC1 in component contraction (mm).
 * CONSTANTS[11] is k_sh in component contraction (dimensionless).
 * CONSTANTS[12] is L_m in component contraction (mm).
 * CONSTANTS[13] is F_im in component contraction (newton).
 * CONSTANTS[14] is tau in component contraction (second).
 * ALGEBRAIC[5] is v_cc in component contraction (mm_per_second).
 * ALGEBRAIC[6] is f_v in component contraction (dimensionless).
 * ALGEBRAIC[1] is f_L in component contraction (newton).
 * ALGEBRAIC[4] is f_sec in component contraction (newton).
 * ALGEBRAIC[0] is f_pec in component contraction (newton).
 * ALGEBRAIC[3] is delta_L_sec in component contraction (mm).
 * CONSTANTS[15] is delta_L_sec1 in component contraction (mm).
 * STATES[0] is delta_L_pec in component contraction (mm).
 * CONSTANTS[16] is k in component contraction (newton_per_mm).
 * ALGEBRAIC[2] is L_mtc in component contraction (mm).
 * CONSTANTS[17] is A in component contraction (dimensionless).
 * CONSTANTS[18] is L_mslack in component contraction (mm).
 * RATES[0] is d/dt delta_L_pec in component contraction (mm).
 */
void
initConsts(double* CONSTANTS, double* RATES, double *STATES)
{
CONSTANTS[0] = -23;
CONSTANTS[1] = -14;
CONSTANTS[2] = 2;
CONSTANTS[3] = 19;
CONSTANTS[4] = 0.49;
CONSTANTS[5] = -141;
CONSTANTS[6] = 5.8;
CONSTANTS[7] = 0.012;
CONSTANTS[8] = 0.317;
CONSTANTS[9] = 4.1;
CONSTANTS[10] = 4.1;
CONSTANTS[11] = 3.3;
CONSTANTS[12] = 0.3;
CONSTANTS[13] = 18.1;
CONSTANTS[14] = 0.006;
CONSTANTS[15] = 4.1;
STATES[0] = 0.2;
CONSTANTS[16] = 3.5;
CONSTANTS[17] = 1;
CONSTANTS[18] = 0.3;
}
void
computeRates(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[1] = (STATES[0]>=CONSTANTS[0]&&STATES[0]<=CONSTANTS[1] ?  (CONSTANTS[4]/(CONSTANTS[1] - CONSTANTS[0]))*(STATES[0] - CONSTANTS[0]) : STATES[0]>CONSTANTS[1]&&STATES[0]<=0.00000 ?  ((1.00000 - CONSTANTS[4])/- CONSTANTS[1])*(STATES[0] - CONSTANTS[1]) : STATES[0]>0.00000&&STATES[0]<=CONSTANTS[2] ? 1.00000 : STATES[0]>CONSTANTS[2]&&STATES[0]<=CONSTANTS[3] ?  (- 1.00000/(CONSTANTS[3] - CONSTANTS[2]))*(STATES[0] - CONSTANTS[2]) : 0.0/0.0);
ALGEBRAIC[2] = (VOI<=1.00000 ? 0.290000 : VOI>1.00000&&VOI<5.00000 ? 0.220000 : 0.190000);
ALGEBRAIC[3] = (ALGEBRAIC[2] - STATES[0]) - CONSTANTS[18];
ALGEBRAIC[4] = (ALGEBRAIC[3]>0.00000&&ALGEBRAIC[3]<CONSTANTS[15] ?  (CONSTANTS[9]/(exp(CONSTANTS[11]) - 1.00000))*(exp(( CONSTANTS[11]*ALGEBRAIC[3])/CONSTANTS[15]) - 1.00000) : ALGEBRAIC[3]<=CONSTANTS[15] ? CONSTANTS[9]+ CONSTANTS[16]*(ALGEBRAIC[3] - CONSTANTS[15]) : 0.0/0.0);
ALGEBRAIC[0] = (STATES[0]>0.00000 ?  CONSTANTS[7]*(exp( CONSTANTS[8]*STATES[0]) - 1.00000) : STATES[0]<=0.00000 ? 0.00000 : 0.0/0.0);
rootfind_0(VOI, CONSTANTS, RATES, STATES, ALGEBRAIC, pret);
RATES[0] = ALGEBRAIC[5];
}
void
computeVariables(double VOI, double* CONSTANTS, double* RATES, double* STATES, double* ALGEBRAIC)
{
ALGEBRAIC[1] = (STATES[0]>=CONSTANTS[0]&&STATES[0]<=CONSTANTS[1] ?  (CONSTANTS[4]/(CONSTANTS[1] - CONSTANTS[0]))*(STATES[0] - CONSTANTS[0]) : STATES[0]>CONSTANTS[1]&&STATES[0]<=0.00000 ?  ((1.00000 - CONSTANTS[4])/- CONSTANTS[1])*(STATES[0] - CONSTANTS[1]) : STATES[0]>0.00000&&STATES[0]<=CONSTANTS[2] ? 1.00000 : STATES[0]>CONSTANTS[2]&&STATES[0]<=CONSTANTS[3] ?  (- 1.00000/(CONSTANTS[3] - CONSTANTS[2]))*(STATES[0] - CONSTANTS[2]) : 0.0/0.0);
ALGEBRAIC[2] = (VOI<=1.00000 ? 0.290000 : VOI>1.00000&&VOI<5.00000 ? 0.220000 : 0.190000);
ALGEBRAIC[3] = (ALGEBRAIC[2] - STATES[0]) - CONSTANTS[18];
ALGEBRAIC[4] = (ALGEBRAIC[3]>0.00000&&ALGEBRAIC[3]<CONSTANTS[15] ?  (CONSTANTS[9]/(exp(CONSTANTS[11]) - 1.00000))*(exp(( CONSTANTS[11]*ALGEBRAIC[3])/CONSTANTS[15]) - 1.00000) : ALGEBRAIC[3]<=CONSTANTS[15] ? CONSTANTS[9]+ CONSTANTS[16]*(ALGEBRAIC[3] - CONSTANTS[15]) : 0.0/0.0);
ALGEBRAIC[0] = (STATES[0]>0.00000 ?  CONSTANTS[7]*(exp( CONSTANTS[8]*STATES[0]) - 1.00000) : STATES[0]<=0.00000 ? 0.00000 : 0.0/0.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] -  (1.00000/ALGEBRAIC[6])*((ALGEBRAIC[4] - ALGEBRAIC[0])/( CONSTANTS[17]*ALGEBRAIC[1]*CONSTANTS[13]));
  hx[1] = ALGEBRAIC[6] - (CONSTANTS[5] - ALGEBRAIC[5])/(CONSTANTS[5]+ ALGEBRAIC[5]*CONSTANTS[6]);
#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];
}