# Size of variable arrays:
sizeAlgebraic = 89
sizeStates = 22
sizeConstants = 151
from math import *
from numpy import *

def createLegends():
    legend_states = [""] * sizeStates
    legend_rates = [""] * sizeStates
    legend_algebraic = [""] * sizeAlgebraic
    legend_voi = ""
    legend_constants = [""] * sizeConstants
    legend_constants[0] = "Buf_C in component cai (mM)"
    legend_states[0] = "Cai in component cai (mM)"
    legend_algebraic[0] = "Cai_bufc in component cai (dimensionless)"
    legend_constants[1] = "Cm in component geom (pF)"
    legend_constants[2] = "F in component phys (C_per_mmol)"
    legend_constants[3] = "Kbuf_C in component cai (mM)"
    legend_constants[144] = "Vc in component geom (um3)"
    legend_algebraic[65] = "i_CaL_Ca in component ical (A_per_F)"
    legend_algebraic[71] = "i_CaT in component icat (A_per_F)"
    legend_algebraic[83] = "i_NaCa in component inaca (A_per_F)"
    legend_algebraic[85] = "i_PCa in component ipca (A_per_F)"
    legend_algebraic[62] = "i_b_Ca in component ibca (A_per_F)"
    legend_algebraic[81] = "i_leak in component ileak (mM_per_ms)"
    legend_algebraic[86] = "i_rel in component irel (mM_per_ms)"
    legend_algebraic[88] = "i_up in component iup (mM_per_ms)"
    legend_voi = "time in component engine (ms)"
    legend_constants[86] = "Buf_SR in component casr (mM)"
    legend_states[1] = "Ca_SR in component casr (mM)"
    legend_algebraic[16] = "Ca_SR_bufSR in component casr (dimensionless)"
    legend_constants[4] = "Kbuf_SR in component casr (mM)"
    legend_constants[126] = "V_SR in component geom (um3)"
    legend_constants[5] = "pace in component engine (dimensionless)"
    legend_constants[6] = "Cao in component extra (mM)"
    legend_algebraic[32] = "E_Ca in component erev (mV)"
    legend_algebraic[46] = "E_K in component erev (mV)"
    legend_algebraic[60] = "E_Na in component erev (mV)"
    legend_states[2] = "Ki in component ki (mM)"
    legend_constants[7] = "Ko in component extra (mM)"
    legend_states[3] = "Nai in component nai (mM)"
    legend_constants[8] = "Nao in component extra (mM)"
    legend_constants[87] = "RTF in component phys (mV)"
    legend_constants[9] = "VSR_tenT in component geom (um3)"
    legend_constants[10] = "V_tot in component geom (um3)"
    legend_constants[85] = "V_tot_tenT in component geom (um3)"
    legend_constants[11] = "Vc_tenT in component geom (um3)"
    legend_states[4] = "V in component membrane (mV)"
    legend_constants[88] = "g_b_Ca in component ibca (mS_per_uF)"
    legend_constants[89] = "g_b_Na in component ibna (mS_per_uF)"
    legend_algebraic[63] = "i_b_Na in component ibna (A_per_F)"
    legend_constants[127] = "FFRT in component phys (s4_A2_per_g_per_m2_per_mol_times_1e3)"
    legend_constants[90] = "FRT in component phys (per_mV)"
    legend_algebraic[3] = "alpha_fCa in component ical (dimensionless)"
    legend_algebraic[19] = "beta_fCa in component ical (dimensionless)"
    legend_states[5] = "d in component ical (dimensionless)"
    legend_constants[12] = "d1 in component ical (mS_per_uF)"
    legend_constants[13] = "d2 in component ical (mV)"
    legend_constants[91] = "d3 in component ical (mS_per_uF)"
    legend_constants[128] = "d4 in component ical (mV)"
    legend_constants[14] = "d5 in component ical (dimensionless)"
    legend_constants[92] = "d6 in component ical (mV)"
    legend_states[6] = "f in component ical (dimensionless)"
    legend_constants[15] = "f1 in component ical (mS_per_uF)"
    legend_constants[93] = "f2 in component ical (mV)"
    legend_constants[94] = "f3 in component ical (mS_per_uF)"
    legend_constants[129] = "f4 in component ical (mV)"
    legend_constants[16] = "f5 in component ical (dimensionless)"
    legend_constants[17] = "f6 in component ical (mV)"
    legend_states[7] = "fCa in component ical (dimensionless)"
    legend_algebraic[49] = "fCa_inf in component ical (dimensionless)"
    legend_algebraic[35] = "gamma_fCa in component ical (dimensionless)"
    legend_algebraic[70] = "i_CaL in component ical (A_per_F)"
    legend_algebraic[67] = "i_CaL_K in component ical (A_per_F)"
    legend_algebraic[69] = "i_CaL_Na in component ical (A_per_F)"
    legend_algebraic[64] = "ibarca in component ical (A_per_F)"
    legend_algebraic[66] = "ibark in component ical (A_per_F)"
    legend_algebraic[68] = "ibarna in component ical (A_per_F)"
    legend_algebraic[1] = "ical_d_a in component ical (mS_per_uF)"
    legend_algebraic[17] = "ical_d_b in component ical (mS_per_uF)"
    legend_algebraic[33] = "ical_d_inf in component ical (dimensionless)"
    legend_algebraic[47] = "ical_d_tau in component ical (ms)"
    legend_algebraic[2] = "ical_f_a in component ical (mS_per_uF)"
    legend_algebraic[18] = "ical_f_b in component ical (mS_per_uF)"
    legend_algebraic[34] = "ical_f_inf in component ical (dimensionless)"
    legend_algebraic[48] = "ical_f_tau in component ical (ms)"
    legend_algebraic[61] = "k_fca in component ical (dimensionless)"
    legend_constants[18] = "p_CaL in component ical (L_per_F_per_ms_times_1e0)"
    legend_constants[145] = "p_CaL_Ca in component ical (L_per_F_per_ms_times_1e0)"
    legend_constants[146] = "p_CaL_K in component ical (L_per_F_per_ms_times_1e0)"
    legend_constants[147] = "p_CaL_Na in component ical (L_per_F_per_ms_times_1e0)"
    legend_constants[19] = "p_CaL_shannonCa in component ical (dimensionless)"
    legend_constants[130] = "p_CaL_shannonCap in component ical (dimensionless)"
    legend_constants[20] = "p_CaL_shannonK in component ical (dimensionless)"
    legend_constants[131] = "p_CaL_shannonKp in component ical (dimensionless)"
    legend_constants[21] = "p_CaL_shannonNa in component ical (dimensionless)"
    legend_constants[132] = "p_CaL_shannonNap in component ical (dimensionless)"
    legend_constants[95] = "p_CaL_shannonTot in component ical (dimensionless)"
    legend_constants[22] = "scale in component ical (dimensionless)"
    legend_constants[23] = "tau_fCa in component ical (ms)"
    legend_constants[24] = "taud_const in component ical (ms)"
    legend_constants[25] = "tauf_const in component ical (ms)"
    legend_states[8] = "d in component icat (dimensionless)"
    legend_states[9] = "f in component icat (dimensionless)"
    legend_constants[26] = "g_CaT in component icat (mS_per_uF)"
    legend_algebraic[4] = "icat_d_inf in component icat (dimensionless)"
    legend_algebraic[20] = "icat_d_tau in component icat (ms)"
    legend_algebraic[5] = "icat_f_inf in component icat (dimensionless)"
    legend_algebraic[21] = "icat_f_tau in component icat (ms)"
    legend_constants[96] = "Na_frac in component ifunny (dimensionless)"
    legend_constants[27] = "NatoK_ratio in component ifunny (dimensionless)"
    legend_states[10] = "Xf in component ifunny (dimensionless)"
    legend_constants[28] = "g_f in component ifunny (mS_per_uF)"
    legend_algebraic[74] = "i_f in component ifunny (A_per_F)"
    legend_algebraic[72] = "i_fK in component ifunny (A_per_F)"
    legend_algebraic[73] = "i_fNa in component ifunny (A_per_F)"
    legend_algebraic[6] = "ifunny_Xf_a in component ifunny (mS_per_uF)"
    legend_algebraic[22] = "ifunny_Xf_b in component ifunny (mS_per_uF)"
    legend_algebraic[36] = "ifunny_Xf_inf in component ifunny (dimensionless)"
    legend_algebraic[50] = "ifunny_Xf_tau in component ifunny (ms)"
    legend_constants[29] = "xF1 in component ifunny (mS_per_uF)"
    legend_constants[97] = "xF2 in component ifunny (mV)"
    legend_constants[98] = "xF3 in component ifunny (mS_per_uF)"
    legend_constants[133] = "xF4 in component ifunny (mV)"
    legend_constants[30] = "xF5 in component ifunny (dimensionless)"
    legend_constants[31] = "xF6 in component ifunny (mV)"
    legend_constants[32] = "xF_const in component ifunny (ms)"
    legend_constants[33] = "g_K1 in component ik1 (mS_per_uF)"
    legend_algebraic[78] = "i_K1 in component ik1 (A_per_F)"
    legend_algebraic[75] = "ik1_inf_a in component ik1 (mS_per_uF)"
    legend_algebraic[76] = "ik1_inf_b in component ik1 (mS_per_uF)"
    legend_algebraic[77] = "inf in component ik1 (dimensionless)"
    legend_constants[34] = "xK11 in component ik1 (mS_per_uF)"
    legend_constants[35] = "xK12 in component ik1 (mV)"
    legend_constants[36] = "xK13 in component ik1 (mV)"
    legend_constants[37] = "xK14 in component ik1 (mV)"
    legend_constants[38] = "xK15 in component ik1 (mV)"
    legend_states[11] = "Xr1 in component ikr (dimensionless)"
    legend_constants[39] = "Xr1_1 in component ikr (mS_per_uF)"
    legend_constants[40] = "Xr1_2 in component ikr (mV)"
    legend_constants[99] = "Xr1_3 in component ikr (mS_per_uF)"
    legend_constants[134] = "Xr1_4 in component ikr (mV)"
    legend_constants[41] = "Xr1_5 in component ikr (dimensionless)"
    legend_constants[100] = "Xr1_6 in component ikr (mV)"
    legend_states[12] = "Xr2 in component ikr (dimensionless)"
    legend_constants[42] = "Xr2_1 in component ikr (mS_per_uF)"
    legend_constants[101] = "Xr2_2 in component ikr (mV)"
    legend_constants[102] = "Xr2_3 in component ikr (mS_per_uF)"
    legend_constants[135] = "Xr2_4 in component ikr (mV)"
    legend_constants[43] = "Xr2_5 in component ikr (dimensionless)"
    legend_constants[44] = "Xr2_6 in component ikr (mV)"
    legend_constants[45] = "g_Kr in component ikr (mS_per_uF)"
    legend_algebraic[79] = "i_Kr in component ikr (A_per_F)"
    legend_algebraic[7] = "ikr_Xr1_a in component ikr (mS_per_uF)"
    legend_algebraic[23] = "ikr_Xr1_b in component ikr (mS_per_uF)"
    legend_algebraic[37] = "ikr_Xr1_inf in component ikr (dimensionless)"
    legend_algebraic[51] = "ikr_Xr1_tau in component ikr (ms)"
    legend_algebraic[8] = "ikr_Xr2_a in component ikr (mS_per_uF)"
    legend_algebraic[24] = "ikr_Xr2_b in component ikr (mS_per_uF)"
    legend_algebraic[38] = "ikr_Xr2_inf in component ikr (dimensionless)"
    legend_algebraic[52] = "ikr_Xr2_tau in component ikr (ms)"
    legend_constants[46] = "tau_1_offset in component ikr (ms)"
    legend_constants[47] = "tau_2_offset in component ikr (ms)"
    legend_states[13] = "Xs in component iks (dimensionless)"
    legend_constants[48] = "g_Ks in component iks (mS_per_uF)"
    legend_algebraic[80] = "i_Ks in component iks (A_per_F)"
    legend_algebraic[9] = "iks_Xs_a in component iks (mS_per_uF)"
    legend_algebraic[25] = "iks_Xs_b in component iks (mS_per_uF)"
    legend_algebraic[39] = "iks_Xs_inf in component iks (dimensionless)"
    legend_algebraic[53] = "iks_Xs_tau in component iks (ms)"
    legend_constants[49] = "ks1 in component iks (mS_per_uF)"
    legend_constants[50] = "ks2 in component iks (mV)"
    legend_constants[103] = "ks3 in component iks (mS_per_uF)"
    legend_constants[136] = "ks4 in component iks (mV)"
    legend_constants[51] = "ks5 in component iks (dimensionless)"
    legend_constants[104] = "ks6 in component iks (mV)"
    legend_constants[52] = "tauks_const in component iks (ms)"
    legend_constants[105] = "V_leak in component ileak (mS_per_uF)"
    legend_constants[53] = "g_Na in component ina (mS_per_uF)"
    legend_states[14] = "h in component ina (dimensionless)"
    legend_constants[54] = "h1 in component ina (mS_per_uF)"
    legend_constants[106] = "h2 in component ina (mV)"
    legend_constants[107] = "h3 in component ina (mS_per_uF)"
    legend_constants[137] = "h4 in component ina (mV)"
    legend_constants[55] = "h5 in component ina (dimensionless)"
    legend_constants[56] = "h6 in component ina (mV)"
    legend_algebraic[82] = "i_Na in component ina (A_per_F)"
    legend_algebraic[10] = "ina_h_a in component ina (mS_per_uF)"
    legend_algebraic[26] = "ina_h_b in component ina (mS_per_uF)"
    legend_algebraic[40] = "ina_h_inf in component ina (dimensionless)"
    legend_algebraic[54] = "ina_h_tau in component ina (ms)"
    legend_algebraic[11] = "ina_j_a in component ina (mS_per_uF)"
    legend_algebraic[27] = "ina_j_b in component ina (mS_per_uF)"
    legend_algebraic[41] = "ina_j_inf in component ina (dimensionless)"
    legend_algebraic[55] = "ina_j_tau in component ina (ms)"
    legend_algebraic[12] = "ina_m_a in component ina (mS_per_uF)"
    legend_algebraic[28] = "ina_m_b in component ina (mS_per_uF)"
    legend_algebraic[42] = "ina_m_inf in component ina (dimensionless)"
    legend_algebraic[56] = "ina_m_tau in component ina (ms)"
    legend_states[15] = "j in component ina (dimensionless)"
    legend_constants[57] = "j1 in component ina (mS_per_uF)"
    legend_constants[108] = "j2 in component ina (mV)"
    legend_constants[138] = "j3 in component ina (mS_per_uF)"
    legend_constants[148] = "j4 in component ina (mV)"
    legend_constants[109] = "j5 in component ina (dimensionless)"
    legend_constants[139] = "j6 in component ina (mV)"
    legend_states[16] = "m in component ina (dimensionless)"
    legend_constants[58] = "m1 in component ina (mS_per_uF)"
    legend_constants[59] = "m2 in component ina (mV)"
    legend_constants[110] = "m3 in component ina (mS_per_uF)"
    legend_constants[140] = "m4 in component ina (mV)"
    legend_constants[60] = "m5 in component ina (dimensionless)"
    legend_constants[111] = "m6 in component ina (mV)"
    legend_constants[61] = "tau_h_const in component ina (ms)"
    legend_constants[62] = "tau_j_const in component ina (ms)"
    legend_constants[63] = "tau_m_const in component ina (ms)"
    legend_constants[64] = "KmCa in component inaca (mM)"
    legend_constants[65] = "KmNai in component inaca (mM)"
    legend_constants[66] = "Ksat in component inaca (dimensionless)"
    legend_constants[112] = "alpha in component inaca (dimensionless)"
    legend_constants[113] = "gamma in component inaca (dimensionless)"
    legend_constants[114] = "kNaCa in component inaca (A_per_F)"
    legend_constants[67] = "Km_K in component inak (mM)"
    legend_constants[68] = "Km_Na in component inak (mM)"
    legend_constants[115] = "PNaK in component inak (A_per_F)"
    legend_algebraic[84] = "i_NaK in component inak (A_per_F)"
    legend_constants[69] = "KPCa in component ipca (mM)"
    legend_constants[116] = "g_PCa in component ipca (A_per_F)"
    legend_states[17] = "I in component irel (dimensionless)"
    legend_constants[70] = "MaxSR in component irel (dimensionless)"
    legend_constants[71] = "MinSR in component irel (dimensionless)"
    legend_states[18] = "O in component irel (dimensionless)"
    legend_states[19] = "R in component irel (dimensionless)"
    legend_algebraic[13] = "RI in component irel (dimensionless)"
    legend_constants[72] = "ec50SR in component irel (mM)"
    legend_algebraic[29] = "kCaSR in component irel (dimensionless)"
    legend_constants[117] = "kiCa in component irel (per_mM_per_ms)"
    legend_algebraic[43] = "kiSRCa in component irel (per_mM_per_ms)"
    legend_constants[118] = "kim in component irel (mS_per_uF)"
    legend_constants[141] = "koCa in component irel (per_mM2_per_ms)"
    legend_algebraic[57] = "koSRCa in component irel (per_mM2_per_ms)"
    legend_constants[119] = "kom in component irel (mS_per_uF)"
    legend_constants[73] = "ks in component irel (mS_per_uF)"
    legend_constants[74] = "g_to in component ito (mS_per_uF)"
    legend_algebraic[87] = "i_to in component ito (A_per_F)"
    legend_algebraic[14] = "ito_r_a in component ito (mS_per_uF)"
    legend_algebraic[30] = "ito_r_b in component ito (mS_per_uF)"
    legend_algebraic[44] = "ito_r_inf in component ito (dimensionless)"
    legend_algebraic[58] = "ito_r_tau in component ito (ms)"
    legend_algebraic[15] = "ito_s_a in component ito (mS_per_uF)"
    legend_algebraic[31] = "ito_s_b in component ito (mS_per_uF)"
    legend_algebraic[45] = "ito_s_inf in component ito (dimensionless)"
    legend_algebraic[59] = "ito_s_tau in component ito (ms)"
    legend_states[20] = "r in component ito (dimensionless)"
    legend_constants[75] = "r1 in component ito (mS_per_uF)"
    legend_constants[76] = "r2 in component ito (mV)"
    legend_constants[120] = "r3 in component ito (mS_per_uF)"
    legend_constants[142] = "r4 in component ito (mV)"
    legend_constants[77] = "r5 in component ito (dimensionless)"
    legend_constants[121] = "r6 in component ito (mV)"
    legend_states[21] = "s in component ito (dimensionless)"
    legend_constants[78] = "s1 in component ito (mS_per_uF)"
    legend_constants[122] = "s2 in component ito (mV)"
    legend_constants[123] = "s3 in component ito (mS_per_uF)"
    legend_constants[143] = "s4 in component ito (mV)"
    legend_constants[79] = "s5 in component ito (dimensionless)"
    legend_constants[80] = "s6 in component ito (mV)"
    legend_constants[81] = "tau_r_const in component ito (ms)"
    legend_constants[82] = "tau_s_const in component ito (ms)"
    legend_constants[124] = "Kup in component iup (mM)"
    legend_constants[125] = "VmaxUp in component iup (mM_per_ms)"
    legend_constants[150] = "i_stim in component stimulus (A_per_F)"
    legend_constants[83] = "R in component phys (J_per_mol_per_K)"
    legend_constants[84] = "T in component phys (kelvin)"
    legend_constants[149] = "amplitude in component stimulus (A_per_F)"
    legend_rates[0] = "d/dt Cai in component cai (mM)"
    legend_rates[1] = "d/dt Ca_SR in component casr (mM)"
    legend_rates[5] = "d/dt d in component ical (dimensionless)"
    legend_rates[6] = "d/dt f in component ical (dimensionless)"
    legend_rates[7] = "d/dt fCa in component ical (dimensionless)"
    legend_rates[8] = "d/dt d in component icat (dimensionless)"
    legend_rates[9] = "d/dt f in component icat (dimensionless)"
    legend_rates[10] = "d/dt Xf in component ifunny (dimensionless)"
    legend_rates[11] = "d/dt Xr1 in component ikr (dimensionless)"
    legend_rates[12] = "d/dt Xr2 in component ikr (dimensionless)"
    legend_rates[13] = "d/dt Xs in component iks (dimensionless)"
    legend_rates[14] = "d/dt h in component ina (dimensionless)"
    legend_rates[15] = "d/dt j in component ina (dimensionless)"
    legend_rates[16] = "d/dt m in component ina (dimensionless)"
    legend_rates[17] = "d/dt I in component irel (dimensionless)"
    legend_rates[18] = "d/dt O in component irel (dimensionless)"
    legend_rates[19] = "d/dt R in component irel (dimensionless)"
    legend_rates[20] = "d/dt r in component ito (dimensionless)"
    legend_rates[21] = "d/dt s in component ito (dimensionless)"
    legend_rates[2] = "d/dt Ki in component ki (mM)"
    legend_rates[4] = "d/dt V in component membrane (mV)"
    legend_rates[3] = "d/dt Nai in component nai (mM)"
    return (legend_states, legend_algebraic, legend_voi, legend_constants)

def initConsts():
    constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
    constants[0] = 0.06
    states[0] = 2.19191642424964390e-04
    constants[1] = 60.0
    constants[2] = 9.64853415000000041e+01
    constants[3] = 0.0006
    states[1] = 3.35086796732326109e-01
    constants[4] = 0.3
    constants[5] = 0.0
    constants[6] = 1.8
    states[2] = 1.04748824394112106e+02
    constants[7] = 5.4
    states[3] = 7.16928091250999167
    constants[8] = 140.0
    constants[9] = 1094.0
    constants[10] = 3960.0
    constants[11] = 16404.0
    states[4] = -7.55966016388546791e+01
    states[5] = 3.94925342652924281e-04
    constants[12] = 1.29662941897219994e+01
    constants[13] = 7.07914596471100044
    constants[14] = 4.49094155069999987e-02
    states[6] = 1.70990105585540286e-01
    constants[15] = 5.12589825999999987e-04
    constants[16] = 1.93121122351431995e+03
    constants[17] = 5.73002749969900016
    states[7] = 8.77798946134088598e-01
    constants[18] = 3.08027691378999990e-01
    constants[19] = 0.00054
    constants[20] = 2.7e-07
    constants[21] = 1.5e-08
    constants[22] = 1.2
    constants[23] = 2.0
    constants[24] = 1.65824694683000007
    constants[25] = 1.00462559171102995e+02
    states[8] = 2.70195573471577175e-04
    states[9] = 7.56032904368393432e-01
    constants[26] = 0.185
    constants[27] = 0.491
    states[10] = 6.40338504912615469e-03
    constants[28] = 0.0435
    constants[29] = 5.78970000000000002e-07
    constants[30] = 2.00866502378844016e+04
    constants[31] = 1.02023528452800001e+01
    constants[32] = 2.39452913465299986e+01
    constants[33] = 1.33785777797606004e-01
    constants[34] = 4.77994972217041014e-01
    constants[35] = 2.72427558793486995e+01
    constants[36] = 4.92502331781412028
    constants[37] = 8.72223760006881932
    constants[38] = 5.66361974998243980e+01
    states[11] = 3.09767485715433222e-01
    constants[39] = 5.74885237435000026e-03
    constants[40] = 1.36234926362576001e+01
    constants[41] = 4.76305711818360011e-02
    states[12] = 4.50577185148518577e-01
    constants[42] = 1.24566405268270002e-02
    constants[43] = 3.73426331501040991e+01
    constants[44] = 2.20919642353902006e+01
    constants[45] = 0.218025
    constants[46] = 50.0
    constants[47] = 0.0
    states[13] = 1.53788281650948710e-01
    constants[48] = 0.0077
    constants[49] = 1.16558447999999992e-03
    constants[50] = 6.67268386758935958e+04
    constants[51] = 2.80458908250000027e-01
    constants[52] = 4.74115000000000034e-06
    constants[53] = 9.72061340924100037
    states[14] = 7.39543607812429227e-01
    constants[54] = 3.62659886399999999e-03
    constants[55] = 9.66329497711473959e+03
    constants[56] = 7.39550356461299963
    states[15] = 1.24515982574504899e-01
    constants[57] = 5.12257182000000044e-04
    states[16] = 2.97549962926413614e-02
    constants[58] = 1.08045846384818006e+02
    constants[59] = 1.31070157339409992e+01
    constants[60] = 2.32691436700000007e-03
    constants[61] = 1.67331502516000014e-01
    constants[62] = 9.51088724962000032e-01
    constants[63] = 3.19775803839999970e-02
    constants[64] = 1.38
    constants[65] = 87.5
    constants[66] = 0.1
    constants[67] = 1.0
    constants[68] = 40.0
    constants[69] = 0.0005
    states[17] = 1.42153622323011597e-02
    constants[70] = 15.0
    constants[71] = 1.0
    states[18] = 1.65045105312396393e-04
    states[19] = 1.13120363433751106e-02
    constants[72] = 0.45
    constants[73] = 12.5
    constants[74] = 1.17833333333299997e-01
    states[20] = 2.67597833344160611e-04
    constants[75] = 5.53614181712999975e-02
    constants[76] = 1.16842023429669002e+01
    constants[77] = 3.98918108037750008
    states[21] = 7.46802810614006107e-01
    constants[78] = 3.44230944300000013e-04
    constants[79] = 1.86760536909694991e+02
    constants[80] = 8.18093387332270083
    constants[81] = 6.96758421171499998e-01
    constants[82] = 1.12244577239468999e+01
    constants[83] = 8.314472
    constants[84] = 310.0
    constants[85] = constants[11]+constants[9]
    constants[86] = 10.0000*1.20000
    constants[87] = (constants[83]*constants[84])/constants[2]
    constants[88] = 0.000592000*0.620000
    constants[89] = 0.000290000*1.50000
    constants[90] = constants[2]/(constants[83]*constants[84])
    constants[91] = constants[14]*constants[12]
    constants[92] = -6.90988
    constants[93] = -49.5057
    constants[94] = constants[16]*constants[15]
    constants[95] = (constants[19]+constants[21])+constants[20]
    constants[96] = constants[27]/(constants[27]+1.00000)
    constants[97] = -14.5897
    constants[98] = constants[30]*constants[29]
    constants[99] = constants[41]*constants[39]
    constants[100] = -7.06809
    constants[101] = -25.9945
    constants[102] = constants[43]*constants[42]
    constants[103] = constants[51]*constants[49]
    constants[104] = -18.8670
    constants[105] = 8.00000e-05*0.0200000
    constants[106] = -19.8394
    constants[107] = constants[55]*constants[54]
    constants[108] = -66.5838
    constants[109] = constants[55]
    constants[110] = constants[60]*constants[58]
    constants[111] = -7.91773
    constants[112] = 2.50000*1.10000
    constants[113] = 0.350000*2.00000
    constants[114] = 1000.00*1.10000
    constants[115] = 1.36200*1.81800
    constants[116] = 0.0250000*10.5000
    constants[117] = 54.0000*0.342500
    constants[118] = 0.00100000*0.557100
    constants[119] = 1.50000*0.142900
    constants[120] = constants[77]*constants[75]
    constants[121] = -11.0471
    constants[122] = -17.6345
    constants[123] = constants[79]*constants[78]
    constants[124] = 0.000250000*0.702000
    constants[125] = 0.000425000*0.260000
    constants[126] = constants[10]*(constants[9]/constants[85])
    constants[127] = constants[2]*constants[90]
    constants[128] = 1.00000/(1.00000/constants[13]+1.00000/constants[92])
    constants[129] = 1.00000/(1.00000/constants[93]+1.00000/constants[17])
    constants[130] = constants[19]/constants[95]
    constants[131] = constants[20]/constants[95]
    constants[132] = constants[21]/constants[95]
    constants[133] = 1.00000/(1.00000/constants[97]+1.00000/constants[31])
    constants[134] = 1.00000/(1.00000/constants[40]+1.00000/constants[100])
    constants[135] = 1.00000/(1.00000/constants[101]+1.00000/constants[44])
    constants[136] = 1.00000/(1.00000/constants[50]+1.00000/constants[104])
    constants[137] = 1.00000/(1.00000/constants[106]+1.00000/constants[56])
    constants[138] = constants[109]*constants[57]
    constants[139] = constants[56]
    constants[140] = 1.00000/(1.00000/constants[59]+1.00000/constants[111])
    constants[141] = 56320.0*11.4302
    constants[142] = 1.00000/(1.00000/constants[76]+1.00000/constants[121])
    constants[143] = 1.00000/(1.00000/constants[122]+1.00000/constants[80])
    constants[144] = constants[10]*(constants[11]/constants[85])
    constants[145] = constants[130]*constants[18]
    constants[146] = constants[131]*constants[18]
    constants[147] = constants[132]*constants[18]
    constants[148] = 1.00000/(1.00000/constants[108]+1.00000/constants[139])
    constants[149] = -3.00000
    constants[150] = constants[5]*constants[149]
    return (states, constants)

def computeRates(voi, states, constants):
    rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
    algebraic[4] = 1.00000/(1.00000+exp((states[4]+26.3000)/-6.00000))
    algebraic[20] = 1.00000/(1.06800*exp((states[4]+26.3000)/30.0000)+1.06800*exp((states[4]+26.3000)/-30.0000))
    rates[8] = (algebraic[4]-states[8])/algebraic[20]
    algebraic[5] = 1.00000/(1.00000+exp((states[4]+61.7000)/5.60000))
    algebraic[21] = 1.00000/(0.0153000*exp(-(states[4]+61.7000)/83.3000)+0.0150000*exp((states[4]+61.7000)/15.3800))
    rates[9] = (algebraic[5]-states[9])/algebraic[21]
    algebraic[1] = constants[12]*exp(states[4]/constants[13])
    algebraic[17] = constants[91]*exp(states[4]/constants[128])
    algebraic[33] = algebraic[1]/(algebraic[1]+algebraic[17])
    algebraic[47] = 1.00000/(algebraic[1]+algebraic[17])+constants[24]
    rates[5] = (algebraic[33]-states[5])/algebraic[47]
    algebraic[2] = constants[15]*exp(states[4]/constants[93])
    algebraic[18] = constants[94]*exp(states[4]/constants[129])
    algebraic[34] = algebraic[2]/(algebraic[2]+algebraic[18])
    algebraic[48] = 1.00000/(algebraic[2]+algebraic[18])+constants[25]
    rates[6] = (algebraic[34]-states[6])/algebraic[48]
    algebraic[6] = constants[29]*exp(states[4]/constants[97])
    algebraic[22] = constants[98]*exp(states[4]/constants[133])
    algebraic[36] = algebraic[6]/(algebraic[6]+algebraic[22])
    algebraic[50] = 1.00000/(algebraic[6]+algebraic[22])+constants[32]
    rates[10] = (algebraic[36]-states[10])/algebraic[50]
    algebraic[7] = constants[39]*exp(states[4]/constants[40])
    algebraic[23] = constants[99]*exp(states[4]/constants[134])
    algebraic[37] = algebraic[7]/(algebraic[7]+algebraic[23])
    algebraic[51] = 1.00000/(algebraic[7]+algebraic[23])+constants[46]
    rates[11] = (algebraic[37]-states[11])/algebraic[51]
    algebraic[8] = constants[42]*exp(states[4]/constants[101])
    algebraic[24] = constants[102]*exp(states[4]/constants[135])
    algebraic[38] = algebraic[8]/(algebraic[8]+algebraic[24])
    algebraic[52] = 1.00000/(algebraic[8]+algebraic[24])+constants[47]
    rates[12] = (algebraic[38]-states[12])/algebraic[52]
    algebraic[9] = constants[49]*exp(states[4]/constants[50])
    algebraic[25] = constants[103]*exp(states[4]/constants[136])
    algebraic[39] = algebraic[9]/(algebraic[9]+algebraic[25])
    algebraic[53] = 1.00000/(algebraic[9]+algebraic[25])+constants[52]
    rates[13] = (algebraic[39]-states[13])/algebraic[53]
    algebraic[10] = constants[54]*exp(states[4]/constants[106])
    algebraic[26] = constants[107]*exp(states[4]/constants[137])
    algebraic[40] = algebraic[10]/(algebraic[10]+algebraic[26])
    algebraic[54] = 1.00000/(algebraic[10]+algebraic[26])+constants[61]
    rates[14] = (algebraic[40]-states[14])/algebraic[54]
    algebraic[11] = constants[57]*exp(states[4]/constants[108])
    algebraic[27] = constants[138]*exp(states[4]/constants[148])
    algebraic[41] = algebraic[11]/(algebraic[11]+algebraic[27])
    algebraic[55] = 1.00000/(algebraic[11]+algebraic[27])+constants[62]
    rates[15] = (algebraic[41]-states[15])/algebraic[55]
    algebraic[12] = constants[58]*exp(states[4]/constants[59])
    algebraic[28] = constants[110]*exp(states[4]/constants[140])
    algebraic[42] = algebraic[12]/(algebraic[12]+algebraic[28])
    algebraic[56] = 1.00000/(algebraic[12]+algebraic[28])+constants[63]
    rates[16] = (algebraic[42]-states[16])/algebraic[56]
    algebraic[13] = ((1.00000-states[19])-states[18])-states[17]
    algebraic[29] = constants[70]-(constants[70]-constants[71])/(1.00000+power(constants[72]/states[1], 2.50000))
    algebraic[43] = constants[117]*algebraic[29]
    algebraic[57] = constants[141]/algebraic[29]
    rates[17] = (((algebraic[43]*states[0])*states[18]-constants[118]*states[17])-constants[119]*states[17])+(algebraic[57]*(power(states[0], 2.00000)))*algebraic[13]
    rates[18] = (((algebraic[57]*(power(states[0], 2.00000)))*states[19]-constants[119]*states[18])-(algebraic[43]*states[0])*states[18])+constants[118]*states[17]
    rates[19] = ((constants[118]*algebraic[13]-(algebraic[43]*states[0])*states[19])-(algebraic[57]*(power(states[0], 2.00000)))*states[19])+constants[119]*states[18]
    algebraic[14] = constants[75]*exp(states[4]/constants[76])
    algebraic[30] = constants[120]*exp(states[4]/constants[142])
    algebraic[44] = algebraic[14]/(algebraic[14]+algebraic[30])
    algebraic[58] = 1.00000/(algebraic[14]+algebraic[30])+constants[81]
    rates[20] = (algebraic[44]-states[20])/algebraic[58]
    algebraic[15] = constants[78]*exp(states[4]/constants[122])
    algebraic[31] = constants[123]*exp(states[4]/constants[143])
    algebraic[45] = algebraic[15]/(algebraic[15]+algebraic[31])
    algebraic[59] = 1.00000/(algebraic[15]+algebraic[31])+constants[82]
    rates[21] = (algebraic[45]-states[21])/algebraic[59]
    algebraic[3] = 1.00000/(1.00000+power((constants[22]*states[0])/0.000325000, 8.00000))
    algebraic[19] = 0.100000/(1.00000+exp((constants[22]*states[0]-0.000500000)/0.000100000))
    algebraic[35] = 0.200000/(1.00000+exp((constants[22]*states[0]-0.000750000)/0.000800000))
    algebraic[49] = (((algebraic[3]+algebraic[19])+algebraic[35])+0.230000)/1.46000
    algebraic[61] = custom_piecewise([greater(algebraic[49] , states[7]) & greater(states[4] , -60.0000), 0.00000 , True, 1.00000])
    rates[7] = (algebraic[61]*(algebraic[49]-states[7]))/constants[23]
    algebraic[83] = (constants[114]*((exp((constants[113]*states[4])*constants[90])*(power(states[3], 3.00000)))*constants[6]-((exp(((constants[113]-1.00000)*states[4])*constants[90])*(power(constants[8], 3.00000)))*states[0])*constants[112]))/(((power(constants[65], 3.00000)+power(constants[8], 3.00000))*(constants[64]+constants[6]))*(1.00000+constants[66]*exp(((constants[113]-1.00000)*states[4])*constants[90])))
    algebraic[60] = constants[87]*log(constants[8]/states[3])
    algebraic[63] = constants[89]*(states[4]-algebraic[60])
    algebraic[68] = (((constants[147]*states[4])*constants[127])*((0.750000*states[3])*exp(states[4]*constants[90])-0.750000*constants[8]))/(exp(states[4]*constants[90])-1.00000)
    algebraic[69] = ((algebraic[68]*states[5])*states[6])*states[7]
    algebraic[73] = ((constants[96]*constants[28])*states[10])*(states[4]-algebraic[60])
    algebraic[82] = (((constants[53]*(power(states[16], 3.00000)))*states[14])*states[15])*(states[4]-algebraic[60])
    algebraic[84] = ((constants[115]*constants[7])*states[3])/(((constants[7]+constants[67])*(states[3]+constants[68]))*((1.00000+0.124500*exp((-0.100000*states[4])*constants[90]))+0.0353000*exp(-states[4]*constants[90])))
    rates[3] = (-constants[1]/(constants[2]*constants[144]))*((((((algebraic[82])+algebraic[63])+algebraic[73])+3.00000*algebraic[84])+3.00000*algebraic[83])+algebraic[69])
    algebraic[66] = (((constants[146]*states[4])*constants[127])*((0.750000*states[2])*exp(states[4]*constants[90])-0.750000*constants[7]))/(exp(states[4]*constants[90])-1.00000)
    algebraic[67] = ((algebraic[66]*states[5])*states[6])*states[7]
    algebraic[46] = constants[87]*log(constants[7]/states[2])
    algebraic[72] = (((1.00000-constants[96])*constants[28])*states[10])*(states[4]-algebraic[46])
    algebraic[75] = constants[34]*exp((states[4]+constants[36])/constants[35])
    algebraic[76] = 1.00000*exp((states[4]+constants[38])/constants[37])
    algebraic[77] = algebraic[75]/(algebraic[75]+algebraic[76])
    algebraic[78] = ((constants[33]*(power(constants[7]/5.40000, 1.0/2)))*algebraic[77])*(states[4]-algebraic[46])
    algebraic[79] = (((constants[45]*(power(constants[7]/5.40000, 1.0/2)))*states[11])*states[12])*(states[4]-algebraic[46])
    algebraic[80] = (constants[48]*(power(states[13], 2.00000)))*(states[4]-algebraic[46])
    algebraic[87] = ((constants[74]*states[20])*states[21])*(states[4]-algebraic[46])
    rates[2] = (-constants[1]/(constants[2]*constants[144]))*(((((((algebraic[78])+algebraic[87])+algebraic[79])+algebraic[80])+algebraic[72])-2.00000*algebraic[84])+algebraic[67])
    algebraic[32] = (0.500000*constants[87])*log(constants[6]/states[0])
    algebraic[71] = ((constants[26]*states[8])*states[9])*(states[4]-algebraic[32])
    algebraic[85] = (constants[116]*states[0])/(states[0]+constants[69])
    algebraic[62] = constants[88]*(states[4]-algebraic[32])
    algebraic[64] = ((((constants[145]*4.00000)*states[4])*constants[127])*((0.341000*states[0])*exp((2.00000*states[4])*constants[90])-0.341000*constants[6]))/(exp((2.00000*states[4])*constants[90])-1.00000)
    algebraic[65] = ((algebraic[64]*states[5])*states[6])*states[7]
    algebraic[70] = (algebraic[65]+algebraic[69])+algebraic[67]
    algebraic[74] = algebraic[73]+algebraic[72]
    rates[4] = -((((((((((((((algebraic[78])+algebraic[87])+algebraic[79])+algebraic[80])+algebraic[70])+algebraic[71])+algebraic[84])+algebraic[82])+algebraic[83])+algebraic[85])+algebraic[74])+algebraic[63])+algebraic[62])+constants[150])
    algebraic[0] = 1.00000/(1.00000+(constants[0]*constants[3])/(power(states[0]+constants[3], 2.00000)))
    algebraic[81] = (states[1]-states[0])*constants[105]
    algebraic[86] = ((constants[73]*states[18])*(states[1]-states[0]))*(constants[126]/constants[144])
    algebraic[88] = constants[125]/(1.00000+(power(constants[124], 2.00000))/(power(states[0], 2.00000)))
    rates[0] = algebraic[0]*(((-algebraic[88]+algebraic[81])+algebraic[86])-(constants[1]/((2.00000*constants[144])*constants[2]))*(((((algebraic[65])+algebraic[71])+algebraic[62])+algebraic[85])-2.00000*algebraic[83]))
    algebraic[16] = 1.00000/(1.00000+(constants[86]*constants[4])/(power(states[1]+constants[4], 2.00000)))
    rates[1] = ((algebraic[16]*constants[144])/constants[126])*((algebraic[88]-algebraic[86])-algebraic[81])
    return(rates)

def computeAlgebraic(constants, states, voi):
    algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
    states = array(states)
    voi = array(voi)
    algebraic[4] = 1.00000/(1.00000+exp((states[4]+26.3000)/-6.00000))
    algebraic[20] = 1.00000/(1.06800*exp((states[4]+26.3000)/30.0000)+1.06800*exp((states[4]+26.3000)/-30.0000))
    algebraic[5] = 1.00000/(1.00000+exp((states[4]+61.7000)/5.60000))
    algebraic[21] = 1.00000/(0.0153000*exp(-(states[4]+61.7000)/83.3000)+0.0150000*exp((states[4]+61.7000)/15.3800))
    algebraic[1] = constants[12]*exp(states[4]/constants[13])
    algebraic[17] = constants[91]*exp(states[4]/constants[128])
    algebraic[33] = algebraic[1]/(algebraic[1]+algebraic[17])
    algebraic[47] = 1.00000/(algebraic[1]+algebraic[17])+constants[24]
    algebraic[2] = constants[15]*exp(states[4]/constants[93])
    algebraic[18] = constants[94]*exp(states[4]/constants[129])
    algebraic[34] = algebraic[2]/(algebraic[2]+algebraic[18])
    algebraic[48] = 1.00000/(algebraic[2]+algebraic[18])+constants[25]
    algebraic[6] = constants[29]*exp(states[4]/constants[97])
    algebraic[22] = constants[98]*exp(states[4]/constants[133])
    algebraic[36] = algebraic[6]/(algebraic[6]+algebraic[22])
    algebraic[50] = 1.00000/(algebraic[6]+algebraic[22])+constants[32]
    algebraic[7] = constants[39]*exp(states[4]/constants[40])
    algebraic[23] = constants[99]*exp(states[4]/constants[134])
    algebraic[37] = algebraic[7]/(algebraic[7]+algebraic[23])
    algebraic[51] = 1.00000/(algebraic[7]+algebraic[23])+constants[46]
    algebraic[8] = constants[42]*exp(states[4]/constants[101])
    algebraic[24] = constants[102]*exp(states[4]/constants[135])
    algebraic[38] = algebraic[8]/(algebraic[8]+algebraic[24])
    algebraic[52] = 1.00000/(algebraic[8]+algebraic[24])+constants[47]
    algebraic[9] = constants[49]*exp(states[4]/constants[50])
    algebraic[25] = constants[103]*exp(states[4]/constants[136])
    algebraic[39] = algebraic[9]/(algebraic[9]+algebraic[25])
    algebraic[53] = 1.00000/(algebraic[9]+algebraic[25])+constants[52]
    algebraic[10] = constants[54]*exp(states[4]/constants[106])
    algebraic[26] = constants[107]*exp(states[4]/constants[137])
    algebraic[40] = algebraic[10]/(algebraic[10]+algebraic[26])
    algebraic[54] = 1.00000/(algebraic[10]+algebraic[26])+constants[61]
    algebraic[11] = constants[57]*exp(states[4]/constants[108])
    algebraic[27] = constants[138]*exp(states[4]/constants[148])
    algebraic[41] = algebraic[11]/(algebraic[11]+algebraic[27])
    algebraic[55] = 1.00000/(algebraic[11]+algebraic[27])+constants[62]
    algebraic[12] = constants[58]*exp(states[4]/constants[59])
    algebraic[28] = constants[110]*exp(states[4]/constants[140])
    algebraic[42] = algebraic[12]/(algebraic[12]+algebraic[28])
    algebraic[56] = 1.00000/(algebraic[12]+algebraic[28])+constants[63]
    algebraic[13] = ((1.00000-states[19])-states[18])-states[17]
    algebraic[29] = constants[70]-(constants[70]-constants[71])/(1.00000+power(constants[72]/states[1], 2.50000))
    algebraic[43] = constants[117]*algebraic[29]
    algebraic[57] = constants[141]/algebraic[29]
    algebraic[14] = constants[75]*exp(states[4]/constants[76])
    algebraic[30] = constants[120]*exp(states[4]/constants[142])
    algebraic[44] = algebraic[14]/(algebraic[14]+algebraic[30])
    algebraic[58] = 1.00000/(algebraic[14]+algebraic[30])+constants[81]
    algebraic[15] = constants[78]*exp(states[4]/constants[122])
    algebraic[31] = constants[123]*exp(states[4]/constants[143])
    algebraic[45] = algebraic[15]/(algebraic[15]+algebraic[31])
    algebraic[59] = 1.00000/(algebraic[15]+algebraic[31])+constants[82]
    algebraic[3] = 1.00000/(1.00000+power((constants[22]*states[0])/0.000325000, 8.00000))
    algebraic[19] = 0.100000/(1.00000+exp((constants[22]*states[0]-0.000500000)/0.000100000))
    algebraic[35] = 0.200000/(1.00000+exp((constants[22]*states[0]-0.000750000)/0.000800000))
    algebraic[49] = (((algebraic[3]+algebraic[19])+algebraic[35])+0.230000)/1.46000
    algebraic[61] = custom_piecewise([greater(algebraic[49] , states[7]) & greater(states[4] , -60.0000), 0.00000 , True, 1.00000])
    algebraic[83] = (constants[114]*((exp((constants[113]*states[4])*constants[90])*(power(states[3], 3.00000)))*constants[6]-((exp(((constants[113]-1.00000)*states[4])*constants[90])*(power(constants[8], 3.00000)))*states[0])*constants[112]))/(((power(constants[65], 3.00000)+power(constants[8], 3.00000))*(constants[64]+constants[6]))*(1.00000+constants[66]*exp(((constants[113]-1.00000)*states[4])*constants[90])))
    algebraic[60] = constants[87]*log(constants[8]/states[3])
    algebraic[63] = constants[89]*(states[4]-algebraic[60])
    algebraic[68] = (((constants[147]*states[4])*constants[127])*((0.750000*states[3])*exp(states[4]*constants[90])-0.750000*constants[8]))/(exp(states[4]*constants[90])-1.00000)
    algebraic[69] = ((algebraic[68]*states[5])*states[6])*states[7]
    algebraic[73] = ((constants[96]*constants[28])*states[10])*(states[4]-algebraic[60])
    algebraic[82] = (((constants[53]*(power(states[16], 3.00000)))*states[14])*states[15])*(states[4]-algebraic[60])
    algebraic[84] = ((constants[115]*constants[7])*states[3])/(((constants[7]+constants[67])*(states[3]+constants[68]))*((1.00000+0.124500*exp((-0.100000*states[4])*constants[90]))+0.0353000*exp(-states[4]*constants[90])))
    algebraic[66] = (((constants[146]*states[4])*constants[127])*((0.750000*states[2])*exp(states[4]*constants[90])-0.750000*constants[7]))/(exp(states[4]*constants[90])-1.00000)
    algebraic[67] = ((algebraic[66]*states[5])*states[6])*states[7]
    algebraic[46] = constants[87]*log(constants[7]/states[2])
    algebraic[72] = (((1.00000-constants[96])*constants[28])*states[10])*(states[4]-algebraic[46])
    algebraic[75] = constants[34]*exp((states[4]+constants[36])/constants[35])
    algebraic[76] = 1.00000*exp((states[4]+constants[38])/constants[37])
    algebraic[77] = algebraic[75]/(algebraic[75]+algebraic[76])
    algebraic[78] = ((constants[33]*(power(constants[7]/5.40000, 1.0/2)))*algebraic[77])*(states[4]-algebraic[46])
    algebraic[79] = (((constants[45]*(power(constants[7]/5.40000, 1.0/2)))*states[11])*states[12])*(states[4]-algebraic[46])
    algebraic[80] = (constants[48]*(power(states[13], 2.00000)))*(states[4]-algebraic[46])
    algebraic[87] = ((constants[74]*states[20])*states[21])*(states[4]-algebraic[46])
    algebraic[32] = (0.500000*constants[87])*log(constants[6]/states[0])
    algebraic[71] = ((constants[26]*states[8])*states[9])*(states[4]-algebraic[32])
    algebraic[85] = (constants[116]*states[0])/(states[0]+constants[69])
    algebraic[62] = constants[88]*(states[4]-algebraic[32])
    algebraic[64] = ((((constants[145]*4.00000)*states[4])*constants[127])*((0.341000*states[0])*exp((2.00000*states[4])*constants[90])-0.341000*constants[6]))/(exp((2.00000*states[4])*constants[90])-1.00000)
    algebraic[65] = ((algebraic[64]*states[5])*states[6])*states[7]
    algebraic[70] = (algebraic[65]+algebraic[69])+algebraic[67]
    algebraic[74] = algebraic[73]+algebraic[72]
    algebraic[0] = 1.00000/(1.00000+(constants[0]*constants[3])/(power(states[0]+constants[3], 2.00000)))
    algebraic[81] = (states[1]-states[0])*constants[105]
    algebraic[86] = ((constants[73]*states[18])*(states[1]-states[0]))*(constants[126]/constants[144])
    algebraic[88] = constants[125]/(1.00000+(power(constants[124], 2.00000))/(power(states[0], 2.00000)))
    algebraic[16] = 1.00000/(1.00000+(constants[86]*constants[4])/(power(states[1]+constants[4], 2.00000)))
    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)