# Size of variable arrays:
sizeAlgebraic = 124
sizeStates = 19
sizeConstants = 58
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] = "R in component environment (J_per_K_per_mol)"
    legend_constants[1] = "T in component environment (kelvin)"
    legend_constants[2] = "F in component environment (C_per_mol)"
    legend_constants[3] = "C_m in component environment (fF)"
    legend_states[0] = "q_Na_o in component environment (fmol)"
    legend_states[1] = "q_Na_i in component environment (fmol)"
    legend_states[2] = "q_S000_Na in component environment (fmol)"
    legend_states[3] = "q_S010_Na in component environment (fmol)"
    legend_states[4] = "q_S100_Na in component environment (fmol)"
    legend_states[5] = "q_S110_Na in component environment (fmol)"
    legend_states[6] = "q_S200_Na in component environment (fmol)"
    legend_states[7] = "q_S210_Na in component environment (fmol)"
    legend_states[8] = "q_S300_Na in component environment (fmol)"
    legend_states[9] = "q_S310_Na in component environment (fmol)"
    legend_states[10] = "q_S001_Na in component environment (fmol)"
    legend_states[11] = "q_S011_Na in component environment (fmol)"
    legend_states[12] = "q_S101_Na in component environment (fmol)"
    legend_states[13] = "q_S111_Na in component environment (fmol)"
    legend_states[14] = "q_S201_Na in component environment (fmol)"
    legend_states[15] = "q_S211_Na in component environment (fmol)"
    legend_states[16] = "q_S301_Na in component environment (fmol)"
    legend_states[17] = "q_S311_Na in component environment (fmol)"
    legend_algebraic[120] = "v_Na in component fast_Na (fmol_per_sec)"
    legend_states[18] = "q_mem in component environment (fC)"
    legend_algebraic[122] = "I_mem_Na in component fast_Na (fA)"
    legend_constants[4] = "K_000_Na in component Na_parameters (per_fmol)"
    legend_constants[5] = "K_001_Na in component Na_parameters (per_fmol)"
    legend_constants[6] = "K_010_Na in component Na_parameters (per_fmol)"
    legend_constants[7] = "K_011_Na in component Na_parameters (per_fmol)"
    legend_constants[8] = "K_100_Na in component Na_parameters (per_fmol)"
    legend_constants[9] = "K_101_Na in component Na_parameters (per_fmol)"
    legend_constants[10] = "K_110_Na in component Na_parameters (per_fmol)"
    legend_constants[11] = "K_111_Na in component Na_parameters (per_fmol)"
    legend_constants[12] = "K_200_Na in component Na_parameters (per_fmol)"
    legend_constants[13] = "K_201_Na in component Na_parameters (per_fmol)"
    legend_constants[14] = "K_210_Na in component Na_parameters (per_fmol)"
    legend_constants[15] = "K_211_Na in component Na_parameters (per_fmol)"
    legend_constants[16] = "K_300_Na in component Na_parameters (per_fmol)"
    legend_constants[17] = "K_301_Na in component Na_parameters (per_fmol)"
    legend_constants[18] = "K_310_Na in component Na_parameters (per_fmol)"
    legend_constants[19] = "K_311_Na in component Na_parameters (per_fmol)"
    legend_constants[20] = "K_Na_o in component Na_parameters (per_fmol)"
    legend_constants[21] = "K_Na_i in component Na_parameters (per_fmol)"
    legend_constants[22] = "kappa_Na in component Na_parameters (fmol_per_sec)"
    legend_constants[23] = "kappa_h000 in component Na_parameters (fmol_per_sec)"
    legend_constants[24] = "kappa_h001 in component Na_parameters (fmol_per_sec)"
    legend_constants[25] = "kappa_h100 in component Na_parameters (fmol_per_sec)"
    legend_constants[26] = "kappa_h101 in component Na_parameters (fmol_per_sec)"
    legend_constants[27] = "kappa_h200 in component Na_parameters (fmol_per_sec)"
    legend_constants[28] = "kappa_h201 in component Na_parameters (fmol_per_sec)"
    legend_constants[29] = "kappa_h300 in component Na_parameters (fmol_per_sec)"
    legend_constants[30] = "kappa_h301 in component Na_parameters (fmol_per_sec)"
    legend_constants[31] = "kappa_j000 in component Na_parameters (fmol_per_sec)"
    legend_constants[32] = "kappa_j010 in component Na_parameters (fmol_per_sec)"
    legend_constants[33] = "kappa_j100 in component Na_parameters (fmol_per_sec)"
    legend_constants[34] = "kappa_j110 in component Na_parameters (fmol_per_sec)"
    legend_constants[35] = "kappa_j200 in component Na_parameters (fmol_per_sec)"
    legend_constants[36] = "kappa_j210 in component Na_parameters (fmol_per_sec)"
    legend_constants[37] = "kappa_j300 in component Na_parameters (fmol_per_sec)"
    legend_constants[38] = "kappa_j310 in component Na_parameters (fmol_per_sec)"
    legend_constants[39] = "kappa_m000 in component Na_parameters (fmol_per_sec)"
    legend_constants[40] = "kappa_m001 in component Na_parameters (fmol_per_sec)"
    legend_constants[41] = "kappa_m010 in component Na_parameters (fmol_per_sec)"
    legend_constants[42] = "kappa_m011 in component Na_parameters (fmol_per_sec)"
    legend_constants[43] = "kappa_m100 in component Na_parameters (fmol_per_sec)"
    legend_constants[44] = "kappa_m101 in component Na_parameters (fmol_per_sec)"
    legend_constants[45] = "kappa_m110 in component Na_parameters (fmol_per_sec)"
    legend_constants[46] = "kappa_m111 in component Na_parameters (fmol_per_sec)"
    legend_constants[47] = "kappa_m200 in component Na_parameters (fmol_per_sec)"
    legend_constants[48] = "kappa_m201 in component Na_parameters (fmol_per_sec)"
    legend_constants[49] = "kappa_m210 in component Na_parameters (fmol_per_sec)"
    legend_constants[50] = "kappa_m211 in component Na_parameters (fmol_per_sec)"
    legend_constants[51] = "zNa in component Na_parameters (dimensionless)"
    legend_constants[52] = "z_fh in component Na_parameters (dimensionless)"
    legend_constants[53] = "z_fj in component Na_parameters (dimensionless)"
    legend_constants[54] = "z_fm in component Na_parameters (dimensionless)"
    legend_constants[55] = "z_rh in component Na_parameters (dimensionless)"
    legend_constants[56] = "z_rj in component Na_parameters (dimensionless)"
    legend_constants[57] = "z_rm in component Na_parameters (dimensionless)"
    legend_algebraic[16] = "mu_Na_o in component fast_Na (J_per_mol)"
    legend_algebraic[17] = "mu_Na_i in component fast_Na (J_per_mol)"
    legend_algebraic[1] = "mu_S000_Na in component fast_Na (J_per_mol)"
    legend_algebraic[2] = "mu_S010_Na in component fast_Na (J_per_mol)"
    legend_algebraic[3] = "mu_S100_Na in component fast_Na (J_per_mol)"
    legend_algebraic[4] = "mu_S110_Na in component fast_Na (J_per_mol)"
    legend_algebraic[5] = "mu_S200_Na in component fast_Na (J_per_mol)"
    legend_algebraic[6] = "mu_S210_Na in component fast_Na (J_per_mol)"
    legend_algebraic[7] = "mu_S300_Na in component fast_Na (J_per_mol)"
    legend_algebraic[8] = "mu_S310_Na in component fast_Na (J_per_mol)"
    legend_algebraic[9] = "mu_S001_Na in component fast_Na (J_per_mol)"
    legend_algebraic[10] = "mu_S011_Na in component fast_Na (J_per_mol)"
    legend_algebraic[11] = "mu_S101_Na in component fast_Na (J_per_mol)"
    legend_algebraic[12] = "mu_S111_Na in component fast_Na (J_per_mol)"
    legend_algebraic[13] = "mu_S201_Na in component fast_Na (J_per_mol)"
    legend_algebraic[14] = "mu_S211_Na in component fast_Na (J_per_mol)"
    legend_algebraic[15] = "mu_S301_Na in component fast_Na (J_per_mol)"
    legend_algebraic[18] = "mu_S311_Na in component fast_Na (J_per_mol)"
    legend_algebraic[48] = "Am_Na in component fast_Na (J_per_mol)"
    legend_algebraic[19] = "Af_Na in component fast_Na (J_per_mol)"
    legend_algebraic[49] = "Ar_Na in component fast_Na (J_per_mol)"
    legend_algebraic[20] = "Af_h000 in component fast_Na (J_per_mol)"
    legend_algebraic[50] = "Ar_h000 in component fast_Na (J_per_mol)"
    legend_algebraic[51] = "v_h000 in component fast_Na (fmol_per_sec)"
    legend_algebraic[21] = "Af_h001 in component fast_Na (J_per_mol)"
    legend_algebraic[52] = "Ar_h001 in component fast_Na (J_per_mol)"
    legend_algebraic[53] = "v_h001 in component fast_Na (fmol_per_sec)"
    legend_algebraic[22] = "Af_h100 in component fast_Na (J_per_mol)"
    legend_algebraic[54] = "Ar_h100 in component fast_Na (J_per_mol)"
    legend_algebraic[55] = "v_h100 in component fast_Na (fmol_per_sec)"
    legend_algebraic[23] = "Af_h101 in component fast_Na (J_per_mol)"
    legend_algebraic[56] = "Ar_h101 in component fast_Na (J_per_mol)"
    legend_algebraic[57] = "v_h101 in component fast_Na (fmol_per_sec)"
    legend_algebraic[24] = "Af_h200 in component fast_Na (J_per_mol)"
    legend_algebraic[58] = "Ar_h200 in component fast_Na (J_per_mol)"
    legend_algebraic[59] = "v_h200 in component fast_Na (fmol_per_sec)"
    legend_algebraic[25] = "Af_h201 in component fast_Na (J_per_mol)"
    legend_algebraic[60] = "Ar_h201 in component fast_Na (J_per_mol)"
    legend_algebraic[61] = "v_h201 in component fast_Na (fmol_per_sec)"
    legend_algebraic[26] = "Af_h300 in component fast_Na (J_per_mol)"
    legend_algebraic[62] = "Ar_h300 in component fast_Na (J_per_mol)"
    legend_algebraic[63] = "v_h300 in component fast_Na (fmol_per_sec)"
    legend_algebraic[27] = "Af_h301 in component fast_Na (J_per_mol)"
    legend_algebraic[64] = "Ar_h301 in component fast_Na (J_per_mol)"
    legend_algebraic[65] = "v_h301 in component fast_Na (fmol_per_sec)"
    legend_algebraic[28] = "Af_j000 in component fast_Na (J_per_mol)"
    legend_algebraic[66] = "Ar_j000 in component fast_Na (J_per_mol)"
    legend_algebraic[67] = "v_j000 in component fast_Na (fmol_per_sec)"
    legend_algebraic[29] = "Af_j010 in component fast_Na (J_per_mol)"
    legend_algebraic[68] = "Ar_j010 in component fast_Na (J_per_mol)"
    legend_algebraic[69] = "v_j010 in component fast_Na (fmol_per_sec)"
    legend_algebraic[30] = "Af_j100 in component fast_Na (J_per_mol)"
    legend_algebraic[70] = "Ar_j100 in component fast_Na (J_per_mol)"
    legend_algebraic[71] = "v_j100 in component fast_Na (fmol_per_sec)"
    legend_algebraic[31] = "Af_j110 in component fast_Na (J_per_mol)"
    legend_algebraic[72] = "Ar_j110 in component fast_Na (J_per_mol)"
    legend_algebraic[73] = "v_j110 in component fast_Na (fmol_per_sec)"
    legend_algebraic[32] = "Af_j200 in component fast_Na (J_per_mol)"
    legend_algebraic[74] = "Ar_j200 in component fast_Na (J_per_mol)"
    legend_algebraic[75] = "v_j200 in component fast_Na (fmol_per_sec)"
    legend_algebraic[33] = "Af_j210 in component fast_Na (J_per_mol)"
    legend_algebraic[76] = "Ar_j210 in component fast_Na (J_per_mol)"
    legend_algebraic[77] = "v_j210 in component fast_Na (fmol_per_sec)"
    legend_algebraic[34] = "Af_j300 in component fast_Na (J_per_mol)"
    legend_algebraic[78] = "Ar_j300 in component fast_Na (J_per_mol)"
    legend_algebraic[79] = "v_j300 in component fast_Na (fmol_per_sec)"
    legend_algebraic[35] = "Af_j310 in component fast_Na (J_per_mol)"
    legend_algebraic[80] = "Ar_j310 in component fast_Na (J_per_mol)"
    legend_algebraic[81] = "v_j310 in component fast_Na (fmol_per_sec)"
    legend_algebraic[36] = "Af_m000 in component fast_Na (J_per_mol)"
    legend_algebraic[82] = "Ar_m000 in component fast_Na (J_per_mol)"
    legend_algebraic[83] = "v_m000 in component fast_Na (fmol_per_sec)"
    legend_algebraic[37] = "Af_m001 in component fast_Na (J_per_mol)"
    legend_algebraic[84] = "Ar_m001 in component fast_Na (J_per_mol)"
    legend_algebraic[86] = "v_m001 in component fast_Na (fmol_per_sec)"
    legend_algebraic[38] = "Af_m010 in component fast_Na (J_per_mol)"
    legend_algebraic[87] = "Ar_m010 in component fast_Na (J_per_mol)"
    legend_algebraic[89] = "v_m010 in component fast_Na (fmol_per_sec)"
    legend_algebraic[39] = "Af_m011 in component fast_Na (J_per_mol)"
    legend_algebraic[90] = "Ar_m011 in component fast_Na (J_per_mol)"
    legend_algebraic[92] = "v_m011 in component fast_Na (fmol_per_sec)"
    legend_algebraic[40] = "Af_m100 in component fast_Na (J_per_mol)"
    legend_algebraic[93] = "Ar_m100 in component fast_Na (J_per_mol)"
    legend_algebraic[95] = "v_m100 in component fast_Na (fmol_per_sec)"
    legend_algebraic[41] = "Af_m101 in component fast_Na (J_per_mol)"
    legend_algebraic[96] = "Ar_m101 in component fast_Na (J_per_mol)"
    legend_algebraic[98] = "v_m101 in component fast_Na (fmol_per_sec)"
    legend_algebraic[42] = "Af_m110 in component fast_Na (J_per_mol)"
    legend_algebraic[99] = "Ar_m110 in component fast_Na (J_per_mol)"
    legend_algebraic[101] = "v_m110 in component fast_Na (fmol_per_sec)"
    legend_algebraic[43] = "Af_m111 in component fast_Na (J_per_mol)"
    legend_algebraic[102] = "Ar_m111 in component fast_Na (J_per_mol)"
    legend_algebraic[104] = "v_m111 in component fast_Na (fmol_per_sec)"
    legend_algebraic[44] = "Af_m200 in component fast_Na (J_per_mol)"
    legend_algebraic[105] = "Ar_m200 in component fast_Na (J_per_mol)"
    legend_algebraic[107] = "v_m200 in component fast_Na (fmol_per_sec)"
    legend_algebraic[45] = "Af_m201 in component fast_Na (J_per_mol)"
    legend_algebraic[108] = "Ar_m201 in component fast_Na (J_per_mol)"
    legend_algebraic[111] = "v_m201 in component fast_Na (fmol_per_sec)"
    legend_algebraic[46] = "Af_m210 in component fast_Na (J_per_mol)"
    legend_algebraic[112] = "Ar_m210 in component fast_Na (J_per_mol)"
    legend_algebraic[115] = "v_m210 in component fast_Na (fmol_per_sec)"
    legend_algebraic[47] = "Af_m211 in component fast_Na (J_per_mol)"
    legend_algebraic[116] = "Ar_m211 in component fast_Na (J_per_mol)"
    legend_algebraic[119] = "v_m211 in component fast_Na (fmol_per_sec)"
    legend_algebraic[85] = "v_S000_Na in component fast_Na (fmol_per_sec)"
    legend_algebraic[91] = "v_S010_Na in component fast_Na (fmol_per_sec)"
    legend_algebraic[97] = "v_S100_Na in component fast_Na (fmol_per_sec)"
    legend_algebraic[103] = "v_S110_Na in component fast_Na (fmol_per_sec)"
    legend_algebraic[109] = "v_S200_Na in component fast_Na (fmol_per_sec)"
    legend_algebraic[117] = "v_S210_Na in component fast_Na (fmol_per_sec)"
    legend_algebraic[110] = "v_S300_Na in component fast_Na (fmol_per_sec)"
    legend_algebraic[118] = "v_S310_Na in component fast_Na (fmol_per_sec)"
    legend_algebraic[88] = "v_S001_Na in component fast_Na (fmol_per_sec)"
    legend_algebraic[94] = "v_S011_Na in component fast_Na (fmol_per_sec)"
    legend_algebraic[100] = "v_S101_Na in component fast_Na (fmol_per_sec)"
    legend_algebraic[106] = "v_S111_Na in component fast_Na (fmol_per_sec)"
    legend_algebraic[113] = "v_S201_Na in component fast_Na (fmol_per_sec)"
    legend_algebraic[121] = "v_S211_Na in component fast_Na (fmol_per_sec)"
    legend_algebraic[114] = "v_S301_Na in component fast_Na (fmol_per_sec)"
    legend_algebraic[123] = "v_S311_Na in component fast_Na (fmol_per_sec)"
    legend_algebraic[0] = "V_mem in component fast_Na (volt)"
    legend_rates[0] = "d/dt q_Na_o in component environment (fmol)"
    legend_rates[1] = "d/dt q_Na_i in component environment (fmol)"
    legend_rates[18] = "d/dt q_mem in component environment (fC)"
    legend_rates[2] = "d/dt q_S000_Na in component environment (fmol)"
    legend_rates[3] = "d/dt q_S010_Na in component environment (fmol)"
    legend_rates[4] = "d/dt q_S100_Na in component environment (fmol)"
    legend_rates[5] = "d/dt q_S110_Na in component environment (fmol)"
    legend_rates[6] = "d/dt q_S200_Na in component environment (fmol)"
    legend_rates[7] = "d/dt q_S210_Na in component environment (fmol)"
    legend_rates[8] = "d/dt q_S300_Na in component environment (fmol)"
    legend_rates[9] = "d/dt q_S310_Na in component environment (fmol)"
    legend_rates[10] = "d/dt q_S001_Na in component environment (fmol)"
    legend_rates[11] = "d/dt q_S011_Na in component environment (fmol)"
    legend_rates[12] = "d/dt q_S101_Na in component environment (fmol)"
    legend_rates[13] = "d/dt q_S111_Na in component environment (fmol)"
    legend_rates[14] = "d/dt q_S201_Na in component environment (fmol)"
    legend_rates[15] = "d/dt q_S211_Na in component environment (fmol)"
    legend_rates[16] = "d/dt q_S301_Na in component environment (fmol)"
    legend_rates[17] = "d/dt q_S311_Na 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] = 8.314
    constants[1] = 310
    constants[2] = 96485
    constants[3] = 153400
    states[0] = 725.48
    states[1] = 380
    states[2] = 1.72572E-05
    states[3] = 1.91747E-06
    states[4] = 5.75256E-06
    states[5] = 6.39176E-07
    states[6] = 6.39176E-07
    states[7] = 7.10196E-08
    states[8] = 2.36737E-08
    states[9] = 2.63039E-09
    states[10] = 1.91747E-06
    states[11] = 2.13059E-07
    states[12] = 6.39176E-07
    states[13] = 7.10196E-08
    states[14] = 7.10196E-08
    states[15] = 7.89103E-09
    states[16] = 2.63039E-09
    states[17] = 2.92267E-10
    states[18] = -13039
    constants[4] = 0.0073376
    constants[5] = 23718.0222
    constants[6] = 17499.8767
    constants[7] = 56566767727.4036
    constants[8] = 1.5633e-05
    constants[9] = 50.5317
    constants[10] = 37.2839
    constants[11] = 120516694.0958
    constants[12] = 9.9918e-08
    constants[13] = 0.32298
    constants[14] = 0.2383
    constants[15] = 770290.0911
    constants[16] = 1.9159e-09
    constants[17] = 0.006193
    constants[18] = 0.0045694
    constants[19] = 14770.0739
    constants[20] = 0.017678
    constants[21] = 0.0024108
    constants[22] = 32.858
    constants[23] = 0.045709
    constants[24] = 1.4141e-08
    constants[25] = 21.4544
    constants[26] = 6.6373e-06
    constants[27] = 3356.675
    constants[28] = 0.0010384
    constants[29] = 175057.5863
    constants[30] = 0.054157
    constants[31] = 0.017824
    constants[32] = 7.4736e-09
    constants[33] = 8.3662
    constants[34] = 3.5079e-06
    constants[35] = 1308.9427
    constants[36] = 0.00054883
    constants[37] = 68264.0876
    constants[38] = 0.028623
    constants[39] = 5117405.6443
    constants[40] = 1.5832
    constants[41] = 2.1457
    constants[42] = 6.6381e-07
    constants[43] = 1601300127.6001
    constants[44] = 495.389
    constants[45] = 671.4132
    constants[46] = 0.00020771
    constants[47] = 125266701374.7363
    constants[48] = 38753.3515
    constants[49] = 52523.3901
    constants[50] = 0.016249
    constants[51] = 1
    constants[52] = -4.1892
    constants[53] = -4.0381
    constants[54] = 0.49541
    constants[55] = 1.2995
    constants[56] = 1.4281
    constants[57] = -2.4284
    return (states, constants)

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

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