# Size of variable arrays:
sizeAlgebraic = 71
sizeStates = 41
sizeConstants = 73
from math import *
from numpy import *

def createLegends():
    legend_states = [""] * sizeStates
    legend_rates = [""] * sizeStates
    legend_algebraic = [""] * sizeAlgebraic
    legend_voi = ""
    legend_constants = [""] * sizeConstants
    legend_voi = "time in component environment (millisecond)"
    legend_states[0] = "V in component membrane (millivolt)"
    legend_constants[0] = "Cm in component membrane (microF_per_cm2)"
    legend_constants[1] = "Vmyo in component membrane (microlitre)"
    legend_constants[2] = "VJSR in component membrane (microlitre)"
    legend_constants[3] = "VNSR in component membrane (microlitre)"
    legend_constants[4] = "Vss in component membrane (microlitre)"
    legend_constants[5] = "Acap in component membrane (cm2)"
    legend_constants[6] = "Ko in component membrane (micromolar)"
    legend_constants[7] = "Nao in component membrane (micromolar)"
    legend_constants[8] = "Cao in component membrane (micromolar)"
    legend_constants[9] = "R in component membrane (joule_per_mole_kelvin)"
    legend_constants[10] = "T in component membrane (kelvin)"
    legend_constants[11] = "F in component membrane (coulomb_per_millimole)"
    legend_algebraic[0] = "i_stim in component membrane (picoA_per_picoF)"
    legend_algebraic[46] = "i_CaL in component L_type_calcium_current (picoA_per_picoF)"
    legend_algebraic[48] = "i_pCa in component calcium_pump_current (picoA_per_picoF)"
    legend_algebraic[50] = "i_NaCa in component sodium_calcium_exchange_current (picoA_per_picoF)"
    legend_algebraic[54] = "i_Cab in component calcium_background_current (picoA_per_picoF)"
    legend_algebraic[57] = "i_Na in component fast_sodium_current (picoA_per_picoF)"
    legend_algebraic[58] = "i_Nab in component sodium_background_current (picoA_per_picoF)"
    legend_algebraic[68] = "i_NaK in component sodium_potassium_pump_current (picoA_per_picoF)"
    legend_algebraic[60] = "i_Kto_f in component fast_transient_outward_potassium_current (picoA_per_picoF)"
    legend_algebraic[61] = "i_Kto_s in component slow_transient_outward_potassium_current (picoA_per_picoF)"
    legend_algebraic[62] = "i_K1 in component time_independent_potassium_current (picoA_per_picoF)"
    legend_algebraic[63] = "i_Ks in component slow_delayed_rectifier_potassium_current (picoA_per_picoF)"
    legend_algebraic[64] = "i_Kur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (picoA_per_picoF)"
    legend_algebraic[65] = "i_Kss in component non_inactivating_steady_state_potassium_current (picoA_per_picoF)"
    legend_algebraic[70] = "i_ClCa in component calcium_activated_chloride_current (picoA_per_picoF)"
    legend_algebraic[66] = "i_Kr in component rapid_delayed_rectifier_potassium_current (picoA_per_picoF)"
    legend_constants[12] = "stim_start in component membrane (millisecond)"
    legend_constants[13] = "stim_end in component membrane (millisecond)"
    legend_constants[14] = "stim_period in component membrane (millisecond)"
    legend_constants[15] = "stim_duration in component membrane (millisecond)"
    legend_constants[16] = "stim_amplitude in component membrane (picoA_per_picoF)"
    legend_states[1] = "Cai in component calcium_concentration (micromolar)"
    legend_states[2] = "Cass in component calcium_concentration (micromolar)"
    legend_states[3] = "CaJSR in component calcium_concentration (micromolar)"
    legend_states[4] = "CaNSR in component calcium_concentration (micromolar)"
    legend_algebraic[11] = "Bi in component calcium_concentration (dimensionless)"
    legend_algebraic[24] = "Bss in component calcium_concentration (dimensionless)"
    legend_algebraic[29] = "BJSR in component calcium_concentration (dimensionless)"
    legend_constants[17] = "CMDN_tot in component calcium_concentration (micromolar)"
    legend_constants[18] = "CSQN_tot in component calcium_concentration (micromolar)"
    legend_constants[19] = "Km_CMDN in component calcium_concentration (micromolar)"
    legend_constants[20] = "Km_CSQN in component calcium_concentration (micromolar)"
    legend_algebraic[40] = "J_leak in component calcium_fluxes (micromolar_per_millisecond)"
    legend_algebraic[33] = "J_rel in component calcium_fluxes (micromolar_per_millisecond)"
    legend_algebraic[42] = "J_up in component calcium_fluxes (micromolar_per_millisecond)"
    legend_algebraic[36] = "J_tr in component calcium_fluxes (micromolar_per_millisecond)"
    legend_algebraic[44] = "J_trpn in component calcium_fluxes (micromolar_per_millisecond)"
    legend_algebraic[38] = "J_xfer in component calcium_fluxes (micromolar_per_millisecond)"
    legend_constants[21] = "k_plus_htrpn in component calcium_fluxes (per_micromolar_millisecond)"
    legend_constants[22] = "k_minus_htrpn in component calcium_fluxes (per_millisecond)"
    legend_constants[23] = "k_plus_ltrpn in component calcium_fluxes (per_micromolar_millisecond)"
    legend_constants[24] = "k_minus_ltrpn in component calcium_fluxes (per_millisecond)"
    legend_states[5] = "P_RyR in component calcium_fluxes (dimensionless)"
    legend_constants[25] = "v1 in component calcium_fluxes (per_millisecond)"
    legend_constants[26] = "tau_tr in component calcium_fluxes (millisecond)"
    legend_constants[27] = "v2 in component calcium_fluxes (per_millisecond)"
    legend_constants[28] = "tau_xfer in component calcium_fluxes (millisecond)"
    legend_constants[29] = "v3 in component calcium_fluxes (micromolar_per_millisecond)"
    legend_constants[30] = "Km_up in component calcium_fluxes (micromolar)"
    legend_constants[31] = "LTRPN_tot in component calcium_buffering (micromolar)"
    legend_constants[32] = "HTRPN_tot in component calcium_buffering (micromolar)"
    legend_states[6] = "LTRPN_Ca in component calcium_buffering (micromolar)"
    legend_states[7] = "HTRPN_Ca in component calcium_buffering (micromolar)"
    legend_constants[33] = "i_CaL_max in component L_type_calcium_current (picoA_per_picoF)"
    legend_states[8] = "P_O1 in component ryanodine_receptors (dimensionless)"
    legend_states[9] = "P_O2 in component ryanodine_receptors (dimensionless)"
    legend_algebraic[1] = "P_C1 in component ryanodine_receptors (dimensionless)"
    legend_states[10] = "P_C2 in component ryanodine_receptors (dimensionless)"
    legend_constants[34] = "k_plus_a in component ryanodine_receptors (micromolar4_per_millisecond)"
    legend_constants[35] = "k_minus_a in component ryanodine_receptors (per_millisecond)"
    legend_constants[36] = "k_plus_b in component ryanodine_receptors (micromolar3_per_millisecond)"
    legend_constants[37] = "k_minus_b in component ryanodine_receptors (per_millisecond)"
    legend_constants[38] = "k_plus_c in component ryanodine_receptors (per_millisecond)"
    legend_constants[39] = "k_minus_c in component ryanodine_receptors (per_millisecond)"
    legend_constants[40] = "m in component ryanodine_receptors (dimensionless)"
    legend_constants[41] = "n in component ryanodine_receptors (dimensionless)"
    legend_constants[42] = "E_CaL in component L_type_calcium_current (millivolt)"
    legend_constants[43] = "g_CaL in component L_type_calcium_current (milliS_per_microF)"
    legend_states[11] = "O in component L_type_calcium_current (dimensionless)"
    legend_algebraic[2] = "C1 in component L_type_calcium_current (dimensionless)"
    legend_states[12] = "C2 in component L_type_calcium_current (dimensionless)"
    legend_states[13] = "C3 in component L_type_calcium_current (dimensionless)"
    legend_states[14] = "C4 in component L_type_calcium_current (dimensionless)"
    legend_states[15] = "I1 in component L_type_calcium_current (dimensionless)"
    legend_states[16] = "I2 in component L_type_calcium_current (dimensionless)"
    legend_states[17] = "I3 in component L_type_calcium_current (dimensionless)"
    legend_algebraic[12] = "alpha in component L_type_calcium_current (per_millisecond)"
    legend_algebraic[25] = "beta in component L_type_calcium_current (per_millisecond)"
    legend_algebraic[30] = "gamma in component L_type_calcium_current (per_millisecond)"
    legend_algebraic[34] = "Kpcf in component L_type_calcium_current (per_millisecond)"
    legend_constants[44] = "Kpcb in component L_type_calcium_current (per_millisecond)"
    legend_constants[45] = "Kpc_max in component L_type_calcium_current (per_millisecond)"
    legend_constants[46] = "Kpc_half in component L_type_calcium_current (micromolar)"
    legend_constants[47] = "i_pCa_max in component calcium_pump_current (picoA_per_picoF)"
    legend_constants[48] = "Km_pCa in component calcium_pump_current (micromolar)"
    legend_constants[49] = "k_NaCa in component sodium_calcium_exchange_current (picoA_per_picoF)"
    legend_constants[50] = "K_mNa in component sodium_calcium_exchange_current (micromolar)"
    legend_constants[51] = "K_mCa in component sodium_calcium_exchange_current (micromolar)"
    legend_constants[52] = "k_sat in component sodium_calcium_exchange_current (dimensionless)"
    legend_constants[53] = "eta in component sodium_calcium_exchange_current (dimensionless)"
    legend_states[18] = "Nai in component sodium_concentration (micromolar)"
    legend_constants[54] = "g_Cab in component calcium_background_current (milliS_per_microF)"
    legend_algebraic[52] = "E_CaN in component calcium_background_current (millivolt)"
    legend_algebraic[56] = "E_Na in component fast_sodium_current (millivolt)"
    legend_constants[55] = "g_Na in component fast_sodium_current (milliS_per_microF)"
    legend_states[19] = "O_Na in component fast_sodium_current (dimensionless)"
    legend_states[20] = "C_Na1 in component fast_sodium_current (dimensionless)"
    legend_states[21] = "C_Na2 in component fast_sodium_current (dimensionless)"
    legend_algebraic[3] = "C_Na3 in component fast_sodium_current (dimensionless)"
    legend_states[22] = "I1_Na in component fast_sodium_current (dimensionless)"
    legend_states[23] = "I2_Na in component fast_sodium_current (dimensionless)"
    legend_states[24] = "IF_Na in component fast_sodium_current (dimensionless)"
    legend_states[25] = "IC_Na2 in component fast_sodium_current (dimensionless)"
    legend_states[26] = "IC_Na3 in component fast_sodium_current (dimensionless)"
    legend_algebraic[13] = "alpha_Na11 in component fast_sodium_current (per_millisecond)"
    legend_algebraic[35] = "beta_Na11 in component fast_sodium_current (per_millisecond)"
    legend_algebraic[26] = "alpha_Na12 in component fast_sodium_current (per_millisecond)"
    legend_algebraic[37] = "beta_Na12 in component fast_sodium_current (per_millisecond)"
    legend_algebraic[31] = "alpha_Na13 in component fast_sodium_current (per_millisecond)"
    legend_algebraic[39] = "beta_Na13 in component fast_sodium_current (per_millisecond)"
    legend_algebraic[41] = "alpha_Na3 in component fast_sodium_current (per_millisecond)"
    legend_algebraic[43] = "beta_Na3 in component fast_sodium_current (per_millisecond)"
    legend_algebraic[45] = "alpha_Na2 in component fast_sodium_current (per_millisecond)"
    legend_algebraic[47] = "beta_Na2 in component fast_sodium_current (per_millisecond)"
    legend_algebraic[49] = "alpha_Na4 in component fast_sodium_current (per_millisecond)"
    legend_algebraic[51] = "beta_Na4 in component fast_sodium_current (per_millisecond)"
    legend_algebraic[53] = "alpha_Na5 in component fast_sodium_current (per_millisecond)"
    legend_algebraic[55] = "beta_Na5 in component fast_sodium_current (per_millisecond)"
    legend_states[27] = "Ki in component potassium_concentration (micromolar)"
    legend_constants[56] = "g_Nab in component sodium_background_current (milliS_per_microF)"
    legend_algebraic[59] = "E_K in component fast_transient_outward_potassium_current (millivolt)"
    legend_constants[57] = "g_Kto_f in component fast_transient_outward_potassium_current (milliS_per_microF)"
    legend_states[28] = "ato_f in component fast_transient_outward_potassium_current (dimensionless)"
    legend_states[29] = "ito_f in component fast_transient_outward_potassium_current (dimensionless)"
    legend_algebraic[4] = "alpha_a in component fast_transient_outward_potassium_current (per_millisecond)"
    legend_algebraic[14] = "beta_a in component fast_transient_outward_potassium_current (per_millisecond)"
    legend_algebraic[5] = "alpha_i in component fast_transient_outward_potassium_current (per_millisecond)"
    legend_algebraic[15] = "beta_i in component fast_transient_outward_potassium_current (per_millisecond)"
    legend_algebraic[6] = "ass in component slow_transient_outward_potassium_current (dimensionless)"
    legend_algebraic[7] = "iss in component slow_transient_outward_potassium_current (dimensionless)"
    legend_constants[58] = "g_Kto_s in component slow_transient_outward_potassium_current (milliS_per_microF)"
    legend_states[30] = "ato_s in component slow_transient_outward_potassium_current (dimensionless)"
    legend_states[31] = "ito_s in component slow_transient_outward_potassium_current (dimensionless)"
    legend_algebraic[16] = "tau_ta_s in component slow_transient_outward_potassium_current (millisecond)"
    legend_algebraic[17] = "tau_ti_s in component slow_transient_outward_potassium_current (millisecond)"
    legend_constants[59] = "g_Ks in component slow_delayed_rectifier_potassium_current (milliS_per_microF)"
    legend_states[32] = "nKs in component slow_delayed_rectifier_potassium_current (dimensionless)"
    legend_algebraic[8] = "alpha_n in component slow_delayed_rectifier_potassium_current (per_millisecond)"
    legend_algebraic[18] = "beta_n in component slow_delayed_rectifier_potassium_current (per_millisecond)"
    legend_constants[60] = "g_Kur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (milliS_per_microF)"
    legend_states[33] = "aur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)"
    legend_states[34] = "iur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)"
    legend_algebraic[19] = "tau_aur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (millisecond)"
    legend_algebraic[20] = "tau_iur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (millisecond)"
    legend_constants[61] = "g_Kss in component non_inactivating_steady_state_potassium_current (milliS_per_microF)"
    legend_states[35] = "aKss in component non_inactivating_steady_state_potassium_current (dimensionless)"
    legend_states[36] = "iKss in component non_inactivating_steady_state_potassium_current (dimensionless)"
    legend_algebraic[21] = "tau_Kss in component non_inactivating_steady_state_potassium_current (millisecond)"
    legend_constants[62] = "g_Kr in component rapid_delayed_rectifier_potassium_current (milliS_per_microF)"
    legend_states[37] = "O_K in component rapid_delayed_rectifier_potassium_current (dimensionless)"
    legend_states[38] = "C_K1 in component rapid_delayed_rectifier_potassium_current (dimensionless)"
    legend_states[39] = "C_K2 in component rapid_delayed_rectifier_potassium_current (dimensionless)"
    legend_algebraic[9] = "C_K0 in component rapid_delayed_rectifier_potassium_current (dimensionless)"
    legend_states[40] = "I_K in component rapid_delayed_rectifier_potassium_current (dimensionless)"
    legend_algebraic[22] = "alpha_a0 in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
    legend_algebraic[27] = "beta_a0 in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
    legend_constants[63] = "kb in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
    legend_constants[64] = "kf in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
    legend_algebraic[10] = "alpha_a1 in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
    legend_algebraic[23] = "beta_a1 in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
    legend_algebraic[28] = "alpha_i in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
    legend_algebraic[32] = "beta_i in component rapid_delayed_rectifier_potassium_current (per_millisecond)"
    legend_constants[65] = "i_NaK_max in component sodium_potassium_pump_current (picoA_per_picoF)"
    legend_constants[66] = "Km_Nai in component sodium_potassium_pump_current (micromolar)"
    legend_constants[67] = "Km_Ko in component sodium_potassium_pump_current (micromolar)"
    legend_algebraic[67] = "f_NaK in component sodium_potassium_pump_current (dimensionless)"
    legend_constants[71] = "sigma in component sodium_potassium_pump_current (dimensionless)"
    legend_constants[68] = "g_ClCa in component calcium_activated_chloride_current (milliS_per_microF)"
    legend_algebraic[69] = "O_ClCa in component calcium_activated_chloride_current (dimensionless)"
    legend_constants[69] = "E_Cl in component calcium_activated_chloride_current (millivolt)"
    legend_constants[70] = "Km_Cl in component calcium_activated_chloride_current (micromolar)"
    legend_rates[0] = "d/dt V in component membrane (millivolt)"
    legend_rates[1] = "d/dt Cai in component calcium_concentration (micromolar)"
    legend_rates[2] = "d/dt Cass in component calcium_concentration (micromolar)"
    legend_rates[3] = "d/dt CaJSR in component calcium_concentration (micromolar)"
    legend_rates[4] = "d/dt CaNSR in component calcium_concentration (micromolar)"
    legend_rates[5] = "d/dt P_RyR in component calcium_fluxes (dimensionless)"
    legend_rates[6] = "d/dt LTRPN_Ca in component calcium_buffering (micromolar)"
    legend_rates[7] = "d/dt HTRPN_Ca in component calcium_buffering (micromolar)"
    legend_rates[8] = "d/dt P_O1 in component ryanodine_receptors (dimensionless)"
    legend_rates[9] = "d/dt P_O2 in component ryanodine_receptors (dimensionless)"
    legend_rates[10] = "d/dt P_C2 in component ryanodine_receptors (dimensionless)"
    legend_rates[11] = "d/dt O in component L_type_calcium_current (dimensionless)"
    legend_rates[12] = "d/dt C2 in component L_type_calcium_current (dimensionless)"
    legend_rates[13] = "d/dt C3 in component L_type_calcium_current (dimensionless)"
    legend_rates[14] = "d/dt C4 in component L_type_calcium_current (dimensionless)"
    legend_rates[15] = "d/dt I1 in component L_type_calcium_current (dimensionless)"
    legend_rates[16] = "d/dt I2 in component L_type_calcium_current (dimensionless)"
    legend_rates[17] = "d/dt I3 in component L_type_calcium_current (dimensionless)"
    legend_rates[18] = "d/dt Nai in component sodium_concentration (micromolar)"
    legend_rates[21] = "d/dt C_Na2 in component fast_sodium_current (dimensionless)"
    legend_rates[20] = "d/dt C_Na1 in component fast_sodium_current (dimensionless)"
    legend_rates[19] = "d/dt O_Na in component fast_sodium_current (dimensionless)"
    legend_rates[24] = "d/dt IF_Na in component fast_sodium_current (dimensionless)"
    legend_rates[22] = "d/dt I1_Na in component fast_sodium_current (dimensionless)"
    legend_rates[23] = "d/dt I2_Na in component fast_sodium_current (dimensionless)"
    legend_rates[25] = "d/dt IC_Na2 in component fast_sodium_current (dimensionless)"
    legend_rates[26] = "d/dt IC_Na3 in component fast_sodium_current (dimensionless)"
    legend_rates[27] = "d/dt Ki in component potassium_concentration (micromolar)"
    legend_rates[28] = "d/dt ato_f in component fast_transient_outward_potassium_current (dimensionless)"
    legend_rates[29] = "d/dt ito_f in component fast_transient_outward_potassium_current (dimensionless)"
    legend_rates[30] = "d/dt ato_s in component slow_transient_outward_potassium_current (dimensionless)"
    legend_rates[31] = "d/dt ito_s in component slow_transient_outward_potassium_current (dimensionless)"
    legend_rates[32] = "d/dt nKs in component slow_delayed_rectifier_potassium_current (dimensionless)"
    legend_rates[33] = "d/dt aur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)"
    legend_rates[34] = "d/dt iur in component ultra_rapidly_activating_delayed_rectifier_potassium_current (dimensionless)"
    legend_rates[35] = "d/dt aKss in component non_inactivating_steady_state_potassium_current (dimensionless)"
    legend_rates[36] = "d/dt iKss in component non_inactivating_steady_state_potassium_current (dimensionless)"
    legend_rates[39] = "d/dt C_K2 in component rapid_delayed_rectifier_potassium_current (dimensionless)"
    legend_rates[38] = "d/dt C_K1 in component rapid_delayed_rectifier_potassium_current (dimensionless)"
    legend_rates[37] = "d/dt O_K in component rapid_delayed_rectifier_potassium_current (dimensionless)"
    legend_rates[40] = "d/dt I_K in component rapid_delayed_rectifier_potassium_current (dimensionless)"
    return (legend_states, legend_algebraic, legend_voi, legend_constants)

def initConsts():
    constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
    states[0] = -82.4202
    constants[0] = 1
    constants[1] = 25.84e-6
    constants[2] = 0.12e-6
    constants[3] = 2.098e-6
    constants[4] = 1.485e-9
    constants[5] = 1.534e-4
    constants[6] = 5400
    constants[7] = 140000
    constants[8] = 1800
    constants[9] = 8.314
    constants[10] = 298
    constants[11] = 96.5
    constants[12] = 20
    constants[13] = 100000
    constants[14] = 71.43
    constants[15] = 0.5
    constants[16] = -80
    states[1] = 0.115001
    states[2] = 0.115001
    states[3] = 1299.5
    states[4] = 1299.5
    constants[17] = 50
    constants[18] = 15000
    constants[19] = 0.238
    constants[20] = 800
    constants[21] = 0.00237
    constants[22] = 3.2e-5
    constants[23] = 0.0327
    constants[24] = 0.0196
    states[5] = 0
    constants[25] = 4.5
    constants[26] = 20
    constants[27] = 1.74e-5
    constants[28] = 8
    constants[29] = 0.45
    constants[30] = 0.5
    constants[31] = 70
    constants[32] = 140
    states[6] = 11.2684
    states[7] = 125.29
    constants[33] = 7
    states[8] = 0.149102e-4
    states[9] = 0.951726e-10
    states[10] = 0.16774e-3
    constants[34] = 0.006075
    constants[35] = 0.07125
    constants[36] = 0.00405
    constants[37] = 0.965
    constants[38] = 0.009
    constants[39] = 0.0008
    constants[40] = 3
    constants[41] = 4
    constants[42] = 63
    constants[43] = 0.1729
    states[11] = 0.930308e-18
    states[12] = 0.124216e-3
    states[13] = 0.578679e-8
    states[14] = 0.119816e-12
    states[15] = 0.497923e-18
    states[16] = 0.345847e-13
    states[17] = 0.185106e-13
    constants[44] = 0.0005
    constants[45] = 0.23324
    constants[46] = 20
    constants[47] = 1
    constants[48] = 0.5
    constants[49] = 292.8
    constants[50] = 87500
    constants[51] = 1380
    constants[52] = 0.1
    constants[53] = 0.35
    states[18] = 14237.1
    constants[54] = 0.000367
    constants[55] = 13
    states[19] = 0.713483e-6
    states[20] = 0.279132e-3
    states[21] = 0.020752
    states[22] = 0.673345e-6
    states[23] = 0.155787e-8
    states[24] = 0.153176e-3
    states[25] = 0.0113879
    states[26] = 0.34278
    states[27] = 143720
    constants[56] = 0.0026
    constants[57] = 0.0798
    states[28] = 0.265563e-2
    states[29] = 0.999977
    constants[58] = 0.0629
    states[30] = 0.417069e-3
    states[31] = 0.998543
    constants[59] = 0.00575
    states[32] = 0.262753e-3
    constants[60] = 0.0975
    states[33] = 0.417069e-3
    states[34] = 0.998543
    constants[61] = 0.0324
    states[35] = 0.417069e-3
    states[36] = 1
    constants[62] = 0.078
    states[37] = 0.175298e-3
    states[38] = 0.992513e-3
    states[39] = 0.641229e-3
    states[40] = 0.319129e-4
    constants[63] = 0.036778
    constants[64] = 0.023761
    constants[65] = 0.88
    constants[66] = 21000
    constants[67] = 1500
    constants[68] = 10
    constants[69] = -40
    constants[70] = 10
    constants[71] = (1.00000/7.00000)*(exp(constants[7]/67300.0)-1.00000)
    constants[72] = 0.00000
    return (states, constants)

def computeRates(voi, states, constants):
    rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
    rates[36] = constants[72]
    rates[6] = constants[23]*states[1]*(constants[31]-states[6])-constants[24]*states[6]
    rates[7] = constants[21]*states[1]*(constants[32]-states[7])-constants[22]*states[7]
    rates[9] = constants[36]*(power(states[2], constants[40]))*states[8]-constants[37]*states[9]
    rates[10] = constants[38]*states[8]-constants[39]*states[10]
    algebraic[1] = 1.00000-(states[10]+states[8]+states[9])
    rates[8] = (constants[34]*(power(states[2], constants[41]))*algebraic[1]+constants[37]*states[9]+constants[39]*states[10])-(constants[35]*states[8]+constants[36]*(power(states[2], constants[40]))*states[8]+constants[38]*states[8])
    algebraic[4] = 0.180640*exp(0.0357700*(states[0]+30.0000))
    algebraic[14] = 0.395600*exp(-0.0623700*(states[0]+30.0000))
    rates[28] = algebraic[4]*(1.00000-states[28])-algebraic[14]*states[28]
    algebraic[5] = (0.000152000*exp(-(states[0]+13.5000)/7.00000))/(0.00670830*exp(-(states[0]+33.5000)/7.00000)+1.00000)
    algebraic[15] = (0.000950000*exp((states[0]+33.5000)/7.00000))/(0.0513350*exp((states[0]+33.5000)/7.00000)+1.00000)
    rates[29] = algebraic[5]*(1.00000-states[29])-algebraic[15]*states[29]
    algebraic[6] = 1.00000/(1.00000+exp(-(states[0]+22.5000)/7.70000))
    algebraic[16] = 0.493000*exp(-0.0629000*states[0])+2.05800
    rates[30] = (algebraic[6]-states[30])/algebraic[16]
    algebraic[7] = 1.00000/(1.00000+exp((states[0]+45.2000)/5.70000))
    algebraic[17] = 270.000+1050.00/(1.00000+exp((states[0]+45.2000)/5.70000))
    rates[31] = (algebraic[7]-states[31])/algebraic[17]
    algebraic[8] = (4.81333e-06*(states[0]+26.5000))/(1.00000-exp(-0.128000*(states[0]+26.5000)))
    algebraic[18] = 9.53333e-05*exp(-0.0380000*(states[0]+26.5000))
    rates[32] = algebraic[8]*(1.00000-states[32])-algebraic[18]*states[32]
    algebraic[19] = 0.493000*exp(-0.0629000*states[0])+2.05800
    rates[33] = (algebraic[6]-states[33])/algebraic[19]
    algebraic[20] = 1200.00-170.000/(1.00000+exp((states[0]+45.2000)/5.70000))
    rates[34] = (algebraic[7]-states[34])/algebraic[20]
    algebraic[21] = 39.3000*exp(-0.0862000*states[0])+13.1700
    rates[35] = (algebraic[6]-states[35])/algebraic[21]
    algebraic[10] = 0.0137330*exp(0.0381980*states[0])
    algebraic[23] = 6.89000e-05*exp(-0.0417800*states[0])
    rates[39] = (constants[64]*states[38]+algebraic[23]*states[37])-(constants[63]*states[39]+algebraic[10]*states[39])
    algebraic[2] = 1.00000-(states[11]+states[12]+states[13]+states[14]+states[15]+states[16]+states[17])
    algebraic[12] = (0.400000*exp((states[0]+12.0000)/10.0000)*((1.00000+0.700000*exp(-(power(states[0]+40.0000, 2.00000))/10.0000))-0.750000*exp(-(power(states[0]+20.0000, 2.00000))/400.000)))/(1.00000+0.120000*exp((states[0]+12.0000)/10.0000))
    algebraic[25] = 0.0500000*exp(-(states[0]+12.0000)/13.0000)
    rates[12] = (4.00000*algebraic[12]*algebraic[2]+2.00000*algebraic[25]*states[13])-(algebraic[25]*states[12]+3.00000*algebraic[12]*states[12])
    rates[13] = (3.00000*algebraic[12]*states[12]+3.00000*algebraic[25]*states[14])-(2.00000*algebraic[25]*states[13]+2.00000*algebraic[12]*states[13])
    algebraic[9] = 1.00000-(states[38]+states[39]+states[37]+states[40])
    algebraic[22] = 0.0223480*exp(0.0117600*states[0])
    algebraic[27] = 0.0470020*exp(-0.0631000*states[0])
    rates[38] = (algebraic[22]*algebraic[9]+constants[63]*states[39])-(algebraic[27]*states[38]+constants[64]*states[38])
    algebraic[28] = 0.0908210*exp(0.0233910*(states[0]+5.00000))
    algebraic[32] = 0.00649700*exp(-0.0326800*(states[0]+5.00000))
    rates[37] = (algebraic[10]*states[39]+algebraic[32]*states[40])-(algebraic[23]*states[37]+algebraic[28]*states[37])
    rates[40] = algebraic[28]*states[37]-algebraic[32]*states[40]
    algebraic[30] = (constants[45]*states[2])/(constants[46]+states[2])
    algebraic[34] = 13.0000*(1.00000-exp(-(power(states[0]+14.5000, 2.00000))/100.000))
    rates[11] = (algebraic[12]*states[14]+constants[44]*states[15]+0.00100000*(algebraic[12]*states[16]-algebraic[34]*states[11]))-(4.00000*algebraic[25]*states[11]+algebraic[30]*states[11])
    rates[14] = (2.00000*algebraic[12]*states[13]+4.00000*algebraic[25]*states[11]+0.0100000*(4.00000*constants[44]*algebraic[25]*states[15]-algebraic[12]*algebraic[30]*states[14])+0.00200000*(4.00000*algebraic[25]*states[16]-algebraic[34]*states[14])+4.00000*algebraic[25]*constants[44]*states[17])-(3.00000*algebraic[25]*states[14]+algebraic[12]*states[14]+1.00000*algebraic[30]*algebraic[34]*states[14])
    rates[15] = (algebraic[30]*states[11]+0.00100000*(algebraic[12]*states[17]-algebraic[34]*states[15])+0.0100000*(algebraic[12]*algebraic[30]*states[14]-4.00000*algebraic[25]*algebraic[34]*states[15]))-constants[44]*states[15]
    rates[16] = (0.00100000*(algebraic[34]*states[11]-algebraic[12]*states[16])+constants[44]*states[17]+0.00200000*(algebraic[34]*states[14]-4.00000*algebraic[25]*states[16]))-algebraic[30]*states[16]
    rates[17] = (0.00100000*(algebraic[34]*states[15]-algebraic[12]*states[17])+algebraic[30]*states[16]+1.00000*algebraic[30]*algebraic[34]*states[14])-(4.00000*algebraic[25]*constants[44]*states[17]+constants[44]*states[17])
    algebraic[29] = power(1.00000+(constants[18]*constants[20])/(power(constants[20]+states[3], 2.00000)), -1.00000)
    algebraic[33] = constants[25]*(states[8]+states[9])*(states[3]-states[2])*states[5]
    algebraic[36] = (states[4]-states[3])/constants[26]
    rates[3] = algebraic[29]*(algebraic[36]-algebraic[33])
    algebraic[40] = constants[27]*(states[4]-states[1])
    algebraic[42] = (constants[29]*(power(states[1], 2.00000)))/(power(constants[30], 2.00000)+power(states[1], 2.00000))
    rates[4] = ((algebraic[42]-algebraic[40])*constants[1])/constants[3]-(algebraic[36]*constants[2])/constants[3]
    algebraic[3] = 1.00000-(states[19]+states[20]+states[21]+states[24]+states[22]+states[23]+states[25]+states[26])
    algebraic[13] = 3.80200/(0.102700*exp(-(states[0]+2.50000)/17.0000)+0.200000*exp(-(states[0]+2.50000)/150.000))
    algebraic[35] = 0.191700*exp(-(states[0]+2.50000)/20.3000)
    algebraic[26] = 3.80200/(0.102700*exp(-(states[0]+2.50000)/15.0000)+0.230000*exp(-(states[0]+2.50000)/150.000))
    algebraic[37] = 0.200000*exp(-(states[0]-2.50000)/20.3000)
    algebraic[41] = 7.00000e-07*exp(-(states[0]+7.00000)/7.70000)
    algebraic[43] = 0.00840000+2.00000e-05*(states[0]+7.00000)
    rates[21] = (algebraic[13]*algebraic[3]+algebraic[37]*states[20]+algebraic[41]*states[25])-(algebraic[35]*states[21]+algebraic[26]*states[21]+algebraic[43]*states[21])
    algebraic[31] = 3.80200/(0.102700*exp(-(states[0]+2.50000)/12.0000)+0.250000*exp(-(states[0]+2.50000)/150.000))
    algebraic[39] = 0.220000*exp(-(states[0]-7.50000)/20.3000)
    rates[20] = (algebraic[26]*states[21]+algebraic[39]*states[19]+algebraic[41]*states[24])-(algebraic[37]*states[20]+algebraic[31]*states[20]+algebraic[43]*states[20])
    rates[25] = (algebraic[13]*states[26]+algebraic[37]*states[24]+algebraic[43]*states[21])-(algebraic[35]*states[25]+algebraic[26]*states[25]+algebraic[41]*states[25])
    rates[26] = (algebraic[35]*states[25]+algebraic[43]*algebraic[3])-(algebraic[13]*states[26]+algebraic[41]*states[26])
    algebraic[46] = constants[43]*states[11]*(states[0]-constants[42])
    algebraic[24] = power(1.00000+(constants[17]*constants[19])/(power(constants[19]+states[2], 2.00000)), -1.00000)
    algebraic[38] = (states[2]-states[1])/constants[28]
    rates[2] = algebraic[24]*((algebraic[33]*constants[2])/constants[4]-((algebraic[38]*constants[1])/constants[4]+(algebraic[46]*constants[5]*constants[0])/(2.00000*constants[4]*constants[11])))
    rates[5] = -0.0400000*states[5]-((0.100000*algebraic[46])/constants[33])*exp(-(power(states[0]-5.00000, 2.00000))/648.000)
    algebraic[45] = 1.00000/(0.188495*exp(-(states[0]+7.00000)/16.6000)+0.393956)
    algebraic[47] = (algebraic[31]*algebraic[45]*algebraic[41])/(algebraic[39]*algebraic[43])
    rates[19] = (algebraic[31]*states[20]+algebraic[47]*states[24])-(algebraic[39]*states[19]+algebraic[45]*states[19])
    algebraic[49] = algebraic[45]/1000.00
    algebraic[51] = algebraic[41]
    rates[24] = (algebraic[45]*states[19]+algebraic[43]*states[20]+algebraic[51]*states[22]+algebraic[26]*states[25])-(algebraic[47]*states[24]+algebraic[41]*states[24]+algebraic[49]*states[24]+algebraic[37]*states[24])
    algebraic[48] = (constants[47]*(power(states[1], 2.00000)))/(power(constants[48], 2.00000)+power(states[1], 2.00000))
    algebraic[50] = ((((((constants[49]*1.00000)/(power(constants[50], 3.00000)+power(constants[7], 3.00000)))*1.00000)/(constants[51]+constants[8]))*1.00000)/(1.00000+constants[52]*exp(((constants[53]-1.00000)*states[0]*constants[11])/(constants[9]*constants[10]))))*(exp((constants[53]*states[0]*constants[11])/(constants[9]*constants[10]))*(power(states[18], 3.00000))*constants[8]-exp(((constants[53]-1.00000)*states[0]*constants[11])/(constants[9]*constants[10]))*(power(constants[7], 3.00000))*states[1])
    algebraic[52] = ((constants[9]*constants[10])/(2.00000*constants[11]))*log(constants[8]/states[1])
    algebraic[54] = constants[54]*(states[0]-algebraic[52])
    algebraic[11] = power(1.00000+(constants[17]*constants[19])/(power(constants[19]+states[1], 2.00000)), -1.00000)
    algebraic[44] = (constants[21]*states[1]*(constants[32]-states[7])+constants[23]*states[1]*(constants[31]-states[6]))-(constants[22]*states[7]+constants[24]*states[6])
    rates[1] = algebraic[11]*((algebraic[40]+algebraic[38])-(algebraic[42]+algebraic[44]+(((algebraic[54]+algebraic[48])-2.00000*algebraic[50])*constants[5]*constants[0])/(2.00000*constants[1]*constants[11])))
    algebraic[53] = algebraic[45]/95000.0
    algebraic[55] = algebraic[41]/50.0000
    rates[22] = (algebraic[49]*states[24]+algebraic[55]*states[23])-(algebraic[51]*states[22]+algebraic[53]*states[22])
    rates[23] = algebraic[53]*states[22]-algebraic[55]*states[23]
    algebraic[56] = ((constants[9]*constants[10])/constants[11])*log((0.900000*constants[7]+0.100000*constants[6])/(0.900000*states[18]+0.100000*states[27]))
    algebraic[57] = constants[55]*states[19]*(states[0]-algebraic[56])
    algebraic[58] = constants[56]*(states[0]-algebraic[56])
    algebraic[67] = 1.00000/(1.00000+0.124500*exp((-0.100000*states[0]*constants[11])/(constants[9]*constants[10]))+0.0365000*constants[71]*exp((-states[0]*constants[11])/(constants[9]*constants[10])))
    algebraic[68] = (((constants[65]*algebraic[67]*1.00000)/(1.00000+power(constants[66]/states[18], 1.50000)))*constants[6])/(constants[6]+constants[67])
    rates[18] = (-(algebraic[57]+algebraic[58]+3.00000*algebraic[68]+3.00000*algebraic[50])*constants[5]*constants[0])/(constants[1]*constants[11])
    algebraic[59] = ((constants[9]*constants[10])/constants[11])*log(constants[6]/states[27])
    algebraic[60] = constants[57]*(power(states[28], 3.00000))*states[29]*(states[0]-algebraic[59])
    algebraic[61] = constants[58]*states[30]*states[31]*(states[0]-algebraic[59])
    algebraic[62] = (((0.293800*constants[6])/(constants[6]+210.000))*(states[0]-algebraic[59]))/(1.00000+exp(0.0896000*(states[0]-algebraic[59])))
    algebraic[63] = constants[59]*(power(states[32], 2.00000))*(states[0]-algebraic[59])
    algebraic[64] = constants[60]*states[33]*states[34]*(states[0]-algebraic[59])
    algebraic[65] = constants[61]*states[35]*states[36]*(states[0]-algebraic[59])
    algebraic[66] = constants[62]*states[37]*(states[0]-((constants[9]*constants[10])/constants[11])*log((0.980000*constants[6]+0.0200000*constants[7])/(0.980000*states[27]+0.0200000*states[18])))
    rates[27] = (-((algebraic[60]+algebraic[61]+algebraic[62]+algebraic[63]+algebraic[65]+algebraic[64]+algebraic[66])-2.00000*algebraic[68])*constants[5]*constants[0])/(constants[1]*constants[11])
    algebraic[0] = custom_piecewise([greater_equal(voi , constants[12]) & less_equal(voi , constants[13]) & less_equal((voi-constants[12])-floor((voi-constants[12])/constants[14])*constants[14] , constants[15]), constants[16] , True, 0.00000])
    algebraic[69] = 0.200000/(1.00000+exp(-(states[0]-46.7000)/7.80000))
    algebraic[70] = ((constants[68]*algebraic[69]*states[1])/(states[1]+constants[70]))*(states[0]-constants[69])
    rates[0] = -(algebraic[46]+algebraic[48]+algebraic[50]+algebraic[54]+algebraic[57]+algebraic[58]+algebraic[68]+algebraic[60]+algebraic[61]+algebraic[62]+algebraic[63]+algebraic[64]+algebraic[65]+algebraic[66]+algebraic[70]+algebraic[0])
    return(rates)

def computeAlgebraic(constants, states, voi):
    algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
    states = array(states)
    voi = array(voi)
    algebraic[1] = 1.00000-(states[10]+states[8]+states[9])
    algebraic[4] = 0.180640*exp(0.0357700*(states[0]+30.0000))
    algebraic[14] = 0.395600*exp(-0.0623700*(states[0]+30.0000))
    algebraic[5] = (0.000152000*exp(-(states[0]+13.5000)/7.00000))/(0.00670830*exp(-(states[0]+33.5000)/7.00000)+1.00000)
    algebraic[15] = (0.000950000*exp((states[0]+33.5000)/7.00000))/(0.0513350*exp((states[0]+33.5000)/7.00000)+1.00000)
    algebraic[6] = 1.00000/(1.00000+exp(-(states[0]+22.5000)/7.70000))
    algebraic[16] = 0.493000*exp(-0.0629000*states[0])+2.05800
    algebraic[7] = 1.00000/(1.00000+exp((states[0]+45.2000)/5.70000))
    algebraic[17] = 270.000+1050.00/(1.00000+exp((states[0]+45.2000)/5.70000))
    algebraic[8] = (4.81333e-06*(states[0]+26.5000))/(1.00000-exp(-0.128000*(states[0]+26.5000)))
    algebraic[18] = 9.53333e-05*exp(-0.0380000*(states[0]+26.5000))
    algebraic[19] = 0.493000*exp(-0.0629000*states[0])+2.05800
    algebraic[20] = 1200.00-170.000/(1.00000+exp((states[0]+45.2000)/5.70000))
    algebraic[21] = 39.3000*exp(-0.0862000*states[0])+13.1700
    algebraic[10] = 0.0137330*exp(0.0381980*states[0])
    algebraic[23] = 6.89000e-05*exp(-0.0417800*states[0])
    algebraic[2] = 1.00000-(states[11]+states[12]+states[13]+states[14]+states[15]+states[16]+states[17])
    algebraic[12] = (0.400000*exp((states[0]+12.0000)/10.0000)*((1.00000+0.700000*exp(-(power(states[0]+40.0000, 2.00000))/10.0000))-0.750000*exp(-(power(states[0]+20.0000, 2.00000))/400.000)))/(1.00000+0.120000*exp((states[0]+12.0000)/10.0000))
    algebraic[25] = 0.0500000*exp(-(states[0]+12.0000)/13.0000)
    algebraic[9] = 1.00000-(states[38]+states[39]+states[37]+states[40])
    algebraic[22] = 0.0223480*exp(0.0117600*states[0])
    algebraic[27] = 0.0470020*exp(-0.0631000*states[0])
    algebraic[28] = 0.0908210*exp(0.0233910*(states[0]+5.00000))
    algebraic[32] = 0.00649700*exp(-0.0326800*(states[0]+5.00000))
    algebraic[30] = (constants[45]*states[2])/(constants[46]+states[2])
    algebraic[34] = 13.0000*(1.00000-exp(-(power(states[0]+14.5000, 2.00000))/100.000))
    algebraic[29] = power(1.00000+(constants[18]*constants[20])/(power(constants[20]+states[3], 2.00000)), -1.00000)
    algebraic[33] = constants[25]*(states[8]+states[9])*(states[3]-states[2])*states[5]
    algebraic[36] = (states[4]-states[3])/constants[26]
    algebraic[40] = constants[27]*(states[4]-states[1])
    algebraic[42] = (constants[29]*(power(states[1], 2.00000)))/(power(constants[30], 2.00000)+power(states[1], 2.00000))
    algebraic[3] = 1.00000-(states[19]+states[20]+states[21]+states[24]+states[22]+states[23]+states[25]+states[26])
    algebraic[13] = 3.80200/(0.102700*exp(-(states[0]+2.50000)/17.0000)+0.200000*exp(-(states[0]+2.50000)/150.000))
    algebraic[35] = 0.191700*exp(-(states[0]+2.50000)/20.3000)
    algebraic[26] = 3.80200/(0.102700*exp(-(states[0]+2.50000)/15.0000)+0.230000*exp(-(states[0]+2.50000)/150.000))
    algebraic[37] = 0.200000*exp(-(states[0]-2.50000)/20.3000)
    algebraic[41] = 7.00000e-07*exp(-(states[0]+7.00000)/7.70000)
    algebraic[43] = 0.00840000+2.00000e-05*(states[0]+7.00000)
    algebraic[31] = 3.80200/(0.102700*exp(-(states[0]+2.50000)/12.0000)+0.250000*exp(-(states[0]+2.50000)/150.000))
    algebraic[39] = 0.220000*exp(-(states[0]-7.50000)/20.3000)
    algebraic[46] = constants[43]*states[11]*(states[0]-constants[42])
    algebraic[24] = power(1.00000+(constants[17]*constants[19])/(power(constants[19]+states[2], 2.00000)), -1.00000)
    algebraic[38] = (states[2]-states[1])/constants[28]
    algebraic[45] = 1.00000/(0.188495*exp(-(states[0]+7.00000)/16.6000)+0.393956)
    algebraic[47] = (algebraic[31]*algebraic[45]*algebraic[41])/(algebraic[39]*algebraic[43])
    algebraic[49] = algebraic[45]/1000.00
    algebraic[51] = algebraic[41]
    algebraic[48] = (constants[47]*(power(states[1], 2.00000)))/(power(constants[48], 2.00000)+power(states[1], 2.00000))
    algebraic[50] = ((((((constants[49]*1.00000)/(power(constants[50], 3.00000)+power(constants[7], 3.00000)))*1.00000)/(constants[51]+constants[8]))*1.00000)/(1.00000+constants[52]*exp(((constants[53]-1.00000)*states[0]*constants[11])/(constants[9]*constants[10]))))*(exp((constants[53]*states[0]*constants[11])/(constants[9]*constants[10]))*(power(states[18], 3.00000))*constants[8]-exp(((constants[53]-1.00000)*states[0]*constants[11])/(constants[9]*constants[10]))*(power(constants[7], 3.00000))*states[1])
    algebraic[52] = ((constants[9]*constants[10])/(2.00000*constants[11]))*log(constants[8]/states[1])
    algebraic[54] = constants[54]*(states[0]-algebraic[52])
    algebraic[11] = power(1.00000+(constants[17]*constants[19])/(power(constants[19]+states[1], 2.00000)), -1.00000)
    algebraic[44] = (constants[21]*states[1]*(constants[32]-states[7])+constants[23]*states[1]*(constants[31]-states[6]))-(constants[22]*states[7]+constants[24]*states[6])
    algebraic[53] = algebraic[45]/95000.0
    algebraic[55] = algebraic[41]/50.0000
    algebraic[56] = ((constants[9]*constants[10])/constants[11])*log((0.900000*constants[7]+0.100000*constants[6])/(0.900000*states[18]+0.100000*states[27]))
    algebraic[57] = constants[55]*states[19]*(states[0]-algebraic[56])
    algebraic[58] = constants[56]*(states[0]-algebraic[56])
    algebraic[67] = 1.00000/(1.00000+0.124500*exp((-0.100000*states[0]*constants[11])/(constants[9]*constants[10]))+0.0365000*constants[71]*exp((-states[0]*constants[11])/(constants[9]*constants[10])))
    algebraic[68] = (((constants[65]*algebraic[67]*1.00000)/(1.00000+power(constants[66]/states[18], 1.50000)))*constants[6])/(constants[6]+constants[67])
    algebraic[59] = ((constants[9]*constants[10])/constants[11])*log(constants[6]/states[27])
    algebraic[60] = constants[57]*(power(states[28], 3.00000))*states[29]*(states[0]-algebraic[59])
    algebraic[61] = constants[58]*states[30]*states[31]*(states[0]-algebraic[59])
    algebraic[62] = (((0.293800*constants[6])/(constants[6]+210.000))*(states[0]-algebraic[59]))/(1.00000+exp(0.0896000*(states[0]-algebraic[59])))
    algebraic[63] = constants[59]*(power(states[32], 2.00000))*(states[0]-algebraic[59])
    algebraic[64] = constants[60]*states[33]*states[34]*(states[0]-algebraic[59])
    algebraic[65] = constants[61]*states[35]*states[36]*(states[0]-algebraic[59])
    algebraic[66] = constants[62]*states[37]*(states[0]-((constants[9]*constants[10])/constants[11])*log((0.980000*constants[6]+0.0200000*constants[7])/(0.980000*states[27]+0.0200000*states[18])))
    algebraic[0] = custom_piecewise([greater_equal(voi , constants[12]) & less_equal(voi , constants[13]) & less_equal((voi-constants[12])-floor((voi-constants[12])/constants[14])*constants[14] , constants[15]), constants[16] , True, 0.00000])
    algebraic[69] = 0.200000/(1.00000+exp(-(states[0]-46.7000)/7.80000))
    algebraic[70] = ((constants[68]*algebraic[69]*states[1])/(states[1]+constants[70]))*(states[0]-constants[69])
    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)