# Size of variable arrays:
sizeAlgebraic = 123
sizeStates = 17
sizeConstants = 59
from math import *
from numpy import *

def createLegends():
    legend_states = [""] * sizeStates
    legend_rates = [""] * sizeStates
    legend_algebraic = [""] * sizeAlgebraic
    legend_voi = ""
    legend_constants = [""] * sizeConstants
    legend_voi = "t in component environment (second)"
    legend_constants[0] = "C_m in component environment (fF)"
    legend_states[0] = "q_Ca_o in component environment (fmol)"
    legend_states[1] = "q_Ca_i in component environment (fmol)"
    legend_states[2] = "q_K_o in component environment (fmol)"
    legend_states[3] = "q_K_i in component environment (fmol)"
    legend_states[4] = "q_S000_LCC in component environment (fmol)"
    legend_states[5] = "q_S010_LCC in component environment (fmol)"
    legend_states[6] = "q_S020_LCC in component environment (fmol)"
    legend_states[7] = "q_S100_LCC in component environment (fmol)"
    legend_states[8] = "q_S110_LCC in component environment (fmol)"
    legend_states[9] = "q_S120_LCC in component environment (fmol)"
    legend_states[10] = "q_S001_LCC in component environment (fmol)"
    legend_states[11] = "q_S011_LCC in component environment (fmol)"
    legend_states[12] = "q_S021_LCC in component environment (fmol)"
    legend_states[13] = "q_S101_LCC in component environment (fmol)"
    legend_states[14] = "q_S111_LCC in component environment (fmol)"
    legend_states[15] = "q_S121_LCC in component environment (fmol)"
    legend_states[16] = "q_mem in component environment (fC)"
    legend_algebraic[119] = "I_mem_LCC in component LCC (fA)"
    legend_algebraic[115] = "v_Ca_o_LCC in component LCC (fmol_per_sec)"
    legend_algebraic[116] = "v_Ca_i_LCC in component LCC (fmol_per_sec)"
    legend_algebraic[120] = "v_K_o_LCC in component LCC (fmol_per_sec)"
    legend_algebraic[121] = "v_K_i_LCC in component LCC (fmol_per_sec)"
    legend_constants[1] = "kappa_LCC_Ca1 in component LCC_parameters (fmol_per_sec)"
    legend_constants[2] = "kappa_LCC_Ca2 in component LCC_parameters (fmol_per_sec)"
    legend_constants[3] = "kappa_LCC_K1 in component LCC_parameters (fmol_per_sec)"
    legend_constants[4] = "kappa_LCC_K2 in component LCC_parameters (fmol_per_sec)"
    legend_constants[5] = "kappa_d000 in component LCC_parameters (fmol_per_sec)"
    legend_constants[6] = "kappa_d010 in component LCC_parameters (fmol_per_sec)"
    legend_constants[7] = "kappa_d020 in component LCC_parameters (fmol_per_sec)"
    legend_constants[8] = "kappa_d001 in component LCC_parameters (fmol_per_sec)"
    legend_constants[9] = "kappa_d011 in component LCC_parameters (fmol_per_sec)"
    legend_constants[10] = "kappa_d021 in component LCC_parameters (fmol_per_sec)"
    legend_constants[11] = "kappa_f1_000 in component LCC_parameters (fmol_per_sec)"
    legend_constants[12] = "kappa_f1_100 in component LCC_parameters (fmol_per_sec)"
    legend_constants[13] = "kappa_f1_001 in component LCC_parameters (fmol_per_sec)"
    legend_constants[14] = "kappa_f1_101 in component LCC_parameters (fmol_per_sec)"
    legend_constants[15] = "kappa_f2_000 in component LCC_parameters (fmol_per_sec)"
    legend_constants[16] = "kappa_f2_100 in component LCC_parameters (fmol_per_sec)"
    legend_constants[17] = "kappa_f2_001 in component LCC_parameters (fmol_per_sec)"
    legend_constants[18] = "kappa_f2_101 in component LCC_parameters (fmol_per_sec)"
    legend_constants[19] = "kappa_f3_010 in component LCC_parameters (fmol_per_sec)"
    legend_constants[20] = "kappa_f3_110 in component LCC_parameters (fmol_per_sec)"
    legend_constants[21] = "kappa_f3_011 in component LCC_parameters (fmol_per_sec)"
    legend_constants[22] = "kappa_f3_111 in component LCC_parameters (fmol_per_sec)"
    legend_constants[23] = "kappa_fCa000 in component LCC_parameters (fmol_per_sec)"
    legend_constants[24] = "kappa_fCa100 in component LCC_parameters (fmol_per_sec)"
    legend_constants[25] = "kappa_fCa010 in component LCC_parameters (fmol_per_sec)"
    legend_constants[26] = "kappa_fCa110 in component LCC_parameters (fmol_per_sec)"
    legend_constants[27] = "kappa_fCa020 in component LCC_parameters (fmol_per_sec)"
    legend_constants[28] = "kappa_fCa120 in component LCC_parameters (fmol_per_sec)"
    legend_constants[29] = "K_Ca_i in component LCC_parameters (per_fmol)"
    legend_constants[30] = "K_Ca_o in component LCC_parameters (per_fmol)"
    legend_constants[31] = "K_K_i in component LCC_parameters (per_fmol)"
    legend_constants[32] = "K_K_o in component LCC_parameters (per_fmol)"
    legend_constants[33] = "K_000_LCC in component LCC_parameters (per_fmol)"
    legend_constants[34] = "K_010_LCC in component LCC_parameters (per_fmol)"
    legend_constants[35] = "K_020_LCC in component LCC_parameters (per_fmol)"
    legend_constants[36] = "K_100_LCC in component LCC_parameters (per_fmol)"
    legend_constants[37] = "K_110_LCC in component LCC_parameters (per_fmol)"
    legend_constants[38] = "K_120_LCC in component LCC_parameters (per_fmol)"
    legend_constants[39] = "K_001_LCC in component LCC_parameters (per_fmol)"
    legend_constants[40] = "K_011_LCC in component LCC_parameters (per_fmol)"
    legend_constants[41] = "K_021_LCC in component LCC_parameters (per_fmol)"
    legend_constants[42] = "K_101_LCC in component LCC_parameters (per_fmol)"
    legend_constants[43] = "K_111_LCC in component LCC_parameters (per_fmol)"
    legend_constants[44] = "K_121_LCC in component LCC_parameters (per_fmol)"
    legend_constants[45] = "zCa in component LCC_parameters (dimensionless)"
    legend_constants[46] = "zK in component LCC_parameters (dimensionless)"
    legend_constants[47] = "z_fd in component LCC_parameters (dimensionless)"
    legend_constants[48] = "z_ff1 in component LCC_parameters (dimensionless)"
    legend_constants[49] = "z_ff2 in component LCC_parameters (dimensionless)"
    legend_constants[50] = "z_ff3 in component LCC_parameters (dimensionless)"
    legend_constants[51] = "z_rCa in component LCC_parameters (dimensionless)"
    legend_constants[52] = "z_rd in component LCC_parameters (dimensionless)"
    legend_constants[53] = "z_rf1 in component LCC_parameters (dimensionless)"
    legend_constants[54] = "z_rf2 in component LCC_parameters (dimensionless)"
    legend_constants[55] = "z_rf3 in component LCC_parameters (dimensionless)"
    legend_constants[56] = "R in component constants (J_per_K_per_mol)"
    legend_constants[57] = "T in component constants (kelvin)"
    legend_constants[58] = "F in component constants (C_per_mol)"
    legend_algebraic[6] = "mu_Ca_o in component LCC (J_per_mol)"
    legend_algebraic[7] = "mu_Ca_i in component LCC (J_per_mol)"
    legend_algebraic[8] = "mu_K_o in component LCC (J_per_mol)"
    legend_algebraic[9] = "mu_K_i in component LCC (J_per_mol)"
    legend_algebraic[0] = "V_mem in component LCC (volt)"
    legend_algebraic[18] = "Af_LCC_Ca1 in component LCC (J_per_mol)"
    legend_algebraic[50] = "Ar_LCC_Ca1 in component LCC (J_per_mol)"
    legend_algebraic[46] = "Am_LCC_Ca1 in component LCC (J_per_mol)"
    legend_algebraic[111] = "v_LCC_Ca1 in component LCC (fmol_per_sec)"
    legend_algebraic[19] = "Af_LCC_Ca2 in component LCC (J_per_mol)"
    legend_algebraic[51] = "Ar_LCC_Ca2 in component LCC (J_per_mol)"
    legend_algebraic[47] = "Am_LCC_Ca2 in component LCC (J_per_mol)"
    legend_algebraic[113] = "v_LCC_Ca2 in component LCC (fmol_per_sec)"
    legend_algebraic[20] = "Af_LCC_K1 in component LCC (J_per_mol)"
    legend_algebraic[52] = "Ar_LCC_K1 in component LCC (J_per_mol)"
    legend_algebraic[48] = "Am_LCC_K1 in component LCC (J_per_mol)"
    legend_algebraic[114] = "v_LCC_K1 in component LCC (fmol_per_sec)"
    legend_algebraic[21] = "Af_LCC_K2 in component LCC (J_per_mol)"
    legend_algebraic[53] = "Ar_LCC_K2 in component LCC (J_per_mol)"
    legend_algebraic[49] = "Am_LCC_K2 in component LCC (J_per_mol)"
    legend_algebraic[117] = "v_LCC_K2 in component LCC (fmol_per_sec)"
    legend_algebraic[40] = "Af_fCa00 in component LCC (J_per_mol)"
    legend_algebraic[90] = "Ar_fCa00 in component LCC (J_per_mol)"
    legend_algebraic[91] = "v_fCa00 in component LCC (fmol_per_sec)"
    legend_algebraic[41] = "Af_fCa01 in component LCC (J_per_mol)"
    legend_algebraic[92] = "Ar_fCa01 in component LCC (J_per_mol)"
    legend_algebraic[95] = "v_fCa01 in component LCC (fmol_per_sec)"
    legend_algebraic[42] = "Af_fCa02 in component LCC (J_per_mol)"
    legend_algebraic[96] = "Ar_fCa02 in component LCC (J_per_mol)"
    legend_algebraic[99] = "v_fCa02 in component LCC (fmol_per_sec)"
    legend_algebraic[43] = "Af_fCa10 in component LCC (J_per_mol)"
    legend_algebraic[100] = "Ar_fCa10 in component LCC (J_per_mol)"
    legend_algebraic[103] = "v_fCa10 in component LCC (fmol_per_sec)"
    legend_algebraic[44] = "Af_fCa11 in component LCC (J_per_mol)"
    legend_algebraic[104] = "Ar_fCa11 in component LCC (J_per_mol)"
    legend_algebraic[107] = "v_fCa11 in component LCC (fmol_per_sec)"
    legend_algebraic[45] = "Af_fCa12 in component LCC (J_per_mol)"
    legend_algebraic[108] = "Ar_fCa12 in component LCC (J_per_mol)"
    legend_algebraic[110] = "v_fCa12 in component LCC (fmol_per_sec)"
    legend_algebraic[10] = "mu_S000_LCC in component LCC (J_per_mol)"
    legend_algebraic[93] = "v_S000_LCC in component LCC (fmol_per_sec)"
    legend_algebraic[11] = "mu_S010_LCC in component LCC (J_per_mol)"
    legend_algebraic[97] = "v_S010_LCC in component LCC (fmol_per_sec)"
    legend_algebraic[12] = "mu_S020_LCC in component LCC (J_per_mol)"
    legend_algebraic[101] = "v_S020_LCC in component LCC (fmol_per_sec)"
    legend_algebraic[13] = "mu_S100_LCC in component LCC (J_per_mol)"
    legend_algebraic[105] = "v_S100_LCC in component LCC (fmol_per_sec)"
    legend_algebraic[14] = "mu_S110_LCC in component LCC (J_per_mol)"
    legend_algebraic[109] = "v_S110_LCC in component LCC (fmol_per_sec)"
    legend_algebraic[16] = "mu_S120_LCC in component LCC (J_per_mol)"
    legend_algebraic[112] = "v_S120_LCC in component LCC (fmol_per_sec)"
    legend_algebraic[1] = "mu_S001_LCC in component LCC (J_per_mol)"
    legend_algebraic[94] = "v_S001_LCC in component LCC (fmol_per_sec)"
    legend_algebraic[2] = "mu_S011_LCC in component LCC (J_per_mol)"
    legend_algebraic[98] = "v_S011_LCC in component LCC (fmol_per_sec)"
    legend_algebraic[3] = "mu_S021_LCC in component LCC (J_per_mol)"
    legend_algebraic[102] = "v_S021_LCC in component LCC (fmol_per_sec)"
    legend_algebraic[4] = "mu_S101_LCC in component LCC (J_per_mol)"
    legend_algebraic[106] = "v_S101_LCC in component LCC (fmol_per_sec)"
    legend_algebraic[15] = "mu_S111_LCC in component LCC (J_per_mol)"
    legend_algebraic[118] = "v_S111_LCC in component LCC (fmol_per_sec)"
    legend_algebraic[17] = "mu_S121_LCC in component LCC (J_per_mol)"
    legend_algebraic[122] = "v_S121_LCC in component LCC (fmol_per_sec)"
    legend_algebraic[22] = "Af_d000 in component LCC (J_per_mol)"
    legend_algebraic[54] = "Ar_d000 in component LCC (J_per_mol)"
    legend_algebraic[55] = "v_d000 in component LCC (fmol_per_sec)"
    legend_algebraic[24] = "Af_d010 in component LCC (J_per_mol)"
    legend_algebraic[58] = "Ar_d010 in component LCC (J_per_mol)"
    legend_algebraic[59] = "v_d010 in component LCC (fmol_per_sec)"
    legend_algebraic[26] = "Af_d020 in component LCC (J_per_mol)"
    legend_algebraic[62] = "Ar_d020 in component LCC (J_per_mol)"
    legend_algebraic[63] = "v_d020 in component LCC (fmol_per_sec)"
    legend_algebraic[28] = "Af_f1_000 in component LCC (J_per_mol)"
    legend_algebraic[66] = "Ar_f1_000 in component LCC (J_per_mol)"
    legend_algebraic[67] = "v_f1_000 in component LCC (fmol_per_sec)"
    legend_algebraic[32] = "Af_f2_000 in component LCC (J_per_mol)"
    legend_algebraic[74] = "Ar_f2_000 in component LCC (J_per_mol)"
    legend_algebraic[75] = "v_f2_000 in component LCC (fmol_per_sec)"
    legend_algebraic[36] = "Af_f3_010 in component LCC (J_per_mol)"
    legend_algebraic[82] = "Ar_f3_010 in component LCC (J_per_mol)"
    legend_algebraic[83] = "v_f3_010 in component LCC (fmol_per_sec)"
    legend_algebraic[30] = "Af_f1_100 in component LCC (J_per_mol)"
    legend_algebraic[70] = "Ar_f1_100 in component LCC (J_per_mol)"
    legend_algebraic[71] = "v_f1_100 in component LCC (fmol_per_sec)"
    legend_algebraic[34] = "Af_f2_100 in component LCC (J_per_mol)"
    legend_algebraic[78] = "Ar_f2_100 in component LCC (J_per_mol)"
    legend_algebraic[79] = "v_f2_100 in component LCC (fmol_per_sec)"
    legend_algebraic[38] = "Af_f3_110 in component LCC (J_per_mol)"
    legend_algebraic[86] = "Ar_f3_110 in component LCC (J_per_mol)"
    legend_algebraic[87] = "v_f3_110 in component LCC (fmol_per_sec)"
    legend_algebraic[23] = "Af_d001 in component LCC (J_per_mol)"
    legend_algebraic[56] = "Ar_d001 in component LCC (J_per_mol)"
    legend_algebraic[57] = "v_d001 in component LCC (fmol_per_sec)"
    legend_algebraic[25] = "Af_d011 in component LCC (J_per_mol)"
    legend_algebraic[60] = "Ar_d011 in component LCC (J_per_mol)"
    legend_algebraic[61] = "v_d011 in component LCC (fmol_per_sec)"
    legend_algebraic[27] = "Af_d021 in component LCC (J_per_mol)"
    legend_algebraic[64] = "Ar_d021 in component LCC (J_per_mol)"
    legend_algebraic[65] = "v_d021 in component LCC (fmol_per_sec)"
    legend_algebraic[29] = "Af_f1_001 in component LCC (J_per_mol)"
    legend_algebraic[68] = "Ar_f1_001 in component LCC (J_per_mol)"
    legend_algebraic[69] = "v_f1_001 in component LCC (fmol_per_sec)"
    legend_algebraic[33] = "Af_f2_001 in component LCC (J_per_mol)"
    legend_algebraic[76] = "Ar_f2_001 in component LCC (J_per_mol)"
    legend_algebraic[77] = "v_f2_001 in component LCC (fmol_per_sec)"
    legend_algebraic[37] = "Af_f3_011 in component LCC (J_per_mol)"
    legend_algebraic[84] = "Ar_f3_011 in component LCC (J_per_mol)"
    legend_algebraic[85] = "v_f3_011 in component LCC (fmol_per_sec)"
    legend_algebraic[31] = "Af_f1_101 in component LCC (J_per_mol)"
    legend_algebraic[72] = "Ar_f1_101 in component LCC (J_per_mol)"
    legend_algebraic[73] = "v_f1_101 in component LCC (fmol_per_sec)"
    legend_algebraic[35] = "Af_f2_101 in component LCC (J_per_mol)"
    legend_algebraic[80] = "Ar_f2_101 in component LCC (J_per_mol)"
    legend_algebraic[81] = "v_f2_101 in component LCC (fmol_per_sec)"
    legend_algebraic[39] = "Af_f3_111 in component LCC (J_per_mol)"
    legend_algebraic[88] = "Ar_f3_111 in component LCC (J_per_mol)"
    legend_algebraic[89] = "v_f3_111 in component LCC (fmol_per_sec)"
    legend_algebraic[5] = "Ca_tot in component LCC (fmol)"
    legend_rates[0] = "d/dt q_Ca_o in component environment (fmol)"
    legend_rates[1] = "d/dt q_Ca_i in component environment (fmol)"
    legend_rates[2] = "d/dt q_K_o in component environment (fmol)"
    legend_rates[3] = "d/dt q_K_i in component environment (fmol)"
    legend_rates[16] = "d/dt q_mem in component environment (fC)"
    legend_rates[4] = "d/dt q_S000_LCC in component environment (fmol)"
    legend_rates[5] = "d/dt q_S010_LCC in component environment (fmol)"
    legend_rates[6] = "d/dt q_S020_LCC in component environment (fmol)"
    legend_rates[7] = "d/dt q_S100_LCC in component environment (fmol)"
    legend_rates[8] = "d/dt q_S110_LCC in component environment (fmol)"
    legend_rates[9] = "d/dt q_S120_LCC in component environment (fmol)"
    legend_rates[10] = "d/dt q_S001_LCC in component environment (fmol)"
    legend_rates[11] = "d/dt q_S011_LCC in component environment (fmol)"
    legend_rates[12] = "d/dt q_S021_LCC in component environment (fmol)"
    legend_rates[13] = "d/dt q_S101_LCC in component environment (fmol)"
    legend_rates[14] = "d/dt q_S111_LCC in component environment (fmol)"
    legend_rates[15] = "d/dt q_S121_LCC in component environment (fmol)"
    return (legend_states, legend_algebraic, legend_voi, legend_constants)

def initConsts():
    constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
    constants[0] = 60000
    states[0] = 0.245463158
    states[1] = 6e-5
    states[2] = 27.9828
    states[3] = 5510
    states[4] = 5.30481E-08
    states[5] = 1.05036E-05
    states[6] = 5.30481E-08
    states[7] = 5.35845E-10
    states[8] = 1.06098E-07
    states[9] = 5.35845E-10
    states[10] = 5.89425E-09
    states[11] = 1.16708E-06
    states[12] = 5.89425E-09
    states[13] = 5.95377E-11
    states[14] = 1.17886E-08
    states[15] = 5.95377E-11
    states[16] = -13039
    constants[1] = 74.6745
    constants[2] = 98.3196
    constants[3] = 0.131085
    constants[4] = 0.172592
    constants[5] = 53.6616
    constants[6] = 2.7099
    constants[7] = 3.56797
    constants[8] = 0.28538
    constants[9] = 0.0144117
    constants[10] = 0.018975
    constants[11] = 5.16892
    constants[12] = 25.6676
    constants[13] = 0.0274891
    constants[14] = 0.136504
    constants[15] = 0.335961
    constants[16] = 1.6683
    constants[17] = 0.00178669
    constants[18] = 0.00887227
    constants[19] = 733.001
    constants[20] = 3639.9
    constants[21] = 3.89821
    constants[22] = 19.3575
    constants[23] = 97713.8
    constants[24] = 485222
    constants[25] = 4934.53
    constants[26] = 24503.7
    constants[27] = 6497.01
    constants[28] = 32262.5
    constants[29] = 0.090536
    constants[30] = 0.663931
    constants[31] = 0.0405426
    constants[32] = 0.297313
    constants[33] = 9.07096
    constants[34] = 179.624
    constants[35] = 136.425
    constants[36] = 1.8267
    constants[37] = 36.1725
    constants[38] = 27.4733
    constants[39] = 1705.66
    constants[40] = 33775.6
    constants[41] = 25652.8
    constants[42] = 343.485
    constants[43] = 6801.71
    constants[44] = 5165.95
    constants[45] = 2
    constants[46] = 1
    constants[47] = 2.1404
    constants[48] = -1.1495
    constants[49] = 0.72162
    constants[50] = 4.2933
    constants[51] = 2
    constants[52] = -2.1404
    constants[53] = 1.8993
    constants[54] = -0.52288
    constants[55] = 0
    constants[56] = 8.31
    constants[57] = 310
    constants[58] = 96485
    return (states, constants)

def computeRates(voi, states, constants):
    rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
    algebraic[10] = constants[56]*constants[57]*log(constants[33]*states[4])
    algebraic[40] = algebraic[10]
    algebraic[7] = constants[56]*constants[57]*log(constants[29]*states[1])
    algebraic[90] = constants[56]*constants[57]*log(constants[39]*states[10])+constants[51]*algebraic[7]
    algebraic[91] = constants[23]*(exp(algebraic[40]/(constants[56]*constants[57]))-exp(algebraic[90]/(constants[56]*constants[57])))
    algebraic[0] = states[16]/constants[0]
    algebraic[22] = constants[56]*constants[57]*log(constants[33]*states[4])+constants[47]*constants[58]*algebraic[0]
    algebraic[54] = constants[56]*constants[57]*log(constants[36]*states[7])+constants[52]*constants[58]*algebraic[0]
    algebraic[55] = constants[5]*(exp(algebraic[22]/(constants[56]*constants[57]))-exp(algebraic[54]/(constants[56]*constants[57])))
    algebraic[28] = constants[56]*constants[57]*log(constants[33]*states[4])+constants[48]*constants[58]*algebraic[0]
    algebraic[66] = constants[56]*constants[57]*log(constants[34]*states[5])+constants[53]*constants[58]*algebraic[0]
    algebraic[67] = constants[11]*(exp(algebraic[28]/(constants[56]*constants[57]))-exp(algebraic[66]/(constants[56]*constants[57])))
    algebraic[32] = constants[56]*constants[57]*log(constants[33]*states[4])+constants[49]*constants[58]*algebraic[0]
    algebraic[74] = constants[56]*constants[57]*log(constants[35]*states[6])+constants[54]*constants[58]*algebraic[0]
    algebraic[75] = constants[15]*(exp(algebraic[32]/(constants[56]*constants[57]))-exp(algebraic[74]/(constants[56]*constants[57])))
    algebraic[93] = ((-algebraic[55]-algebraic[67])-algebraic[75])-algebraic[91]
    rates[4] = algebraic[93]
    algebraic[23] = constants[56]*constants[57]*log(constants[39]*states[10])+constants[47]*constants[58]*algebraic[0]
    algebraic[56] = constants[56]*constants[57]*log(constants[42]*states[13])+constants[52]*constants[58]*algebraic[0]
    algebraic[57] = constants[8]*(exp(algebraic[23]/(constants[56]*constants[57]))-exp(algebraic[56]/(constants[56]*constants[57])))
    algebraic[29] = constants[56]*constants[57]*log(constants[39]*states[10])+constants[48]*constants[58]*algebraic[0]
    algebraic[68] = constants[56]*constants[57]*log(constants[40]*states[11])+constants[53]*constants[58]*algebraic[0]
    algebraic[69] = constants[13]*(exp(algebraic[29]/(constants[56]*constants[57]))-exp(algebraic[68]/(constants[56]*constants[57])))
    algebraic[33] = constants[56]*constants[57]*log(constants[39]*states[10])+constants[49]*constants[58]*algebraic[0]
    algebraic[76] = constants[56]*constants[57]*log(constants[41]*states[12])+constants[54]*constants[58]*algebraic[0]
    algebraic[77] = constants[17]*(exp(algebraic[33]/(constants[56]*constants[57]))-exp(algebraic[76]/(constants[56]*constants[57])))
    algebraic[94] = ((-algebraic[57]-algebraic[69])-algebraic[77])--algebraic[91]
    rates[10] = algebraic[94]
    algebraic[11] = constants[56]*constants[57]*log(constants[34]*states[5])
    algebraic[41] = algebraic[11]
    algebraic[92] = constants[56]*constants[57]*log(constants[40]*states[11])+constants[51]*algebraic[7]
    algebraic[95] = constants[25]*(exp(algebraic[41]/(constants[56]*constants[57]))-exp(algebraic[92]/(constants[56]*constants[57])))
    algebraic[24] = constants[56]*constants[57]*log(constants[34]*states[5])+constants[47]*constants[58]*algebraic[0]
    algebraic[58] = constants[56]*constants[57]*log(constants[37]*states[8])+constants[52]*constants[58]*algebraic[0]
    algebraic[59] = constants[6]*(exp(algebraic[24]/(constants[56]*constants[57]))-exp(algebraic[58]/(constants[56]*constants[57])))
    algebraic[36] = constants[56]*constants[57]*log(constants[34]*states[5])+constants[50]*constants[58]*algebraic[0]
    algebraic[82] = constants[56]*constants[57]*log(constants[35]*states[6])+constants[55]*constants[58]*algebraic[0]
    algebraic[83] = constants[19]*(exp(algebraic[36]/(constants[56]*constants[57]))-exp(algebraic[82]/(constants[56]*constants[57])))
    algebraic[97] = ((-algebraic[59]+algebraic[67])-algebraic[95])-algebraic[83]
    rates[5] = algebraic[97]
    algebraic[25] = constants[56]*constants[57]*log(constants[40]*states[11])+constants[47]*constants[58]*algebraic[0]
    algebraic[60] = constants[56]*constants[57]*log(constants[43]*states[14])+constants[52]*constants[58]*algebraic[0]
    algebraic[61] = constants[9]*(exp(algebraic[25]/(constants[56]*constants[57]))-exp(algebraic[60]/(constants[56]*constants[57])))
    algebraic[37] = constants[56]*constants[57]*log(constants[40]*states[11])+constants[50]*constants[58]*algebraic[0]
    algebraic[84] = constants[56]*constants[57]*log(constants[41]*states[12])+constants[55]*constants[58]*algebraic[0]
    algebraic[85] = constants[21]*(exp(algebraic[37]/(constants[56]*constants[57]))-exp(algebraic[84]/(constants[56]*constants[57])))
    algebraic[98] = ((-algebraic[61]+algebraic[69])--algebraic[95])-algebraic[85]
    rates[11] = algebraic[98]
    algebraic[12] = constants[56]*constants[57]*log(constants[35]*states[6])
    algebraic[42] = algebraic[12]
    algebraic[96] = constants[56]*constants[57]*log(constants[41]*states[12])+constants[51]*algebraic[7]
    algebraic[99] = constants[27]*(exp(algebraic[42]/(constants[56]*constants[57]))-exp(algebraic[96]/(constants[56]*constants[57])))
    algebraic[26] = constants[56]*constants[57]*log(constants[35]*states[6])+constants[47]*constants[58]*algebraic[0]
    algebraic[62] = constants[56]*constants[57]*log(constants[38]*states[9])+constants[52]*constants[58]*algebraic[0]
    algebraic[63] = constants[7]*(exp(algebraic[26]/(constants[56]*constants[57]))-exp(algebraic[62]/(constants[56]*constants[57])))
    algebraic[101] = ((-algebraic[63]+algebraic[75])-algebraic[99])+algebraic[83]
    rates[6] = algebraic[101]
    algebraic[27] = constants[56]*constants[57]*log(constants[41]*states[12])+constants[47]*constants[58]*algebraic[0]
    algebraic[64] = constants[56]*constants[57]*log(constants[44]*states[15])+constants[52]*constants[58]*algebraic[0]
    algebraic[65] = constants[10]*(exp(algebraic[27]/(constants[56]*constants[57]))-exp(algebraic[64]/(constants[56]*constants[57])))
    algebraic[102] = ((-algebraic[65]+algebraic[77])--algebraic[99])+algebraic[85]
    rates[12] = algebraic[102]
    algebraic[13] = constants[56]*constants[57]*log(constants[36]*states[7])
    algebraic[43] = algebraic[13]
    algebraic[100] = constants[56]*constants[57]*log(constants[42]*states[13])+constants[51]*algebraic[7]
    algebraic[103] = constants[24]*(exp(algebraic[43]/(constants[56]*constants[57]))-exp(algebraic[100]/(constants[56]*constants[57])))
    algebraic[30] = constants[56]*constants[57]*log(constants[36]*states[7])+constants[48]*constants[58]*algebraic[0]
    algebraic[70] = constants[56]*constants[57]*log(constants[37]*states[8])+constants[53]*constants[58]*algebraic[0]
    algebraic[71] = constants[12]*(exp(algebraic[30]/(constants[56]*constants[57]))-exp(algebraic[70]/(constants[56]*constants[57])))
    algebraic[34] = constants[56]*constants[57]*log(constants[36]*states[7])+constants[49]*constants[58]*algebraic[0]
    algebraic[78] = constants[56]*constants[57]*log(constants[38]*states[9])+constants[54]*constants[58]*algebraic[0]
    algebraic[79] = constants[16]*(exp(algebraic[34]/(constants[56]*constants[57]))-exp(algebraic[78]/(constants[56]*constants[57])))
    algebraic[105] = ((-algebraic[79]+algebraic[55])-algebraic[71])-algebraic[103]
    rates[7] = algebraic[105]
    algebraic[31] = constants[56]*constants[57]*log(constants[42]*states[13])+constants[48]*constants[58]*algebraic[0]
    algebraic[72] = constants[56]*constants[57]*log(constants[43]*states[14])+constants[53]*constants[58]*algebraic[0]
    algebraic[73] = constants[14]*(exp(algebraic[31]/(constants[56]*constants[57]))-exp(algebraic[72]/(constants[56]*constants[57])))
    algebraic[35] = constants[56]*constants[57]*log(constants[42]*states[13])+constants[49]*constants[58]*algebraic[0]
    algebraic[80] = constants[56]*constants[57]*log(constants[44]*states[15])+constants[54]*constants[58]*algebraic[0]
    algebraic[81] = constants[18]*(exp(algebraic[35]/(constants[56]*constants[57]))-exp(algebraic[80]/(constants[56]*constants[57])))
    algebraic[106] = ((-algebraic[81]+algebraic[57])-algebraic[73])--algebraic[103]
    rates[13] = algebraic[106]
    algebraic[14] = constants[56]*constants[57]*log(constants[37]*states[8])
    algebraic[44] = algebraic[14]
    algebraic[104] = constants[56]*constants[57]*log(constants[43]*states[14])+constants[51]*algebraic[7]
    algebraic[107] = constants[26]*(exp(algebraic[44]/(constants[56]*constants[57]))-exp(algebraic[104]/(constants[56]*constants[57])))
    algebraic[38] = constants[56]*constants[57]*log(constants[37]*states[8])+constants[50]*constants[58]*algebraic[0]
    algebraic[86] = constants[56]*constants[57]*log(constants[38]*states[9])+constants[55]*constants[58]*algebraic[0]
    algebraic[87] = constants[20]*(exp(algebraic[38]/(constants[56]*constants[57]))-exp(algebraic[86]/(constants[56]*constants[57])))
    algebraic[109] = ((algebraic[71]+algebraic[59])-algebraic[107])-algebraic[87]
    rates[8] = algebraic[109]
    algebraic[16] = constants[56]*constants[57]*log(constants[38]*states[9])
    algebraic[45] = algebraic[16]
    algebraic[108] = constants[56]*constants[57]*log(constants[44]*states[15])+constants[51]*algebraic[7]
    algebraic[110] = constants[28]*(exp(algebraic[45]/(constants[56]*constants[57]))-exp(algebraic[108]/(constants[56]*constants[57])))
    algebraic[112] = ((algebraic[79]+algebraic[63])-algebraic[110])+algebraic[87]
    rates[9] = algebraic[112]
    algebraic[15] = constants[56]*constants[57]*log(constants[43]*states[14])
    algebraic[18] = algebraic[7]+constants[45]*constants[58]*algebraic[0]+algebraic[15]
    algebraic[6] = constants[56]*constants[57]*log(constants[30]*states[0])
    algebraic[50] = algebraic[6]+algebraic[15]
    algebraic[46] = constants[45]*constants[58]*algebraic[0]
    algebraic[111] = custom_piecewise([equal(algebraic[46] , 0.00000), constants[1]*(exp(algebraic[18]/(constants[56]*constants[57]))-exp(algebraic[50]/(constants[56]*constants[57]))) , True, (((constants[1]*algebraic[46])/(constants[56]*constants[57]))/(exp(algebraic[46]/(constants[56]*constants[57]))-1.00000))*(exp(algebraic[18]/(constants[56]*constants[57]))-exp(algebraic[50]/(constants[56]*constants[57])))])
    algebraic[17] = constants[56]*constants[57]*log(constants[44]*states[15])
    algebraic[19] = algebraic[7]+constants[45]*constants[58]*algebraic[0]+algebraic[17]
    algebraic[51] = algebraic[6]+algebraic[17]
    algebraic[47] = constants[45]*constants[58]*algebraic[0]
    algebraic[113] = custom_piecewise([equal(algebraic[47] , 0.00000), constants[2]*(exp(algebraic[19]/(constants[56]*constants[57]))-exp(algebraic[51]/(constants[56]*constants[57]))) , True, (((constants[2]*algebraic[47])/(constants[56]*constants[57]))/(exp(algebraic[47]/(constants[56]*constants[57]))-1.00000))*(exp(algebraic[19]/(constants[56]*constants[57]))-exp(algebraic[51]/(constants[56]*constants[57])))])
    algebraic[115] = algebraic[113]+algebraic[111]
    rates[0] = algebraic[115]
    algebraic[116] = (((constants[51])*algebraic[103]+constants[51]*algebraic[107]+constants[51]*algebraic[99]+constants[51]*algebraic[95]+constants[51]*algebraic[91]+constants[51]*algebraic[110])-algebraic[111])-algebraic[113]
    rates[1] = algebraic[116]
    algebraic[9] = constants[56]*constants[57]*log(constants[31]*states[3])
    algebraic[20] = algebraic[9]+constants[46]*constants[58]*algebraic[0]+algebraic[15]
    algebraic[8] = constants[56]*constants[57]*log(constants[32]*states[2])
    algebraic[52] = algebraic[8]+algebraic[15]
    algebraic[48] = constants[46]*constants[58]*algebraic[0]
    algebraic[114] = custom_piecewise([equal(algebraic[48] , 0.00000), constants[3]*(exp(algebraic[20]/(constants[56]*constants[57]))-exp(algebraic[52]/(constants[56]*constants[57]))) , True, (((constants[3]*algebraic[48])/(constants[56]*constants[57]))/(exp(algebraic[48]/(constants[56]*constants[57]))-1.00000))*(exp(algebraic[20]/(constants[56]*constants[57]))-exp(algebraic[52]/(constants[56]*constants[57])))])
    algebraic[39] = constants[56]*constants[57]*log(constants[43]*states[14])+constants[50]*constants[58]*algebraic[0]
    algebraic[88] = constants[56]*constants[57]*log(constants[44]*states[15])+constants[55]*constants[58]*algebraic[0]
    algebraic[89] = constants[22]*(exp(algebraic[39]/(constants[56]*constants[57]))-exp(algebraic[88]/(constants[56]*constants[57])))
    algebraic[118] = ((algebraic[73]+algebraic[61])-((((-algebraic[107]+algebraic[111])-algebraic[111])+algebraic[114])-algebraic[114]))-algebraic[89]
    rates[14] = algebraic[118]
    algebraic[21] = algebraic[9]+constants[46]*constants[58]*algebraic[0]+algebraic[17]
    algebraic[53] = algebraic[8]+algebraic[17]
    algebraic[49] = constants[46]*constants[58]*algebraic[0]
    algebraic[117] = custom_piecewise([equal(algebraic[49] , 0.00000), constants[4]*(exp(algebraic[21]/(constants[56]*constants[57]))-exp(algebraic[53]/(constants[56]*constants[57]))) , True, (((constants[4]*algebraic[49])/(constants[56]*constants[57]))/(exp(algebraic[49]/(constants[56]*constants[57]))-1.00000))*(exp(algebraic[21]/(constants[56]*constants[57]))-exp(algebraic[53]/(constants[56]*constants[57])))])
    algebraic[120] = algebraic[117]+algebraic[114]
    rates[2] = algebraic[120]
    algebraic[121] = -algebraic[114]-algebraic[117]
    rates[3] = algebraic[121]
    algebraic[119] = constants[58]*(((((((((((((((((((((((((((((((((((((constants[54]*algebraic[77]-constants[49]*algebraic[77])+constants[55]*algebraic[85])-constants[50]*algebraic[85])+constants[53]*algebraic[73])-constants[48]*algebraic[73])+constants[54]*algebraic[81])-constants[49]*algebraic[81])+constants[55]*algebraic[89])-constants[50]*algebraic[89])+constants[53]*algebraic[69])-constants[48]*algebraic[69])+constants[52]*algebraic[65])-constants[47]*algebraic[65])+constants[52]*algebraic[61])-constants[47]*algebraic[61])+constants[52]*algebraic[57])-constants[47]*algebraic[57])+constants[55]*algebraic[87])-constants[50]*algebraic[87])+constants[54]*algebraic[79])-constants[49]*algebraic[79])+constants[53]*algebraic[71])-constants[48]*algebraic[71])+constants[55]*algebraic[83])-constants[50]*algebraic[83])+constants[54]*algebraic[75])-constants[49]*algebraic[75])+constants[53]*algebraic[67])-constants[48]*algebraic[67])+constants[52]*algebraic[63])-constants[47]*algebraic[63])+constants[52]*algebraic[59])-constants[47]*algebraic[59])+constants[52]*algebraic[55])-constants[47]*algebraic[55])-(constants[45]*algebraic[111]+constants[45]*algebraic[113]))-(constants[46]*algebraic[114]+constants[46]*algebraic[117]))
    rates[16] = algebraic[119]
    algebraic[122] = ((algebraic[81]+algebraic[65])-((((-algebraic[110]+algebraic[113])-algebraic[113])+algebraic[117])-algebraic[117]))+algebraic[89]
    rates[15] = algebraic[122]
    return(rates)

def computeAlgebraic(constants, states, voi):
    algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
    states = array(states)
    voi = array(voi)
    algebraic[10] = constants[56]*constants[57]*log(constants[33]*states[4])
    algebraic[40] = algebraic[10]
    algebraic[7] = constants[56]*constants[57]*log(constants[29]*states[1])
    algebraic[90] = constants[56]*constants[57]*log(constants[39]*states[10])+constants[51]*algebraic[7]
    algebraic[91] = constants[23]*(exp(algebraic[40]/(constants[56]*constants[57]))-exp(algebraic[90]/(constants[56]*constants[57])))
    algebraic[0] = states[16]/constants[0]
    algebraic[22] = constants[56]*constants[57]*log(constants[33]*states[4])+constants[47]*constants[58]*algebraic[0]
    algebraic[54] = constants[56]*constants[57]*log(constants[36]*states[7])+constants[52]*constants[58]*algebraic[0]
    algebraic[55] = constants[5]*(exp(algebraic[22]/(constants[56]*constants[57]))-exp(algebraic[54]/(constants[56]*constants[57])))
    algebraic[28] = constants[56]*constants[57]*log(constants[33]*states[4])+constants[48]*constants[58]*algebraic[0]
    algebraic[66] = constants[56]*constants[57]*log(constants[34]*states[5])+constants[53]*constants[58]*algebraic[0]
    algebraic[67] = constants[11]*(exp(algebraic[28]/(constants[56]*constants[57]))-exp(algebraic[66]/(constants[56]*constants[57])))
    algebraic[32] = constants[56]*constants[57]*log(constants[33]*states[4])+constants[49]*constants[58]*algebraic[0]
    algebraic[74] = constants[56]*constants[57]*log(constants[35]*states[6])+constants[54]*constants[58]*algebraic[0]
    algebraic[75] = constants[15]*(exp(algebraic[32]/(constants[56]*constants[57]))-exp(algebraic[74]/(constants[56]*constants[57])))
    algebraic[93] = ((-algebraic[55]-algebraic[67])-algebraic[75])-algebraic[91]
    algebraic[23] = constants[56]*constants[57]*log(constants[39]*states[10])+constants[47]*constants[58]*algebraic[0]
    algebraic[56] = constants[56]*constants[57]*log(constants[42]*states[13])+constants[52]*constants[58]*algebraic[0]
    algebraic[57] = constants[8]*(exp(algebraic[23]/(constants[56]*constants[57]))-exp(algebraic[56]/(constants[56]*constants[57])))
    algebraic[29] = constants[56]*constants[57]*log(constants[39]*states[10])+constants[48]*constants[58]*algebraic[0]
    algebraic[68] = constants[56]*constants[57]*log(constants[40]*states[11])+constants[53]*constants[58]*algebraic[0]
    algebraic[69] = constants[13]*(exp(algebraic[29]/(constants[56]*constants[57]))-exp(algebraic[68]/(constants[56]*constants[57])))
    algebraic[33] = constants[56]*constants[57]*log(constants[39]*states[10])+constants[49]*constants[58]*algebraic[0]
    algebraic[76] = constants[56]*constants[57]*log(constants[41]*states[12])+constants[54]*constants[58]*algebraic[0]
    algebraic[77] = constants[17]*(exp(algebraic[33]/(constants[56]*constants[57]))-exp(algebraic[76]/(constants[56]*constants[57])))
    algebraic[94] = ((-algebraic[57]-algebraic[69])-algebraic[77])--algebraic[91]
    algebraic[11] = constants[56]*constants[57]*log(constants[34]*states[5])
    algebraic[41] = algebraic[11]
    algebraic[92] = constants[56]*constants[57]*log(constants[40]*states[11])+constants[51]*algebraic[7]
    algebraic[95] = constants[25]*(exp(algebraic[41]/(constants[56]*constants[57]))-exp(algebraic[92]/(constants[56]*constants[57])))
    algebraic[24] = constants[56]*constants[57]*log(constants[34]*states[5])+constants[47]*constants[58]*algebraic[0]
    algebraic[58] = constants[56]*constants[57]*log(constants[37]*states[8])+constants[52]*constants[58]*algebraic[0]
    algebraic[59] = constants[6]*(exp(algebraic[24]/(constants[56]*constants[57]))-exp(algebraic[58]/(constants[56]*constants[57])))
    algebraic[36] = constants[56]*constants[57]*log(constants[34]*states[5])+constants[50]*constants[58]*algebraic[0]
    algebraic[82] = constants[56]*constants[57]*log(constants[35]*states[6])+constants[55]*constants[58]*algebraic[0]
    algebraic[83] = constants[19]*(exp(algebraic[36]/(constants[56]*constants[57]))-exp(algebraic[82]/(constants[56]*constants[57])))
    algebraic[97] = ((-algebraic[59]+algebraic[67])-algebraic[95])-algebraic[83]
    algebraic[25] = constants[56]*constants[57]*log(constants[40]*states[11])+constants[47]*constants[58]*algebraic[0]
    algebraic[60] = constants[56]*constants[57]*log(constants[43]*states[14])+constants[52]*constants[58]*algebraic[0]
    algebraic[61] = constants[9]*(exp(algebraic[25]/(constants[56]*constants[57]))-exp(algebraic[60]/(constants[56]*constants[57])))
    algebraic[37] = constants[56]*constants[57]*log(constants[40]*states[11])+constants[50]*constants[58]*algebraic[0]
    algebraic[84] = constants[56]*constants[57]*log(constants[41]*states[12])+constants[55]*constants[58]*algebraic[0]
    algebraic[85] = constants[21]*(exp(algebraic[37]/(constants[56]*constants[57]))-exp(algebraic[84]/(constants[56]*constants[57])))
    algebraic[98] = ((-algebraic[61]+algebraic[69])--algebraic[95])-algebraic[85]
    algebraic[12] = constants[56]*constants[57]*log(constants[35]*states[6])
    algebraic[42] = algebraic[12]
    algebraic[96] = constants[56]*constants[57]*log(constants[41]*states[12])+constants[51]*algebraic[7]
    algebraic[99] = constants[27]*(exp(algebraic[42]/(constants[56]*constants[57]))-exp(algebraic[96]/(constants[56]*constants[57])))
    algebraic[26] = constants[56]*constants[57]*log(constants[35]*states[6])+constants[47]*constants[58]*algebraic[0]
    algebraic[62] = constants[56]*constants[57]*log(constants[38]*states[9])+constants[52]*constants[58]*algebraic[0]
    algebraic[63] = constants[7]*(exp(algebraic[26]/(constants[56]*constants[57]))-exp(algebraic[62]/(constants[56]*constants[57])))
    algebraic[101] = ((-algebraic[63]+algebraic[75])-algebraic[99])+algebraic[83]
    algebraic[27] = constants[56]*constants[57]*log(constants[41]*states[12])+constants[47]*constants[58]*algebraic[0]
    algebraic[64] = constants[56]*constants[57]*log(constants[44]*states[15])+constants[52]*constants[58]*algebraic[0]
    algebraic[65] = constants[10]*(exp(algebraic[27]/(constants[56]*constants[57]))-exp(algebraic[64]/(constants[56]*constants[57])))
    algebraic[102] = ((-algebraic[65]+algebraic[77])--algebraic[99])+algebraic[85]
    algebraic[13] = constants[56]*constants[57]*log(constants[36]*states[7])
    algebraic[43] = algebraic[13]
    algebraic[100] = constants[56]*constants[57]*log(constants[42]*states[13])+constants[51]*algebraic[7]
    algebraic[103] = constants[24]*(exp(algebraic[43]/(constants[56]*constants[57]))-exp(algebraic[100]/(constants[56]*constants[57])))
    algebraic[30] = constants[56]*constants[57]*log(constants[36]*states[7])+constants[48]*constants[58]*algebraic[0]
    algebraic[70] = constants[56]*constants[57]*log(constants[37]*states[8])+constants[53]*constants[58]*algebraic[0]
    algebraic[71] = constants[12]*(exp(algebraic[30]/(constants[56]*constants[57]))-exp(algebraic[70]/(constants[56]*constants[57])))
    algebraic[34] = constants[56]*constants[57]*log(constants[36]*states[7])+constants[49]*constants[58]*algebraic[0]
    algebraic[78] = constants[56]*constants[57]*log(constants[38]*states[9])+constants[54]*constants[58]*algebraic[0]
    algebraic[79] = constants[16]*(exp(algebraic[34]/(constants[56]*constants[57]))-exp(algebraic[78]/(constants[56]*constants[57])))
    algebraic[105] = ((-algebraic[79]+algebraic[55])-algebraic[71])-algebraic[103]
    algebraic[31] = constants[56]*constants[57]*log(constants[42]*states[13])+constants[48]*constants[58]*algebraic[0]
    algebraic[72] = constants[56]*constants[57]*log(constants[43]*states[14])+constants[53]*constants[58]*algebraic[0]
    algebraic[73] = constants[14]*(exp(algebraic[31]/(constants[56]*constants[57]))-exp(algebraic[72]/(constants[56]*constants[57])))
    algebraic[35] = constants[56]*constants[57]*log(constants[42]*states[13])+constants[49]*constants[58]*algebraic[0]
    algebraic[80] = constants[56]*constants[57]*log(constants[44]*states[15])+constants[54]*constants[58]*algebraic[0]
    algebraic[81] = constants[18]*(exp(algebraic[35]/(constants[56]*constants[57]))-exp(algebraic[80]/(constants[56]*constants[57])))
    algebraic[106] = ((-algebraic[81]+algebraic[57])-algebraic[73])--algebraic[103]
    algebraic[14] = constants[56]*constants[57]*log(constants[37]*states[8])
    algebraic[44] = algebraic[14]
    algebraic[104] = constants[56]*constants[57]*log(constants[43]*states[14])+constants[51]*algebraic[7]
    algebraic[107] = constants[26]*(exp(algebraic[44]/(constants[56]*constants[57]))-exp(algebraic[104]/(constants[56]*constants[57])))
    algebraic[38] = constants[56]*constants[57]*log(constants[37]*states[8])+constants[50]*constants[58]*algebraic[0]
    algebraic[86] = constants[56]*constants[57]*log(constants[38]*states[9])+constants[55]*constants[58]*algebraic[0]
    algebraic[87] = constants[20]*(exp(algebraic[38]/(constants[56]*constants[57]))-exp(algebraic[86]/(constants[56]*constants[57])))
    algebraic[109] = ((algebraic[71]+algebraic[59])-algebraic[107])-algebraic[87]
    algebraic[16] = constants[56]*constants[57]*log(constants[38]*states[9])
    algebraic[45] = algebraic[16]
    algebraic[108] = constants[56]*constants[57]*log(constants[44]*states[15])+constants[51]*algebraic[7]
    algebraic[110] = constants[28]*(exp(algebraic[45]/(constants[56]*constants[57]))-exp(algebraic[108]/(constants[56]*constants[57])))
    algebraic[112] = ((algebraic[79]+algebraic[63])-algebraic[110])+algebraic[87]
    algebraic[15] = constants[56]*constants[57]*log(constants[43]*states[14])
    algebraic[18] = algebraic[7]+constants[45]*constants[58]*algebraic[0]+algebraic[15]
    algebraic[6] = constants[56]*constants[57]*log(constants[30]*states[0])
    algebraic[50] = algebraic[6]+algebraic[15]
    algebraic[46] = constants[45]*constants[58]*algebraic[0]
    algebraic[111] = custom_piecewise([equal(algebraic[46] , 0.00000), constants[1]*(exp(algebraic[18]/(constants[56]*constants[57]))-exp(algebraic[50]/(constants[56]*constants[57]))) , True, (((constants[1]*algebraic[46])/(constants[56]*constants[57]))/(exp(algebraic[46]/(constants[56]*constants[57]))-1.00000))*(exp(algebraic[18]/(constants[56]*constants[57]))-exp(algebraic[50]/(constants[56]*constants[57])))])
    algebraic[17] = constants[56]*constants[57]*log(constants[44]*states[15])
    algebraic[19] = algebraic[7]+constants[45]*constants[58]*algebraic[0]+algebraic[17]
    algebraic[51] = algebraic[6]+algebraic[17]
    algebraic[47] = constants[45]*constants[58]*algebraic[0]
    algebraic[113] = custom_piecewise([equal(algebraic[47] , 0.00000), constants[2]*(exp(algebraic[19]/(constants[56]*constants[57]))-exp(algebraic[51]/(constants[56]*constants[57]))) , True, (((constants[2]*algebraic[47])/(constants[56]*constants[57]))/(exp(algebraic[47]/(constants[56]*constants[57]))-1.00000))*(exp(algebraic[19]/(constants[56]*constants[57]))-exp(algebraic[51]/(constants[56]*constants[57])))])
    algebraic[115] = algebraic[113]+algebraic[111]
    algebraic[116] = (((constants[51])*algebraic[103]+constants[51]*algebraic[107]+constants[51]*algebraic[99]+constants[51]*algebraic[95]+constants[51]*algebraic[91]+constants[51]*algebraic[110])-algebraic[111])-algebraic[113]
    algebraic[9] = constants[56]*constants[57]*log(constants[31]*states[3])
    algebraic[20] = algebraic[9]+constants[46]*constants[58]*algebraic[0]+algebraic[15]
    algebraic[8] = constants[56]*constants[57]*log(constants[32]*states[2])
    algebraic[52] = algebraic[8]+algebraic[15]
    algebraic[48] = constants[46]*constants[58]*algebraic[0]
    algebraic[114] = custom_piecewise([equal(algebraic[48] , 0.00000), constants[3]*(exp(algebraic[20]/(constants[56]*constants[57]))-exp(algebraic[52]/(constants[56]*constants[57]))) , True, (((constants[3]*algebraic[48])/(constants[56]*constants[57]))/(exp(algebraic[48]/(constants[56]*constants[57]))-1.00000))*(exp(algebraic[20]/(constants[56]*constants[57]))-exp(algebraic[52]/(constants[56]*constants[57])))])
    algebraic[39] = constants[56]*constants[57]*log(constants[43]*states[14])+constants[50]*constants[58]*algebraic[0]
    algebraic[88] = constants[56]*constants[57]*log(constants[44]*states[15])+constants[55]*constants[58]*algebraic[0]
    algebraic[89] = constants[22]*(exp(algebraic[39]/(constants[56]*constants[57]))-exp(algebraic[88]/(constants[56]*constants[57])))
    algebraic[118] = ((algebraic[73]+algebraic[61])-((((-algebraic[107]+algebraic[111])-algebraic[111])+algebraic[114])-algebraic[114]))-algebraic[89]
    algebraic[21] = algebraic[9]+constants[46]*constants[58]*algebraic[0]+algebraic[17]
    algebraic[53] = algebraic[8]+algebraic[17]
    algebraic[49] = constants[46]*constants[58]*algebraic[0]
    algebraic[117] = custom_piecewise([equal(algebraic[49] , 0.00000), constants[4]*(exp(algebraic[21]/(constants[56]*constants[57]))-exp(algebraic[53]/(constants[56]*constants[57]))) , True, (((constants[4]*algebraic[49])/(constants[56]*constants[57]))/(exp(algebraic[49]/(constants[56]*constants[57]))-1.00000))*(exp(algebraic[21]/(constants[56]*constants[57]))-exp(algebraic[53]/(constants[56]*constants[57])))])
    algebraic[120] = algebraic[117]+algebraic[114]
    algebraic[121] = -algebraic[114]-algebraic[117]
    algebraic[119] = constants[58]*(((((((((((((((((((((((((((((((((((((constants[54]*algebraic[77]-constants[49]*algebraic[77])+constants[55]*algebraic[85])-constants[50]*algebraic[85])+constants[53]*algebraic[73])-constants[48]*algebraic[73])+constants[54]*algebraic[81])-constants[49]*algebraic[81])+constants[55]*algebraic[89])-constants[50]*algebraic[89])+constants[53]*algebraic[69])-constants[48]*algebraic[69])+constants[52]*algebraic[65])-constants[47]*algebraic[65])+constants[52]*algebraic[61])-constants[47]*algebraic[61])+constants[52]*algebraic[57])-constants[47]*algebraic[57])+constants[55]*algebraic[87])-constants[50]*algebraic[87])+constants[54]*algebraic[79])-constants[49]*algebraic[79])+constants[53]*algebraic[71])-constants[48]*algebraic[71])+constants[55]*algebraic[83])-constants[50]*algebraic[83])+constants[54]*algebraic[75])-constants[49]*algebraic[75])+constants[53]*algebraic[67])-constants[48]*algebraic[67])+constants[52]*algebraic[63])-constants[47]*algebraic[63])+constants[52]*algebraic[59])-constants[47]*algebraic[59])+constants[52]*algebraic[55])-constants[47]*algebraic[55])-(constants[45]*algebraic[111]+constants[45]*algebraic[113]))-(constants[46]*algebraic[114]+constants[46]*algebraic[117]))
    algebraic[122] = ((algebraic[81]+algebraic[65])-((((-algebraic[110]+algebraic[113])-algebraic[113])+algebraic[117])-algebraic[117]))+algebraic[89]
    algebraic[1] = constants[56]*constants[57]*log(constants[39]*states[10])
    algebraic[2] = constants[56]*constants[57]*log(constants[40]*states[11])
    algebraic[3] = constants[56]*constants[57]*log(constants[41]*states[12])
    algebraic[4] = constants[56]*constants[57]*log(constants[42]*states[13])
    algebraic[5] = states[0]+states[1]+states[10]+states[11]+states[12]+states[13]+states[14]+states[15]
    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)