# Size of variable arrays:
sizeAlgebraic = 115
sizeStates = 39
sizeConstants = 125
from math import *
from numpy import *

def createLegends():
    legend_states = [""] * sizeStates
    legend_rates = [""] * sizeStates
    legend_algebraic = [""] * sizeAlgebraic
    legend_voi = ""
    legend_constants = [""] * sizeConstants
    legend_algebraic[22] = "I_app in component Protocol (uA_per_uF)"
    legend_voi = "time in component environment (msec)"
    legend_algebraic[9] = "V_clamp in component Protocol (mV)"
    legend_states[0] = "V_m in component membrane_potential (mV)"
    legend_constants[0] = "R_clamp in component Protocol (uF_per_mS)"
    legend_constants[1] = "epi in component parameters (dimensionless)"
    legend_constants[2] = "R in component parameters (joule_per_kelvin_per_kilomole)"
    legend_constants[96] = "Frdy in component parameters (coulomb_per_mole)"
    legend_constants[106] = "Temp in component parameters (kelvin)"
    legend_constants[110] = "FoRT in component parameters (per_mV)"
    legend_constants[3] = "Cmem in component parameters (farad)"
    legend_constants[111] = "Qpow in component parameters (dimensionless)"
    legend_constants[4] = "cellLength in component parameters (um)"
    legend_constants[97] = "cellRadius in component parameters (um)"
    legend_constants[5] = "junctionLength in component parameters (um)"
    legend_constants[98] = "junctionRadius in component parameters (um)"
    legend_constants[6] = "distSLcyto in component parameters (um)"
    legend_constants[7] = "distJuncSL in component parameters (um)"
    legend_constants[8] = "DcaJuncSL in component parameters (cm2_per_sec)"
    legend_constants[9] = "DcaSLcyto in component parameters (cm2_per_sec)"
    legend_constants[10] = "DnaJuncSL in component parameters (cm2_per_sec)"
    legend_constants[11] = "DnaSLcyto in component parameters (cm2_per_sec)"
    legend_constants[107] = "Vcell in component parameters (liter)"
    legend_constants[112] = "Vmyo in component parameters (liter)"
    legend_constants[113] = "Vsr in component parameters (liter)"
    legend_constants[114] = "Vsl in component parameters (liter)"
    legend_constants[115] = "Vjunc in component parameters (liter)"
    legend_constants[108] = "SAjunc in component parameters (um2)"
    legend_constants[109] = "SAsl in component parameters (um2)"
    legend_constants[12] = "J_ca_juncsl in component parameters (liters_per_msec)"
    legend_constants[13] = "J_ca_slmyo in component parameters (liters_per_msec)"
    legend_constants[14] = "J_na_juncsl in component parameters (liters_per_msec)"
    legend_constants[15] = "J_na_slmyo in component parameters (liters_per_msec)"
    legend_constants[16] = "Fjunc in component parameters (dimensionless)"
    legend_constants[99] = "Fsl in component parameters (dimensionless)"
    legend_constants[17] = "Fjunc_CaL in component parameters (dimensionless)"
    legend_constants[100] = "Fsl_CaL in component parameters (dimensionless)"
    legend_constants[116] = "Cli in component parameters (mM)"
    legend_constants[117] = "Clo in component parameters (mM)"
    legend_constants[18] = "Ko in component parameters (mM)"
    legend_constants[19] = "Nao in component parameters (mM)"
    legend_constants[20] = "Cao in component parameters (mM)"
    legend_constants[21] = "Mgi in component parameters (mM)"
    legend_algebraic[26] = "ena_junc in component parameters (mV)"
    legend_algebraic[30] = "ena_sl in component parameters (mV)"
    legend_algebraic[32] = "ek in component parameters (mV)"
    legend_algebraic[33] = "eca_junc in component parameters (mV)"
    legend_algebraic[34] = "eca_sl in component parameters (mV)"
    legend_constants[123] = "ecl in component parameters (mV)"
    legend_constants[22] = "GNa in component parameters (mS_per_uF)"
    legend_constants[23] = "GNaB in component parameters (mS_per_uF)"
    legend_constants[24] = "IbarNaK in component parameters (uA_per_uF)"
    legend_constants[25] = "KmNaip in component parameters (mM)"
    legend_constants[26] = "KmKo in component parameters (mM)"
    legend_constants[27] = "Q10NaK in component parameters (dimensionless)"
    legend_constants[28] = "Q10KmNai in component parameters (dimensionless)"
    legend_constants[29] = "pNaK in component parameters (dimensionless)"
    legend_constants[30] = "gkp in component parameters (mS_per_uF)"
    legend_constants[31] = "GClCa in component parameters (mS_per_uF)"
    legend_constants[32] = "GClB in component parameters (mS_per_uF)"
    legend_constants[33] = "KdClCa in component parameters (mM)"
    legend_constants[34] = "pNa in component parameters (cm_per_sec)"
    legend_constants[35] = "pCa in component parameters (cm_per_sec)"
    legend_constants[36] = "pK in component parameters (cm_per_sec)"
    legend_constants[37] = "Q10CaL in component parameters (dimensionless)"
    legend_constants[38] = "IbarNCX in component parameters (uA_per_uF)"
    legend_constants[39] = "KmCai in component parameters (mM)"
    legend_constants[40] = "KmCao in component parameters (mM)"
    legend_constants[41] = "KmNai in component parameters (mM)"
    legend_constants[42] = "KmNao in component parameters (mM)"
    legend_constants[43] = "ksat in component parameters (dimensionless)"
    legend_constants[44] = "nu in component parameters (dimensionless)"
    legend_constants[45] = "Kdact in component parameters (mM)"
    legend_constants[46] = "Q10NCX in component parameters (dimensionless)"
    legend_constants[47] = "IbarSLCaP in component parameters (uA_per_uF)"
    legend_constants[48] = "KmPCa in component parameters (mM)"
    legend_constants[49] = "GCaB in component parameters (mS_per_uF)"
    legend_constants[50] = "Q10SLCaP in component parameters (dimensionless)"
    legend_constants[51] = "Q10SRCaP in component parameters (dimensionless)"
    legend_constants[52] = "Vmax_SRCaP in component parameters (mM_per_msec)"
    legend_constants[53] = "Kmf in component parameters (mM)"
    legend_constants[54] = "Kmr in component parameters (mM)"
    legend_constants[55] = "hillSRCaP in component parameters (dimensionless)"
    legend_constants[56] = "ks in component parameters (per_msec)"
    legend_constants[57] = "koCa in component parameters (per_mM2_per_msec)"
    legend_constants[58] = "kom in component parameters (per_msec)"
    legend_constants[105] = "kiCa in component parameters (per_mM_per_msec)"
    legend_constants[59] = "kim in component parameters (per_msec)"
    legend_constants[60] = "ec50SR in component parameters (mM)"
    legend_constants[61] = "Bmax_Naj in component parameters (mM)"
    legend_constants[62] = "Bmax_Nasl in component parameters (mM)"
    legend_constants[63] = "koff_na in component parameters (per_msec)"
    legend_constants[64] = "kon_na in component parameters (per_mM_per_msec)"
    legend_constants[65] = "Bmax_TnClow in component parameters (mM)"
    legend_constants[66] = "koff_tncl in component parameters (per_msec)"
    legend_constants[67] = "kon_tncl in component parameters (per_mM_per_msec)"
    legend_constants[68] = "Bmax_TnChigh in component parameters (mM)"
    legend_constants[69] = "koff_tnchca in component parameters (per_msec)"
    legend_constants[70] = "kon_tnchca in component parameters (per_mM_per_msec)"
    legend_constants[71] = "koff_tnchmg in component parameters (per_msec)"
    legend_constants[72] = "kon_tnchmg in component parameters (per_mM_per_msec)"
    legend_constants[73] = "Bmax_CaM in component parameters (mM)"
    legend_constants[74] = "koff_cam in component parameters (per_msec)"
    legend_constants[75] = "kon_cam in component parameters (per_mM_per_msec)"
    legend_constants[76] = "Bmax_myosin in component parameters (mM)"
    legend_constants[77] = "koff_myoca in component parameters (per_msec)"
    legend_constants[78] = "kon_myoca in component parameters (per_mM_per_msec)"
    legend_constants[79] = "koff_myomg in component parameters (per_msec)"
    legend_constants[80] = "kon_myomg in component parameters (per_mM_per_msec)"
    legend_constants[81] = "Bmax_SR in component parameters (mM)"
    legend_constants[82] = "koff_sr in component parameters (per_msec)"
    legend_constants[83] = "kon_sr in component parameters (per_mM_per_msec)"
    legend_constants[118] = "Bmax_SLlowsl in component parameters (mM)"
    legend_constants[119] = "Bmax_SLlowj in component parameters (mM)"
    legend_constants[84] = "koff_sll in component parameters (per_msec)"
    legend_constants[85] = "kon_sll in component parameters (per_mM_per_msec)"
    legend_constants[120] = "Bmax_SLhighsl in component parameters (mM)"
    legend_constants[121] = "Bmax_SLhighj in component parameters (mM)"
    legend_constants[86] = "koff_slh in component parameters (per_msec)"
    legend_constants[87] = "kon_slh in component parameters (per_mM_per_msec)"
    legend_constants[122] = "Bmax_Csqn in component parameters (mM)"
    legend_constants[88] = "koff_csqn in component parameters (per_msec)"
    legend_constants[89] = "kon_csqn in component parameters (per_mM_per_msec)"
    legend_states[1] = "Na_j in component Na_Concentrations (mM)"
    legend_states[2] = "Na_sl in component Na_Concentrations (mM)"
    legend_states[3] = "K_i in component K_Concentration (mM)"
    legend_states[4] = "Ca_j in component Ca_Concentrations (mM)"
    legend_states[5] = "Ca_sl in component Ca_Concentrations (mM)"
    legend_algebraic[0] = "mss in component I_Na (dimensionless)"
    legend_algebraic[11] = "taum in component I_Na (msec)"
    legend_algebraic[1] = "ah in component I_Na (dimensionless)"
    legend_algebraic[12] = "bh in component I_Na (dimensionless)"
    legend_algebraic[24] = "tauh in component I_Na (msec)"
    legend_algebraic[28] = "hss in component I_Na (dimensionless)"
    legend_algebraic[2] = "aj in component I_Na (dimensionless)"
    legend_algebraic[13] = "bj in component I_Na (dimensionless)"
    legend_algebraic[25] = "tauj in component I_Na (msec)"
    legend_algebraic[29] = "jss in component I_Na (dimensionless)"
    legend_states[6] = "m in component I_Na (dimensionless)"
    legend_states[7] = "h in component I_Na (dimensionless)"
    legend_states[8] = "j in component I_Na (dimensionless)"
    legend_algebraic[35] = "I_Na_junc in component I_Na (uA_per_uF)"
    legend_algebraic[36] = "I_Na_sl in component I_Na (uA_per_uF)"
    legend_algebraic[37] = "I_Na in component I_Na (uA_per_uF)"
    legend_algebraic[38] = "I_nabk_junc in component I_NaBK (uA_per_uF)"
    legend_algebraic[39] = "I_nabk_sl in component I_NaBK (uA_per_uF)"
    legend_algebraic[40] = "I_nabk in component I_NaBK (uA_per_uF)"
    legend_constants[101] = "sigma in component I_NaK (dimensionless)"
    legend_algebraic[41] = "fnak in component I_NaK (dimensionless)"
    legend_algebraic[42] = "I_nak_junc in component I_NaK (uA_per_uF)"
    legend_algebraic[43] = "I_nak_sl in component I_NaK (uA_per_uF)"
    legend_algebraic[44] = "I_nak in component I_NaK (uA_per_uF)"
    legend_constants[102] = "gkr in component I_Kr (mS_per_uF)"
    legend_algebraic[3] = "xrss in component I_Kr (dimensionless)"
    legend_algebraic[14] = "tauxr in component I_Kr (msec)"
    legend_states[9] = "x_kr in component I_Kr (dimensionless)"
    legend_algebraic[45] = "rkr in component I_Kr (dimensionless)"
    legend_algebraic[46] = "I_kr in component I_Kr (uA_per_uF)"
    legend_algebraic[47] = "kp_kp in component I_Kp (dimensionless)"
    legend_algebraic[48] = "I_kp_junc in component I_Kp (uA_per_uF)"
    legend_algebraic[49] = "I_kp_sl in component I_Kp (uA_per_uF)"
    legend_algebraic[50] = "I_kp in component I_Kp (uA_per_uF)"
    legend_algebraic[51] = "eks in component I_Ks (mV)"
    legend_constants[90] = "gks_junc in component I_Ks (mS_per_uF)"
    legend_constants[91] = "gks_sl in component I_Ks (mS_per_uF)"
    legend_algebraic[4] = "xsss in component I_Ks (dimensionless)"
    legend_algebraic[15] = "tauxs in component I_Ks (msec)"
    legend_states[10] = "x_ks in component I_Ks (dimensionless)"
    legend_algebraic[52] = "I_ks_junc in component I_Ks (uA_per_uF)"
    legend_algebraic[53] = "I_ks_sl in component I_Ks (uA_per_uF)"
    legend_algebraic[54] = "I_ks in component I_Ks (uA_per_uF)"
    legend_states[11] = "Na_i in component Na_Concentrations (mM)"
    legend_constants[103] = "GtoSlow in component I_to (mS_per_uF)"
    legend_constants[104] = "GtoFast in component I_to (mS_per_uF)"
    legend_algebraic[5] = "xtoss in component I_to (dimensionless)"
    legend_algebraic[6] = "ytoss in component I_to (dimensionless)"
    legend_algebraic[16] = "tauxtos in component I_to (msec)"
    legend_algebraic[17] = "tauytos in component I_to (msec)"
    legend_states[12] = "x_to_s in component I_to (dimensionless)"
    legend_states[13] = "y_to_s in component I_to (dimensionless)"
    legend_algebraic[55] = "I_tos in component I_to (uA_per_uF)"
    legend_algebraic[18] = "tauxtof in component I_to (msec)"
    legend_algebraic[19] = "tauytof in component I_to (msec)"
    legend_states[14] = "x_to_f in component I_to (dimensionless)"
    legend_states[15] = "y_to_f in component I_to (dimensionless)"
    legend_algebraic[56] = "I_tof in component I_to (uA_per_uF)"
    legend_algebraic[57] = "I_to in component I_to (uA_per_uF)"
    legend_algebraic[58] = "aki in component I_Ki (dimensionless)"
    legend_algebraic[59] = "bki in component I_Ki (dimensionless)"
    legend_algebraic[60] = "kiss in component I_Ki (dimensionless)"
    legend_algebraic[61] = "I_ki in component I_Ki (uA_per_uF)"
    legend_algebraic[62] = "I_ClCa_junc in component I_ClCa (uA_per_uF)"
    legend_algebraic[63] = "I_ClCa_sl in component I_ClCa (uA_per_uF)"
    legend_algebraic[64] = "I_ClCa in component I_ClCa (uA_per_uF)"
    legend_algebraic[65] = "I_Clbk in component I_ClCa (uA_per_uF)"
    legend_algebraic[7] = "fss in component I_Ca (dimensionless)"
    legend_algebraic[8] = "dss in component I_Ca (dimensionless)"
    legend_algebraic[20] = "taud in component I_Ca (msec)"
    legend_algebraic[21] = "tauf in component I_Ca (msec)"
    legend_states[16] = "d in component I_Ca (dimensionless)"
    legend_states[17] = "f in component I_Ca (dimensionless)"
    legend_states[18] = "f_Ca_Bj in component I_Ca (dimensionless)"
    legend_states[19] = "f_Ca_Bsl in component I_Ca (dimensionless)"
    legend_constants[92] = "fcaCaMSL in component I_Ca (dimensionless)"
    legend_constants[93] = "fcaCaj in component I_Ca (dimensionless)"
    legend_algebraic[66] = "ibarca_j in component I_Ca (uA_per_uF)"
    legend_algebraic[67] = "ibarca_sl in component I_Ca (uA_per_uF)"
    legend_algebraic[68] = "ibark in component I_Ca (uA_per_uF)"
    legend_algebraic[69] = "ibarna_j in component I_Ca (uA_per_uF)"
    legend_algebraic[70] = "ibarna_sl in component I_Ca (uA_per_uF)"
    legend_algebraic[71] = "I_Ca_junc in component I_Ca (uA_per_uF)"
    legend_algebraic[72] = "I_Ca_sl in component I_Ca (uA_per_uF)"
    legend_algebraic[73] = "I_Ca in component I_Ca (uA_per_uF)"
    legend_algebraic[74] = "I_CaK in component I_Ca (uA_per_uF)"
    legend_algebraic[75] = "I_CaNa_junc in component I_Ca (uA_per_uF)"
    legend_algebraic[76] = "I_CaNa_sl in component I_Ca (uA_per_uF)"
    legend_algebraic[77] = "I_CaNa in component I_Ca (uA_per_uF)"
    legend_algebraic[79] = "I_Catot in component I_Ca (uA_per_uF)"
    legend_algebraic[78] = "Ka_junc in component I_NCX (dimensionless)"
    legend_algebraic[80] = "Ka_sl in component I_NCX (dimensionless)"
    legend_algebraic[81] = "s1_junc in component I_NCX (mM4)"
    legend_algebraic[82] = "s1_sl in component I_NCX (mM4)"
    legend_algebraic[83] = "s2_junc in component I_NCX (mM4)"
    legend_algebraic[84] = "s3_junc in component I_NCX (mM4)"
    legend_algebraic[85] = "s2_sl in component I_NCX (mM4)"
    legend_algebraic[86] = "s3_sl in component I_NCX (mM4)"
    legend_algebraic[87] = "I_ncx_junc in component I_NCX (uA_per_uF)"
    legend_algebraic[88] = "I_ncx_sl in component I_NCX (uA_per_uF)"
    legend_algebraic[90] = "I_ncx in component I_NCX (uA_per_uF)"
    legend_algebraic[91] = "I_pca_junc in component I_PCa (uA_per_uF)"
    legend_algebraic[93] = "I_pca_sl in component I_PCa (uA_per_uF)"
    legend_algebraic[94] = "I_pca in component I_PCa (uA_per_uF)"
    legend_algebraic[95] = "I_cabk_junc in component I_CaBK (uA_per_uF)"
    legend_algebraic[96] = "I_cabk_sl in component I_CaBK (uA_per_uF)"
    legend_algebraic[97] = "I_cabk in component I_CaBK (uA_per_uF)"
    legend_constants[94] = "MaxSR in component SR_Fluxes (dimensionless)"
    legend_constants[95] = "MinSR in component SR_Fluxes (dimensionless)"
    legend_algebraic[10] = "kCaSR in component SR_Fluxes (dimensionless)"
    legend_algebraic[23] = "koSRCa in component SR_Fluxes (per_mM2_per_msec)"
    legend_algebraic[27] = "kiSRCa in component SR_Fluxes (per_mM_per_msec)"
    legend_algebraic[31] = "RI in component SR_Fluxes (mM)"
    legend_states[20] = "Ry_Rr in component SR_Fluxes (mM)"
    legend_states[21] = "Ry_Ro in component SR_Fluxes (mM)"
    legend_states[22] = "Ry_Ri in component SR_Fluxes (mM)"
    legend_algebraic[98] = "J_SRCarel in component SR_Fluxes (mM_per_msec)"
    legend_algebraic[99] = "J_serca in component SR_Fluxes (mM_per_msec)"
    legend_algebraic[100] = "J_SRleak in component SR_Fluxes (mM_per_msec)"
    legend_states[23] = "Ca_sr in component SR_Ca_Concentrations (mM)"
    legend_states[24] = "Ca_i in component Ca_Concentrations (mM)"
    legend_states[25] = "Na_Bj in component Na_Buffers (mM)"
    legend_states[26] = "Na_Bsl in component Na_Buffers (mM)"
    legend_algebraic[102] = "dNa_Bj_dt in component Na_Buffers (mM_per_msec)"
    legend_algebraic[103] = "dNa_Bsl_dt in component Na_Buffers (mM_per_msec)"
    legend_states[27] = "Tn_CL in component Cytosolic_Ca_Buffers (mM)"
    legend_states[28] = "Tn_CHc in component Cytosolic_Ca_Buffers (mM)"
    legend_states[29] = "Tn_CHm in component Cytosolic_Ca_Buffers (mM)"
    legend_states[30] = "CaM in component Cytosolic_Ca_Buffers (mM)"
    legend_states[31] = "Myo_c in component Cytosolic_Ca_Buffers (mM)"
    legend_states[32] = "Myo_m in component Cytosolic_Ca_Buffers (mM)"
    legend_states[33] = "SRB in component Cytosolic_Ca_Buffers (mM)"
    legend_algebraic[101] = "J_CaB_cytosol in component Cytosolic_Ca_Buffers (mM_per_msec)"
    legend_states[34] = "SLL_j in component Junctional_and_SL_Ca_Buffers (mM)"
    legend_states[35] = "SLL_sl in component Junctional_and_SL_Ca_Buffers (mM)"
    legend_states[36] = "SLH_j in component Junctional_and_SL_Ca_Buffers (mM)"
    legend_states[37] = "SLH_sl in component Junctional_and_SL_Ca_Buffers (mM)"
    legend_algebraic[104] = "J_CaB_junction in component Junctional_and_SL_Ca_Buffers (mM_per_msec)"
    legend_algebraic[105] = "J_CaB_sl in component Junctional_and_SL_Ca_Buffers (mM_per_msec)"
    legend_states[38] = "Csqn_b in component SR_Ca_Concentrations (mM)"
    legend_algebraic[106] = "I_Na_tot_junc in component Na_Concentrations (uA_per_uF)"
    legend_algebraic[107] = "I_Na_tot_sl in component Na_Concentrations (uA_per_uF)"
    legend_algebraic[92] = "I_Na_tot_sl2 in component Na_Concentrations (uA_per_uF)"
    legend_algebraic[89] = "I_Na_tot_junc2 in component Na_Concentrations (uA_per_uF)"
    legend_algebraic[108] = "I_K_tot in component K_Concentration (uA_per_uF)"
    legend_algebraic[109] = "I_Ca_tot_junc in component Ca_Concentrations (uA_per_uF)"
    legend_algebraic[110] = "I_Ca_tot_sl in component Ca_Concentrations (uA_per_uF)"
    legend_algebraic[111] = "I_Na_tot in component membrane_potential (uA_per_uF)"
    legend_algebraic[112] = "I_Cl_tot in component membrane_potential (uA_per_uF)"
    legend_algebraic[113] = "I_Ca_tot in component membrane_potential (uA_per_uF)"
    legend_algebraic[114] = "I_tot in component membrane_potential (uA_per_uF)"
    legend_rates[6] = "d/dt m in component I_Na (dimensionless)"
    legend_rates[7] = "d/dt h in component I_Na (dimensionless)"
    legend_rates[8] = "d/dt j in component I_Na (dimensionless)"
    legend_rates[9] = "d/dt x_kr in component I_Kr (dimensionless)"
    legend_rates[10] = "d/dt x_ks in component I_Ks (dimensionless)"
    legend_rates[12] = "d/dt x_to_s in component I_to (dimensionless)"
    legend_rates[13] = "d/dt y_to_s in component I_to (dimensionless)"
    legend_rates[14] = "d/dt x_to_f in component I_to (dimensionless)"
    legend_rates[15] = "d/dt y_to_f in component I_to (dimensionless)"
    legend_rates[16] = "d/dt d in component I_Ca (dimensionless)"
    legend_rates[17] = "d/dt f in component I_Ca (dimensionless)"
    legend_rates[18] = "d/dt f_Ca_Bj in component I_Ca (dimensionless)"
    legend_rates[19] = "d/dt f_Ca_Bsl in component I_Ca (dimensionless)"
    legend_rates[20] = "d/dt Ry_Rr in component SR_Fluxes (mM)"
    legend_rates[21] = "d/dt Ry_Ro in component SR_Fluxes (mM)"
    legend_rates[22] = "d/dt Ry_Ri in component SR_Fluxes (mM)"
    legend_rates[25] = "d/dt Na_Bj in component Na_Buffers (mM)"
    legend_rates[26] = "d/dt Na_Bsl in component Na_Buffers (mM)"
    legend_rates[27] = "d/dt Tn_CL in component Cytosolic_Ca_Buffers (mM)"
    legend_rates[28] = "d/dt Tn_CHc in component Cytosolic_Ca_Buffers (mM)"
    legend_rates[29] = "d/dt Tn_CHm in component Cytosolic_Ca_Buffers (mM)"
    legend_rates[30] = "d/dt CaM in component Cytosolic_Ca_Buffers (mM)"
    legend_rates[31] = "d/dt Myo_c in component Cytosolic_Ca_Buffers (mM)"
    legend_rates[32] = "d/dt Myo_m in component Cytosolic_Ca_Buffers (mM)"
    legend_rates[33] = "d/dt SRB in component Cytosolic_Ca_Buffers (mM)"
    legend_rates[34] = "d/dt SLL_j in component Junctional_and_SL_Ca_Buffers (mM)"
    legend_rates[35] = "d/dt SLL_sl in component Junctional_and_SL_Ca_Buffers (mM)"
    legend_rates[36] = "d/dt SLH_j in component Junctional_and_SL_Ca_Buffers (mM)"
    legend_rates[37] = "d/dt SLH_sl in component Junctional_and_SL_Ca_Buffers (mM)"
    legend_rates[38] = "d/dt Csqn_b in component SR_Ca_Concentrations (mM)"
    legend_rates[23] = "d/dt Ca_sr in component SR_Ca_Concentrations (mM)"
    legend_rates[1] = "d/dt Na_j in component Na_Concentrations (mM)"
    legend_rates[2] = "d/dt Na_sl in component Na_Concentrations (mM)"
    legend_rates[11] = "d/dt Na_i in component Na_Concentrations (mM)"
    legend_rates[3] = "d/dt K_i in component K_Concentration (mM)"
    legend_rates[4] = "d/dt Ca_j in component Ca_Concentrations (mM)"
    legend_rates[5] = "d/dt Ca_sl in component Ca_Concentrations (mM)"
    legend_rates[24] = "d/dt Ca_i in component Ca_Concentrations (mM)"
    legend_rates[0] = "d/dt V_m in component membrane_potential (mV)"
    return (legend_states, legend_algebraic, legend_voi, legend_constants)

def initConsts():
    constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
    states[0] = -8.09763e1
    constants[0] = 0.02
    states[1] = 9.06
    states[2] = 9.06
    states[3] = 120
    states[4] = 1.737475e-4
    states[5] = 1.031812e-4
    states[6] = 1.405627e-3
    states[7] = 9.867005e-1
    states[8] = 9.91562e-1
    states[9] = 8.641386e-3
    states[10] = 5.412034e-3
    states[11] = 9.06
    states[12] = 4.051574e-3
    states[13] = 9.945511e-1
    states[14] = 4.051574e-3
    states[15] = 9.945511e-1
    states[16] = 7.175662e-6
    states[17] = 1.000681
    states[18] = 2.421991e-2
    states[19] = 1.452605e-2
    states[20] = 8.884332e-1
    states[21] = 8.156628e-7
    states[22] = 1.024274e-7
    states[23] = 0.1e-1
    states[24] = 8.597401e-5
    states[25] = 3.539892
    states[26] = 7.720854e-1
    states[27] = 8.773191e-3
    states[28] = 1.078283e-1
    states[29] = 1.524002e-2
    states[30] = 2.911916e-4
    states[31] = 1.298754e-3
    states[32] = 1.381982e-1
    states[33] = 2.143165e-3
    states[34] = 9.566355e-3
    states[35] = 1.110363e-1
    states[36] = 7.347888e-3
    states[37] = 7.297378e-2
    states[38] = 1.242988
    constants[1] = 1.00000
    constants[2] = 8314.00
    constants[3] = 1.38100e-10
    constants[4] = 100.000
    constants[5] = 0.160000
    constants[6] = 0.450000
    constants[7] = 0.500000
    constants[8] = 1.64000e-06
    constants[9] = 1.22000e-06
    constants[10] = 1.09000e-05
    constants[11] = 1.79000e-05
    constants[12] = 8.24130e-13
    constants[13] = 3.27430e-12
    constants[14] = 1.83130e-14
    constants[15] = 1.63860e-12
    constants[16] = 0.110000
    constants[17] = 0.900000
    constants[18] = 5.40000
    constants[19] = 140.000
    constants[20] = 1.80000
    constants[21] = 1.00000
    constants[22] = 23.0000
    constants[23] = 0.000597000
    constants[24] = 1.00000*1.80000
    constants[25] = 11.0000
    constants[26] = 1.50000
    constants[27] = 1.63000
    constants[28] = 1.39000
    constants[29] = 0.0183300
    constants[30] = 2.00000*0.00100000
    constants[31] = 0.500000*0.109625
    constants[32] = 1.00000*0.00900000
    constants[33] = 0.100000
    constants[34] = 0.500000*1.50000e-08
    constants[35] = 0.500000*0.000540000
    constants[36] = 0.500000*2.70000e-07
    constants[37] = 1.80000
    constants[38] = 1.00000*4.50000
    constants[39] = 0.00359000
    constants[40] = 1.30000
    constants[41] = 12.2900
    constants[42] = 87.5000
    constants[43] = 0.320000
    constants[44] = 0.270000
    constants[45] = 0.000150000
    constants[46] = 1.57000
    constants[47] = 0.0673000
    constants[48] = 0.000500000
    constants[49] = 0.000551300
    constants[50] = 2.35000
    constants[51] = 2.60000
    constants[52] = 0.00531140
    constants[53] = 0.000246000
    constants[54] = 1.70000
    constants[55] = 1.78700
    constants[56] = 25.0000
    constants[57] = 10.0000
    constants[58] = 0.0600000
    constants[59] = 0.00500000
    constants[60] = 0.450000
    constants[61] = 7.56100
    constants[62] = 1.65000
    constants[63] = 0.00100000
    constants[64] = 0.000100000
    constants[65] = 0.0700000
    constants[66] = 0.0196000
    constants[67] = 32.7000
    constants[68] = 0.140000
    constants[69] = 3.20000e-05
    constants[70] = 2.37000
    constants[71] = 0.00333000
    constants[72] = 0.00300000
    constants[73] = 0.0240000
    constants[74] = 0.238000
    constants[75] = 34.0000
    constants[76] = 0.140000
    constants[77] = 0.000460000
    constants[78] = 13.8000
    constants[79] = 5.70000e-05
    constants[80] = 0.0157000
    constants[81] = 19.0000*0.000900000
    constants[82] = 0.0600000
    constants[83] = 100.000
    constants[84] = 1.30000
    constants[85] = 100.000
    constants[86] = 0.0300000
    constants[87] = 100.000
    constants[88] = 65.0000
    constants[89] = 100.000
    constants[90] = 0.00350000
    constants[91] = 0.00350000
    constants[92] = 0.00000
    constants[93] = 0.00000
    constants[94] = 15.0000
    constants[95] = 1.00000
    constants[124] = 0.00000
    constants[96] = 96485.0
    constants[97] = 10.2500
    constants[98] = 0.0150000
    constants[99] = 1.00000-constants[16]
    constants[100] = 1.00000-constants[17]
    constants[101] = (exp(constants[19]/67.3000)-1.00000)/7.00000
    constants[102] = 0.0350000*(power(constants[18]/5.40000, 1.0/2))
    constants[103] = custom_piecewise([equal(constants[1] , 1.00000), 1.00000*0.130000*0.120000 , True, 0.130000*0.300000*0.964000])
    constants[104] = custom_piecewise([equal(constants[1] , 1.00000), 1.00000*0.130000*0.880000 , True, 0.130000*0.300000*0.0360000])
    constants[105] = 0.500000
    constants[106] = 310.000
    constants[107] =  pi*(power(constants[97], 2.00000))*constants[4]*1.00000e-15
    constants[108] = 20150.0* pi*2.00000*constants[5]*constants[98]
    constants[109] =  pi*2.00000*constants[97]*constants[4]
    constants[110] = constants[96]/(constants[2]*constants[106])
    constants[111] = (constants[106]-310.000)/10.0000
    constants[112] = 0.650000*constants[107]
    constants[113] = 0.0350000*constants[107]
    constants[114] = 0.0200000*constants[107]
    constants[115] = 0.0539000*0.0100000*constants[107]
    constants[116] = 15.0000
    constants[117] = 150.000
    constants[118] = (0.0374000*constants[112])/constants[114]
    constants[119] = ((0.00460000*constants[112])/constants[115])*0.100000
    constants[120] = (0.0134000*constants[112])/constants[114]
    constants[121] = ((0.00165000*constants[112])/constants[115])*0.100000
    constants[122] = (0.140000*constants[112])/constants[113]
    constants[123] = (1.00000/constants[110])*log(constants[116]/constants[117])
    return (states, constants)

def computeRates(voi, states, constants):
    rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
    rates[3] = constants[124]
    rates[18] = ((1.70000*states[4])/1.00000)*(1.00000-states[18])-0.0119000*states[18]
    rates[19] = ((1.70000*states[5])/1.00000)*(1.00000-states[19])-0.0119000*states[19]
    rates[27] = constants[67]*states[24]*(constants[65]-states[27])-constants[66]*states[27]
    rates[28] = constants[70]*states[24]*((constants[68]-states[28])-states[29])-constants[69]*states[28]
    rates[29] = constants[72]*constants[21]*((constants[68]-states[28])-states[29])-constants[71]*states[29]
    rates[30] = constants[75]*states[24]*(constants[73]-states[30])-constants[74]*states[30]
    rates[31] = constants[78]*states[24]*((constants[76]-states[31])-states[32])-constants[77]*states[31]
    rates[32] = constants[80]*constants[21]*((constants[76]-states[31])-states[32])-constants[79]*states[32]
    rates[33] = constants[83]*states[24]*(constants[81]-states[33])-constants[82]*states[33]
    rates[34] = constants[85]*states[4]*(constants[119]-states[34])-constants[84]*states[34]
    rates[35] = constants[85]*states[5]*(constants[118]-states[35])-constants[84]*states[35]
    rates[36] = constants[87]*states[4]*(constants[121]-states[36])-constants[86]*states[36]
    rates[37] = constants[87]*states[5]*(constants[120]-states[37])-constants[86]*states[37]
    rates[38] = constants[89]*states[23]*(constants[122]-states[38])-constants[88]*states[38]
    rates[11] = (constants[15]/constants[112])*(states[2]-states[11])
    algebraic[0] = 1.00000/(power(1.00000+exp(-(56.8600+states[0])/9.03000), 2.00000))
    algebraic[11] = 0.129200*exp(-(power((states[0]+45.7900)/15.5400, 2.00000)))+0.0648700*exp(-(power((states[0]-4.82300)/51.1200, 2.00000)))
    rates[6] = (algebraic[0]-states[6])/algebraic[11]
    algebraic[3] = 1.00000/(1.00000+exp(-(states[0]+10.0000)/5.00000))
    algebraic[14] = ((550.000/(1.00000+exp((-22.0000-states[0])/9.00000)))*6.00000)/(1.00000+exp((states[0]--11.0000)/9.00000))+230.000/(1.00000+exp((states[0]--40.0000)/20.0000))
    rates[9] = (algebraic[3]-states[9])/algebraic[14]
    algebraic[4] = 1.00000/(1.00000+exp(-(states[0]+3.80000)/14.2500))
    algebraic[15] = 990.100/(1.00000+exp(-(states[0]+2.43600)/14.1200))
    rates[10] = (algebraic[4]-states[10])/algebraic[15]
    algebraic[5] = 1.00000/(1.00000+exp(-(states[0]-19.0000)/13.0000))
    algebraic[16] = 9.00000/(1.00000+exp((states[0]+3.00000)/15.0000))+0.500000
    rates[12] = (algebraic[5]-states[12])/algebraic[16]
    algebraic[6] = 1.00000/(1.00000+exp((states[0]+19.5000)/5.00000))
    algebraic[17] = 800.000/(1.00000+exp((states[0]+60.0000)/10.0000))+30.0000
    rates[13] = (algebraic[6]-states[13])/algebraic[17]
    algebraic[18] = 8.50000*exp(-(power((states[0]+45.0000)/50.0000, 2.00000)))+0.500000
    rates[14] = (algebraic[5]-states[14])/algebraic[18]
    algebraic[19] = 85.0000*exp(-(power(states[0]+40.0000, 2.00000))/220.000)+7.00000
    rates[15] = (algebraic[6]-states[15])/algebraic[19]
    algebraic[8] = 1.00000/(1.00000+exp(-(states[0]+5.00000)/6.00000))
    algebraic[20] = (1.00000*algebraic[8]*(1.00000-exp(-(states[0]+5.00000)/6.00000)))/(0.0350000*(states[0]+5.00000))
    rates[16] = (algebraic[8]-states[16])/algebraic[20]
    algebraic[7] = 1.00000/(1.00000+exp((states[0]+35.0000)/9.00000))+0.600000/(1.00000+exp((50.0000-states[0])/20.0000))
    algebraic[21] = 1.00000/(0.0197000*exp(-(power(0.0337000*(states[0]+14.5000), 2.00000)))+0.0200000)
    rates[17] = (algebraic[7]-states[17])/algebraic[21]
    algebraic[10] = constants[94]-(constants[94]-constants[95])/(1.00000+power(constants[60]/states[23], 2.50000))
    algebraic[23] = constants[57]/algebraic[10]
    algebraic[27] = constants[105]*algebraic[10]
    rates[21] = (algebraic[23]*(power(states[4], 2.00000))*states[20]-constants[58]*states[21])-(algebraic[27]*states[4]*states[21]-constants[59]*states[22])
    algebraic[1] = custom_piecewise([greater_equal(states[0] , -40.0000), 0.00000 , True, 0.0570000*exp(-(states[0]+80.0000)/6.80000)])
    algebraic[12] = custom_piecewise([greater_equal(states[0] , -40.0000), 0.770000/(0.130000*(1.00000+exp(-(states[0]+10.6600)/11.1000))) , True, 2.70000*exp(0.0790000*states[0])+310000.*exp(0.348500*states[0])])
    algebraic[24] = 1.00000/(algebraic[1]+algebraic[12])
    algebraic[28] = 1.00000/(power(1.00000+exp((states[0]+71.5500)/7.43000), 2.00000))
    rates[7] = (algebraic[28]-states[7])/algebraic[24]
    algebraic[2] = custom_piecewise([greater_equal(states[0] , -40.0000), 0.00000 , True, ((-25428.0*exp(0.244400*states[0])-6.94800e-06*exp(-0.0439100*states[0]))*(states[0]+37.7800))/(1.00000+exp(0.311000*(states[0]+79.2300)))])
    algebraic[13] = custom_piecewise([greater_equal(states[0] , -40.0000), (0.600000*exp(0.0570000*states[0]))/(1.00000+exp(-0.100000*(states[0]+32.0000))) , True, (0.0242400*exp(-0.0105200*states[0]))/(1.00000+exp(-0.137800*(states[0]+40.1400)))])
    algebraic[25] = 1.00000/(algebraic[2]+algebraic[13])
    algebraic[29] = 1.00000/(power(1.00000+exp((states[0]+71.5500)/7.43000), 2.00000))
    rates[8] = (algebraic[29]-states[8])/algebraic[25]
    algebraic[31] = ((1.00000-states[20])-states[21])-states[22]
    rates[20] = (constants[59]*algebraic[31]-algebraic[27]*states[4]*states[20])-(algebraic[23]*(power(states[4], 2.00000))*states[20]-constants[58]*states[21])
    rates[22] = (algebraic[27]*states[4]*states[21]-constants[59]*states[22])-(constants[58]*states[22]-algebraic[23]*(power(states[4], 2.00000))*algebraic[31])
    algebraic[98] = ((constants[56]*states[21])/1.00000)*(states[23]-states[4])
    algebraic[99] = ((power(constants[51], constants[111]))*constants[52]*(power(states[24]/constants[53], constants[55])-power(states[23]/constants[54], constants[55])))/(1.00000+power(states[24]/constants[53], constants[55])+power(states[23]/constants[54], constants[55]))
    algebraic[100] = 5.34800e-06*(states[23]-states[4])
    rates[23] = (algebraic[99]-((algebraic[100]*constants[112])/constants[113]+algebraic[98]))-(constants[89]*states[23]*(constants[122]-states[38])-constants[88]*states[38])
    algebraic[101] = (((((((((((constants[67]*states[24]*(constants[65]-states[27])-constants[66]*states[27])+constants[70]*states[24]*((constants[68]-states[28])-states[29]))-constants[69]*states[28])+constants[72]*constants[21]*((constants[68]-states[28])-states[29]))-constants[71]*states[29])+constants[75]*states[24]*(constants[73]-states[30]))-constants[74]*states[30])+constants[78]*states[24]*((constants[76]-states[31])-states[32]))-constants[77]*states[31])+constants[80]*constants[21]*((constants[76]-states[31])-states[32]))-constants[79]*states[32])+(constants[83]*states[24]*(constants[81]-states[33])-constants[82]*states[33])
    rates[24] = ((-algebraic[99]*constants[113])/constants[112]-algebraic[101])+(constants[13]/constants[112])*(states[5]-states[24])
    algebraic[102] = constants[64]*states[1]*(constants[61]-states[25])-constants[63]*states[25]
    rates[25] = algebraic[102]
    algebraic[103] = constants[64]*states[2]*(constants[62]-states[26])-constants[63]*states[26]
    rates[26] = algebraic[103]
    algebraic[26] = (1.00000/constants[110])*log(constants[19]/states[1])
    algebraic[35] = constants[16]*constants[22]*(power(states[6], 3.00000))*states[7]*states[8]*(states[0]-algebraic[26])
    algebraic[38] = constants[16]*constants[23]*(states[0]-algebraic[26])
    algebraic[41] = 1.00000/(1.00000+0.124500*exp(-0.100000*states[0]*constants[110])+0.0365000*constants[101]*exp(-states[0]*constants[110]))
    algebraic[42] = ((constants[16]*constants[24]*algebraic[41]*constants[18])/(1.00000+power(constants[25]/states[1], 4.00000)))/(constants[18]+constants[26])
    algebraic[69] = (constants[34]*states[0]*constants[96]*constants[110]*(0.750000*states[1]*exp(states[0]*constants[110])-0.750000*constants[19]))/(exp(states[0]*constants[110])-1.00000)
    algebraic[75] = constants[17]*algebraic[69]*states[16]*states[17]*((1.00000-states[18])+constants[93])*(power(constants[37], constants[111]))*0.450000*1.00000
    algebraic[78] = 1.00000/(1.00000+power(constants[45]/states[4], 2.00000))
    algebraic[81] = exp(constants[44]*states[0]*constants[110])*(power(states[1], 3.00000))*constants[20]
    algebraic[83] = exp((constants[44]-1.00000)*states[0]*constants[110])*(power(constants[19], 3.00000))*states[4]
    algebraic[84] = constants[39]*(power(constants[19], 3.00000))*(1.00000+power(states[1]/constants[41], 3.00000))+(power(constants[42], 3.00000))*states[4]*(1.00000+states[4]/constants[39])+constants[40]*(power(states[1], 3.00000))+(power(states[1], 3.00000))*constants[20]+(power(constants[19], 3.00000))*states[4]
    algebraic[87] = ((constants[16]*constants[38]*(power(constants[46], constants[111]))*algebraic[78]*(algebraic[81]-algebraic[83]))/algebraic[84])/(1.00000+constants[43]*exp((constants[44]-1.00000)*states[0]*constants[110]))
    algebraic[106] = algebraic[35]+algebraic[38]+3.00000*algebraic[87]+3.00000*algebraic[42]+algebraic[75]
    rates[1] = ((-algebraic[106]*constants[3])/(constants[115]*constants[96])+(constants[14]/constants[115])*(states[2]-states[1]))-algebraic[102]
    algebraic[30] = (1.00000/constants[110])*log(constants[19]/states[2])
    algebraic[36] = constants[99]*constants[22]*(power(states[6], 3.00000))*states[7]*states[8]*(states[0]-algebraic[30])
    algebraic[39] = constants[99]*constants[23]*(states[0]-algebraic[30])
    algebraic[43] = ((constants[99]*constants[24]*algebraic[41]*constants[18])/(1.00000+power(constants[25]/states[2], 4.00000)))/(constants[18]+constants[26])
    algebraic[70] = (constants[34]*states[0]*constants[96]*constants[110]*(0.750000*states[2]*exp(states[0]*constants[110])-0.750000*constants[19]))/(exp(states[0]*constants[110])-1.00000)
    algebraic[76] = constants[100]*algebraic[70]*states[16]*states[17]*((1.00000-states[19])+constants[92])*(power(constants[37], constants[111]))*0.450000*1.00000
    algebraic[80] = 1.00000/(1.00000+power(constants[45]/states[5], 2.00000))
    algebraic[82] = exp(constants[44]*states[0]*constants[110])*(power(states[2], 3.00000))*constants[20]
    algebraic[85] = exp((constants[44]-1.00000)*states[0]*constants[110])*(power(constants[19], 3.00000))*states[5]
    algebraic[86] = constants[39]*(power(constants[19], 3.00000))*(1.00000+power(states[2]/constants[41], 3.00000))+(power(constants[42], 3.00000))*states[5]*(1.00000+states[5]/constants[39])+constants[40]*(power(states[2], 3.00000))+(power(states[2], 3.00000))*constants[20]+(power(constants[19], 3.00000))*states[5]
    algebraic[88] = ((constants[99]*constants[38]*(power(constants[46], constants[111]))*algebraic[80]*(algebraic[82]-algebraic[85]))/algebraic[86])/(1.00000+constants[43]*exp((constants[44]-1.00000)*states[0]*constants[110]))
    algebraic[107] = algebraic[36]+algebraic[39]+3.00000*algebraic[88]+3.00000*algebraic[43]+algebraic[76]
    rates[2] = ((-algebraic[107]*constants[3])/(constants[114]*constants[96])+(constants[14]/constants[114])*(states[1]-states[2])+(constants[15]/constants[114])*(states[11]-states[2]))-algebraic[103]
    algebraic[104] = (constants[85]*states[4]*(constants[119]-states[34])-constants[84]*states[34])+(constants[87]*states[4]*(constants[121]-states[36])-constants[86]*states[36])
    algebraic[66] = (constants[35]*4.00000*states[0]*constants[96]*constants[110]*(0.341000*states[4]*exp(2.00000*states[0]*constants[110])-0.341000*constants[20]))/(exp(2.00000*states[0]*constants[110])-1.00000)
    algebraic[71] = constants[17]*algebraic[66]*states[16]*states[17]*((1.00000-states[18])+constants[93])*(power(constants[37], constants[111]))*0.450000*1.00000
    algebraic[91] = (constants[16]*(power(constants[50], constants[111]))*constants[47]*(power(states[4], 1.60000)))/(power(constants[48], 1.60000)+power(states[4], 1.60000))
    algebraic[33] = ((1.00000/constants[110])/2.00000)*log(constants[20]/states[4])
    algebraic[95] = constants[16]*constants[49]*(states[0]-algebraic[33])
    algebraic[109] = (algebraic[71]+algebraic[95]+algebraic[91])-2.00000*algebraic[87]
    rates[4] = (((-algebraic[109]*constants[3])/(constants[115]*2.00000*constants[96])+(constants[12]/constants[115])*(states[5]-states[4]))-algebraic[104])+(algebraic[98]*constants[113])/constants[115]+(algebraic[100]*constants[112])/constants[115]
    algebraic[105] = (constants[85]*states[5]*(constants[118]-states[35])-constants[84]*states[35])+(constants[87]*states[5]*(constants[120]-states[37])-constants[86]*states[37])
    algebraic[67] = (constants[35]*4.00000*states[0]*constants[96]*constants[110]*(0.341000*states[5]*exp(2.00000*states[0]*constants[110])-0.341000*constants[20]))/(exp(2.00000*states[0]*constants[110])-1.00000)
    algebraic[72] = constants[100]*algebraic[67]*states[16]*states[17]*((1.00000-states[19])+constants[92])*(power(constants[37], constants[111]))*0.450000*1.00000
    algebraic[93] = (constants[99]*(power(constants[50], constants[111]))*constants[47]*(power(states[5], 1.60000)))/(power(constants[48], 1.60000)+power(states[5], 1.60000))
    algebraic[34] = ((1.00000/constants[110])/2.00000)*log(constants[20]/states[5])
    algebraic[96] = constants[99]*constants[49]*(states[0]-algebraic[34])
    algebraic[110] = (algebraic[72]+algebraic[96]+algebraic[93])-2.00000*algebraic[88]
    rates[5] = ((-algebraic[110]*constants[3])/(constants[114]*2.00000*constants[96])+(constants[12]/constants[114])*(states[4]-states[5])+(constants[13]/constants[114])*(states[24]-states[5]))-algebraic[105]
    algebraic[9] = custom_piecewise([less_equal( voi+1800.00 % 2000.00 , 0.00000), 30.0000 , True, -40.0000])
    algebraic[22] = (algebraic[9]-states[0])/constants[0]
    algebraic[44] = algebraic[42]+algebraic[43]
    algebraic[32] = (1.00000/constants[110])*log(constants[18]/states[3])
    algebraic[45] = 1.00000/(1.00000+exp((states[0]+74.0000)/24.0000))
    algebraic[46] = constants[102]*states[9]*algebraic[45]*(states[0]-algebraic[32])
    algebraic[47] = 1.00000/(1.00000+exp(7.48800-states[0]/5.98000))
    algebraic[48] = constants[16]*constants[30]*algebraic[47]*(states[0]-algebraic[32])
    algebraic[49] = constants[99]*constants[30]*algebraic[47]*(states[0]-algebraic[32])
    algebraic[50] = algebraic[48]+algebraic[49]
    algebraic[51] = (1.00000/constants[110])*log((constants[18]+constants[29]*constants[19])/(states[3]+constants[29]*states[11]))
    algebraic[52] = constants[16]*constants[90]*(power(states[10], 2.00000))*(states[0]-algebraic[51])
    algebraic[53] = constants[99]*constants[91]*(power(states[10], 2.00000))*(states[0]-algebraic[51])
    algebraic[54] = algebraic[52]+algebraic[53]
    algebraic[55] = constants[103]*states[12]*states[13]*(states[0]-algebraic[32])
    algebraic[56] = constants[104]*states[14]*states[15]*(states[0]-algebraic[32])
    algebraic[57] = algebraic[55]+algebraic[56]
    algebraic[58] = 1.02000/(1.00000+exp(0.238500*((states[0]-algebraic[32])-59.2150)))
    algebraic[59] = (0.491240*exp(0.0803200*((states[0]+5.47600)-algebraic[32]))+exp(0.0617500*((states[0]-algebraic[32])-594.310)))/(1.00000+exp(-0.514300*((states[0]-algebraic[32])+4.75300)))
    algebraic[60] = algebraic[58]/(algebraic[58]+algebraic[59])
    algebraic[61] = 1.00000*0.350000*(power(constants[18]/5.40000, 1.0/2))*algebraic[60]*(states[0]-algebraic[32])
    algebraic[68] = (constants[36]*states[0]*constants[96]*constants[110]*(0.750000*states[3]*exp(states[0]*constants[110])-0.750000*constants[18]))/(exp(states[0]*constants[110])-1.00000)
    algebraic[74] = algebraic[68]*states[16]*states[17]*(constants[17]*(constants[93]+(1.00000-states[18]))+constants[100]*(constants[92]+(1.00000-states[19])))*(power(constants[37], constants[111]))*0.450000*1.00000
    algebraic[108] = ((algebraic[57]+algebraic[46]+algebraic[54]+algebraic[61])-2.00000*algebraic[44])+algebraic[74]+algebraic[50]
    algebraic[111] = algebraic[106]+algebraic[107]
    algebraic[62] = ((constants[16]*constants[31])/(1.00000+constants[33]/states[4]))*(states[0]-constants[123])
    algebraic[63] = ((constants[99]*constants[31])/(1.00000+constants[33]/states[5]))*(states[0]-constants[123])
    algebraic[64] = algebraic[62]+algebraic[63]
    algebraic[65] = constants[32]*(states[0]-constants[123])
    algebraic[112] = algebraic[64]+algebraic[65]
    algebraic[113] = algebraic[109]+algebraic[110]
    algebraic[114] = algebraic[111]+algebraic[112]+algebraic[113]+algebraic[108]
    rates[0] = -(algebraic[114]-algebraic[22])
    return(rates)

def computeAlgebraic(constants, states, voi):
    algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
    states = array(states)
    voi = array(voi)
    algebraic[0] = 1.00000/(power(1.00000+exp(-(56.8600+states[0])/9.03000), 2.00000))
    algebraic[11] = 0.129200*exp(-(power((states[0]+45.7900)/15.5400, 2.00000)))+0.0648700*exp(-(power((states[0]-4.82300)/51.1200, 2.00000)))
    algebraic[3] = 1.00000/(1.00000+exp(-(states[0]+10.0000)/5.00000))
    algebraic[14] = ((550.000/(1.00000+exp((-22.0000-states[0])/9.00000)))*6.00000)/(1.00000+exp((states[0]--11.0000)/9.00000))+230.000/(1.00000+exp((states[0]--40.0000)/20.0000))
    algebraic[4] = 1.00000/(1.00000+exp(-(states[0]+3.80000)/14.2500))
    algebraic[15] = 990.100/(1.00000+exp(-(states[0]+2.43600)/14.1200))
    algebraic[5] = 1.00000/(1.00000+exp(-(states[0]-19.0000)/13.0000))
    algebraic[16] = 9.00000/(1.00000+exp((states[0]+3.00000)/15.0000))+0.500000
    algebraic[6] = 1.00000/(1.00000+exp((states[0]+19.5000)/5.00000))
    algebraic[17] = 800.000/(1.00000+exp((states[0]+60.0000)/10.0000))+30.0000
    algebraic[18] = 8.50000*exp(-(power((states[0]+45.0000)/50.0000, 2.00000)))+0.500000
    algebraic[19] = 85.0000*exp(-(power(states[0]+40.0000, 2.00000))/220.000)+7.00000
    algebraic[8] = 1.00000/(1.00000+exp(-(states[0]+5.00000)/6.00000))
    algebraic[20] = (1.00000*algebraic[8]*(1.00000-exp(-(states[0]+5.00000)/6.00000)))/(0.0350000*(states[0]+5.00000))
    algebraic[7] = 1.00000/(1.00000+exp((states[0]+35.0000)/9.00000))+0.600000/(1.00000+exp((50.0000-states[0])/20.0000))
    algebraic[21] = 1.00000/(0.0197000*exp(-(power(0.0337000*(states[0]+14.5000), 2.00000)))+0.0200000)
    algebraic[10] = constants[94]-(constants[94]-constants[95])/(1.00000+power(constants[60]/states[23], 2.50000))
    algebraic[23] = constants[57]/algebraic[10]
    algebraic[27] = constants[105]*algebraic[10]
    algebraic[1] = custom_piecewise([greater_equal(states[0] , -40.0000), 0.00000 , True, 0.0570000*exp(-(states[0]+80.0000)/6.80000)])
    algebraic[12] = custom_piecewise([greater_equal(states[0] , -40.0000), 0.770000/(0.130000*(1.00000+exp(-(states[0]+10.6600)/11.1000))) , True, 2.70000*exp(0.0790000*states[0])+310000.*exp(0.348500*states[0])])
    algebraic[24] = 1.00000/(algebraic[1]+algebraic[12])
    algebraic[28] = 1.00000/(power(1.00000+exp((states[0]+71.5500)/7.43000), 2.00000))
    algebraic[2] = custom_piecewise([greater_equal(states[0] , -40.0000), 0.00000 , True, ((-25428.0*exp(0.244400*states[0])-6.94800e-06*exp(-0.0439100*states[0]))*(states[0]+37.7800))/(1.00000+exp(0.311000*(states[0]+79.2300)))])
    algebraic[13] = custom_piecewise([greater_equal(states[0] , -40.0000), (0.600000*exp(0.0570000*states[0]))/(1.00000+exp(-0.100000*(states[0]+32.0000))) , True, (0.0242400*exp(-0.0105200*states[0]))/(1.00000+exp(-0.137800*(states[0]+40.1400)))])
    algebraic[25] = 1.00000/(algebraic[2]+algebraic[13])
    algebraic[29] = 1.00000/(power(1.00000+exp((states[0]+71.5500)/7.43000), 2.00000))
    algebraic[31] = ((1.00000-states[20])-states[21])-states[22]
    algebraic[98] = ((constants[56]*states[21])/1.00000)*(states[23]-states[4])
    algebraic[99] = ((power(constants[51], constants[111]))*constants[52]*(power(states[24]/constants[53], constants[55])-power(states[23]/constants[54], constants[55])))/(1.00000+power(states[24]/constants[53], constants[55])+power(states[23]/constants[54], constants[55]))
    algebraic[100] = 5.34800e-06*(states[23]-states[4])
    algebraic[101] = (((((((((((constants[67]*states[24]*(constants[65]-states[27])-constants[66]*states[27])+constants[70]*states[24]*((constants[68]-states[28])-states[29]))-constants[69]*states[28])+constants[72]*constants[21]*((constants[68]-states[28])-states[29]))-constants[71]*states[29])+constants[75]*states[24]*(constants[73]-states[30]))-constants[74]*states[30])+constants[78]*states[24]*((constants[76]-states[31])-states[32]))-constants[77]*states[31])+constants[80]*constants[21]*((constants[76]-states[31])-states[32]))-constants[79]*states[32])+(constants[83]*states[24]*(constants[81]-states[33])-constants[82]*states[33])
    algebraic[102] = constants[64]*states[1]*(constants[61]-states[25])-constants[63]*states[25]
    algebraic[103] = constants[64]*states[2]*(constants[62]-states[26])-constants[63]*states[26]
    algebraic[26] = (1.00000/constants[110])*log(constants[19]/states[1])
    algebraic[35] = constants[16]*constants[22]*(power(states[6], 3.00000))*states[7]*states[8]*(states[0]-algebraic[26])
    algebraic[38] = constants[16]*constants[23]*(states[0]-algebraic[26])
    algebraic[41] = 1.00000/(1.00000+0.124500*exp(-0.100000*states[0]*constants[110])+0.0365000*constants[101]*exp(-states[0]*constants[110]))
    algebraic[42] = ((constants[16]*constants[24]*algebraic[41]*constants[18])/(1.00000+power(constants[25]/states[1], 4.00000)))/(constants[18]+constants[26])
    algebraic[69] = (constants[34]*states[0]*constants[96]*constants[110]*(0.750000*states[1]*exp(states[0]*constants[110])-0.750000*constants[19]))/(exp(states[0]*constants[110])-1.00000)
    algebraic[75] = constants[17]*algebraic[69]*states[16]*states[17]*((1.00000-states[18])+constants[93])*(power(constants[37], constants[111]))*0.450000*1.00000
    algebraic[78] = 1.00000/(1.00000+power(constants[45]/states[4], 2.00000))
    algebraic[81] = exp(constants[44]*states[0]*constants[110])*(power(states[1], 3.00000))*constants[20]
    algebraic[83] = exp((constants[44]-1.00000)*states[0]*constants[110])*(power(constants[19], 3.00000))*states[4]
    algebraic[84] = constants[39]*(power(constants[19], 3.00000))*(1.00000+power(states[1]/constants[41], 3.00000))+(power(constants[42], 3.00000))*states[4]*(1.00000+states[4]/constants[39])+constants[40]*(power(states[1], 3.00000))+(power(states[1], 3.00000))*constants[20]+(power(constants[19], 3.00000))*states[4]
    algebraic[87] = ((constants[16]*constants[38]*(power(constants[46], constants[111]))*algebraic[78]*(algebraic[81]-algebraic[83]))/algebraic[84])/(1.00000+constants[43]*exp((constants[44]-1.00000)*states[0]*constants[110]))
    algebraic[106] = algebraic[35]+algebraic[38]+3.00000*algebraic[87]+3.00000*algebraic[42]+algebraic[75]
    algebraic[30] = (1.00000/constants[110])*log(constants[19]/states[2])
    algebraic[36] = constants[99]*constants[22]*(power(states[6], 3.00000))*states[7]*states[8]*(states[0]-algebraic[30])
    algebraic[39] = constants[99]*constants[23]*(states[0]-algebraic[30])
    algebraic[43] = ((constants[99]*constants[24]*algebraic[41]*constants[18])/(1.00000+power(constants[25]/states[2], 4.00000)))/(constants[18]+constants[26])
    algebraic[70] = (constants[34]*states[0]*constants[96]*constants[110]*(0.750000*states[2]*exp(states[0]*constants[110])-0.750000*constants[19]))/(exp(states[0]*constants[110])-1.00000)
    algebraic[76] = constants[100]*algebraic[70]*states[16]*states[17]*((1.00000-states[19])+constants[92])*(power(constants[37], constants[111]))*0.450000*1.00000
    algebraic[80] = 1.00000/(1.00000+power(constants[45]/states[5], 2.00000))
    algebraic[82] = exp(constants[44]*states[0]*constants[110])*(power(states[2], 3.00000))*constants[20]
    algebraic[85] = exp((constants[44]-1.00000)*states[0]*constants[110])*(power(constants[19], 3.00000))*states[5]
    algebraic[86] = constants[39]*(power(constants[19], 3.00000))*(1.00000+power(states[2]/constants[41], 3.00000))+(power(constants[42], 3.00000))*states[5]*(1.00000+states[5]/constants[39])+constants[40]*(power(states[2], 3.00000))+(power(states[2], 3.00000))*constants[20]+(power(constants[19], 3.00000))*states[5]
    algebraic[88] = ((constants[99]*constants[38]*(power(constants[46], constants[111]))*algebraic[80]*(algebraic[82]-algebraic[85]))/algebraic[86])/(1.00000+constants[43]*exp((constants[44]-1.00000)*states[0]*constants[110]))
    algebraic[107] = algebraic[36]+algebraic[39]+3.00000*algebraic[88]+3.00000*algebraic[43]+algebraic[76]
    algebraic[104] = (constants[85]*states[4]*(constants[119]-states[34])-constants[84]*states[34])+(constants[87]*states[4]*(constants[121]-states[36])-constants[86]*states[36])
    algebraic[66] = (constants[35]*4.00000*states[0]*constants[96]*constants[110]*(0.341000*states[4]*exp(2.00000*states[0]*constants[110])-0.341000*constants[20]))/(exp(2.00000*states[0]*constants[110])-1.00000)
    algebraic[71] = constants[17]*algebraic[66]*states[16]*states[17]*((1.00000-states[18])+constants[93])*(power(constants[37], constants[111]))*0.450000*1.00000
    algebraic[91] = (constants[16]*(power(constants[50], constants[111]))*constants[47]*(power(states[4], 1.60000)))/(power(constants[48], 1.60000)+power(states[4], 1.60000))
    algebraic[33] = ((1.00000/constants[110])/2.00000)*log(constants[20]/states[4])
    algebraic[95] = constants[16]*constants[49]*(states[0]-algebraic[33])
    algebraic[109] = (algebraic[71]+algebraic[95]+algebraic[91])-2.00000*algebraic[87]
    algebraic[105] = (constants[85]*states[5]*(constants[118]-states[35])-constants[84]*states[35])+(constants[87]*states[5]*(constants[120]-states[37])-constants[86]*states[37])
    algebraic[67] = (constants[35]*4.00000*states[0]*constants[96]*constants[110]*(0.341000*states[5]*exp(2.00000*states[0]*constants[110])-0.341000*constants[20]))/(exp(2.00000*states[0]*constants[110])-1.00000)
    algebraic[72] = constants[100]*algebraic[67]*states[16]*states[17]*((1.00000-states[19])+constants[92])*(power(constants[37], constants[111]))*0.450000*1.00000
    algebraic[93] = (constants[99]*(power(constants[50], constants[111]))*constants[47]*(power(states[5], 1.60000)))/(power(constants[48], 1.60000)+power(states[5], 1.60000))
    algebraic[34] = ((1.00000/constants[110])/2.00000)*log(constants[20]/states[5])
    algebraic[96] = constants[99]*constants[49]*(states[0]-algebraic[34])
    algebraic[110] = (algebraic[72]+algebraic[96]+algebraic[93])-2.00000*algebraic[88]
    algebraic[9] = custom_piecewise([less_equal( voi+1800.00 % 2000.00 , 0.00000), 30.0000 , True, -40.0000])
    algebraic[22] = (algebraic[9]-states[0])/constants[0]
    algebraic[44] = algebraic[42]+algebraic[43]
    algebraic[32] = (1.00000/constants[110])*log(constants[18]/states[3])
    algebraic[45] = 1.00000/(1.00000+exp((states[0]+74.0000)/24.0000))
    algebraic[46] = constants[102]*states[9]*algebraic[45]*(states[0]-algebraic[32])
    algebraic[47] = 1.00000/(1.00000+exp(7.48800-states[0]/5.98000))
    algebraic[48] = constants[16]*constants[30]*algebraic[47]*(states[0]-algebraic[32])
    algebraic[49] = constants[99]*constants[30]*algebraic[47]*(states[0]-algebraic[32])
    algebraic[50] = algebraic[48]+algebraic[49]
    algebraic[51] = (1.00000/constants[110])*log((constants[18]+constants[29]*constants[19])/(states[3]+constants[29]*states[11]))
    algebraic[52] = constants[16]*constants[90]*(power(states[10], 2.00000))*(states[0]-algebraic[51])
    algebraic[53] = constants[99]*constants[91]*(power(states[10], 2.00000))*(states[0]-algebraic[51])
    algebraic[54] = algebraic[52]+algebraic[53]
    algebraic[55] = constants[103]*states[12]*states[13]*(states[0]-algebraic[32])
    algebraic[56] = constants[104]*states[14]*states[15]*(states[0]-algebraic[32])
    algebraic[57] = algebraic[55]+algebraic[56]
    algebraic[58] = 1.02000/(1.00000+exp(0.238500*((states[0]-algebraic[32])-59.2150)))
    algebraic[59] = (0.491240*exp(0.0803200*((states[0]+5.47600)-algebraic[32]))+exp(0.0617500*((states[0]-algebraic[32])-594.310)))/(1.00000+exp(-0.514300*((states[0]-algebraic[32])+4.75300)))
    algebraic[60] = algebraic[58]/(algebraic[58]+algebraic[59])
    algebraic[61] = 1.00000*0.350000*(power(constants[18]/5.40000, 1.0/2))*algebraic[60]*(states[0]-algebraic[32])
    algebraic[68] = (constants[36]*states[0]*constants[96]*constants[110]*(0.750000*states[3]*exp(states[0]*constants[110])-0.750000*constants[18]))/(exp(states[0]*constants[110])-1.00000)
    algebraic[74] = algebraic[68]*states[16]*states[17]*(constants[17]*(constants[93]+(1.00000-states[18]))+constants[100]*(constants[92]+(1.00000-states[19])))*(power(constants[37], constants[111]))*0.450000*1.00000
    algebraic[108] = ((algebraic[57]+algebraic[46]+algebraic[54]+algebraic[61])-2.00000*algebraic[44])+algebraic[74]+algebraic[50]
    algebraic[111] = algebraic[106]+algebraic[107]
    algebraic[62] = ((constants[16]*constants[31])/(1.00000+constants[33]/states[4]))*(states[0]-constants[123])
    algebraic[63] = ((constants[99]*constants[31])/(1.00000+constants[33]/states[5]))*(states[0]-constants[123])
    algebraic[64] = algebraic[62]+algebraic[63]
    algebraic[65] = constants[32]*(states[0]-constants[123])
    algebraic[112] = algebraic[64]+algebraic[65]
    algebraic[113] = algebraic[109]+algebraic[110]
    algebraic[114] = algebraic[111]+algebraic[112]+algebraic[113]+algebraic[108]
    algebraic[37] = algebraic[35]+algebraic[36]
    algebraic[40] = algebraic[38]+algebraic[39]
    algebraic[73] = algebraic[71]+algebraic[72]
    algebraic[77] = algebraic[75]+algebraic[76]
    algebraic[79] = algebraic[73]+algebraic[74]+algebraic[77]
    algebraic[89] = 3.00000*algebraic[87]+3.00000*algebraic[42]+algebraic[75]
    algebraic[90] = algebraic[87]+algebraic[88]
    algebraic[92] = 3.00000*algebraic[88]+3.00000*algebraic[43]+algebraic[76]
    algebraic[94] = algebraic[91]+algebraic[93]
    algebraic[97] = algebraic[95]+algebraic[96]
    return algebraic

def custom_piecewise(cases):
    """Compute result of a piecewise function"""
    return select(cases[0::2],cases[1::2])

def solve_model():
    """Solve model with ODE solver"""
    from scipy.integrate import ode
    # Initialise constants and state variables
    (init_states, constants) = initConsts()

    # Set timespan to solve over
    voi = linspace(0, 10, 500)

    # Construct ODE object to solve
    r = ode(computeRates)
    r.set_integrator('vode', method='bdf', atol=1e-06, rtol=1e-06, max_step=1)
    r.set_initial_value(init_states, voi[0])
    r.set_f_params(constants)

    # Solve model
    states = array([[0.0] * len(voi)] * sizeStates)
    states[:,0] = init_states
    for (i,t) in enumerate(voi[1:]):
        if r.successful():
            r.integrate(t)
            states[:,i+1] = r.y
        else:
            break

    # Compute algebraic variables
    algebraic = computeAlgebraic(constants, states, voi)
    return (voi, states, algebraic)

def plot_model(voi, states, algebraic):
    """Plot variables against variable of integration"""
    import pylab
    (legend_states, legend_algebraic, legend_voi, legend_constants) = createLegends()
    pylab.figure(1)
    pylab.plot(voi,vstack((states,algebraic)).T)
    pylab.xlabel(legend_voi)
    pylab.legend(legend_states + legend_algebraic, loc='best')
    pylab.show()

if __name__ == "__main__":
    (voi, states, algebraic) = solve_model()
    plot_model(voi, states, algebraic)