# Size of variable arrays:
sizeAlgebraic = 122
sizeStates = 2
sizeConstants = 34
from math import *
from numpy import *

def createLegends():
    legend_states = [""] * sizeStates
    legend_rates = [""] * sizeStates
    legend_algebraic = [""] * sizeAlgebraic
    legend_voi = ""
    legend_constants = [""] * sizeConstants
    legend_constants[0] = "k1_m in component gas_exchange (dimensionless)"
    legend_constants[1] = "k1_81m in component gas_exchange (dimensionless)"
    legend_constants[2] = "k2_m in component gas_exchange (dimensionless)"
    legend_constants[3] = "k3_m in component gas_exchange (dimensionless)"
    legend_constants[4] = "K3_81m in component gas_exchange (dimensionless)"
    legend_constants[5] = "k4_m in component gas_exchange (dimensionless)"
    legend_constants[6] = "pH_m in component gas_exchange (dimensionless)"
    legend_constants[7] = "k1_f in component gas_exchange (dimensionless)"
    legend_constants[8] = "k2_f in component gas_exchange (dimensionless)"
    legend_constants[9] = "k3_f in component gas_exchange (dimensionless)"
    legend_constants[10] = "k4_f in component gas_exchange (dimensionless)"
    legend_constants[11] = "pH_f in component gas_exchange (dimensionless)"
    legend_states[0] = "Pin_m in component gas_exchange (mmHg)"
    legend_states[1] = "Pin_f in component gas_exchange (mmHg)"
    legend_algebraic[0] = "x_m in component gas_exchange (dimensionless)"
    legend_algebraic[1] = "x_m1 in component gas_exchange (dimensionless)"
    legend_algebraic[3] = "x_f in component gas_exchange (dimensionless)"
    legend_algebraic[4] = "x_f1 in component gas_exchange (dimensionless)"
    legend_algebraic[2] = "S_m in component gas_exchange (dimensionless)"
    legend_algebraic[5] = "S_f in component gas_exchange (dimensionless)"
    legend_voi = "t in component gas_exchange (second)"
    legend_algebraic[6] = "O2_m in component gas_exchange (ml_per_ml)"
    legend_algebraic[7] = "O2_f in component gas_exchange (ml_per_ml)"
    legend_constants[12] = "Hb_m in component gas_exchange (dimensionless)"
    legend_constants[13] = "Hb_f in component gas_exchange (dimensionless)"
    legend_constants[14] = "a in component gas_exchange (dimensionless)"
    legend_constants[15] = "ks_m in component gas_exchange (dimensionless)"
    legend_constants[16] = "ks_f in component gas_exchange (dimensionless)"
    legend_constants[17] = "kh in component gas_exchange (dimensionless)"
    legend_constants[18] = "B in component gas_exchange (dimensionless)"
    legend_constants[19] = "V_m in component gas_exchange (ml)"
    legend_constants[20] = "V_f in component gas_exchange (ml)"
    legend_constants[21] = "Q_m in component gas_exchange (ml_per_min)"
    legend_constants[22] = "Q_f in component gas_exchange (ml_per_min)"
    legend_constants[23] = "Dm in component gas_exchange (ml_per_ml_min_mmHg)"
    legend_constants[24] = "c1 in component gas_exchange (dimensionless)"
    legend_constants[25] = "c2 in component gas_exchange (dimensionless)"
    legend_constants[26] = "c3 in component gas_exchange (dimensionless)"
    legend_constants[27] = "c4 in component gas_exchange (dimensionless)"
    legend_algebraic[10] = "theta_m in component gas_exchange (dimensionless)"
    legend_algebraic[13] = "theta_f in component gas_exchange (dimensionless)"
    legend_algebraic[8] = "x1_m in component gas_exchange (dimensionless)"
    legend_algebraic[11] = "x1_f in component gas_exchange (dimensionless)"
    legend_algebraic[14] = "R_m in component gas_exchange (dimensionless)"
    legend_algebraic[15] = "R_f in component gas_exchange (dimensionless)"
    legend_constants[28] = "R_p in component gas_exchange (dimensionless)"
    legend_algebraic[16] = "R_t in component gas_exchange (dimensionless)"
    legend_algebraic[17] = "Dp in component gas_exchange (dimensionless)"
    legend_algebraic[18] = "dodt_m in component gas_exchange (dimensionless)"
    legend_algebraic[19] = "dodt_f in component gas_exchange (dimensionless)"
    legend_algebraic[20] = "O2_mnew in component gas_exchange (ml_per_ml)"
    legend_algebraic[21] = "O2_fnew in component gas_exchange (ml_per_ml)"
    legend_algebraic[22] = "f_m in component gas_exchange (dimensionless)"
    legend_algebraic[23] = "f_f in component gas_exchange (dimensionless)"
    legend_algebraic[24] = "upper_m in component gas_exchange (dimensionless)"
    legend_algebraic[25] = "lower_m in component gas_exchange (dimensionless)"
    legend_algebraic[27] = "upper_f in component gas_exchange (dimensionless)"
    legend_algebraic[28] = "lower_f in component gas_exchange (dimensionless)"
    legend_algebraic[29] = "dfdP_f in component gas_exchange (dimensionless)"
    legend_algebraic[26] = "dfdP_m in component gas_exchange (dimensionless)"
    legend_constants[29] = "tt in component gas_exchange (second)"
    legend_constants[30] = "dt in component gas_exchange (second)"
    legend_algebraic[30] = "Pend1_m in component gas_exchange (mmHg)"
    legend_algebraic[31] = "Pend1_f in component gas_exchange (mmHg)"
    legend_algebraic[9] = "theta_96_m in component gas_exchange (dimensionless)"
    legend_algebraic[12] = "theta_96_f in component gas_exchange (dimensionless)"
    legend_algebraic[32] = "x_m_sec in component gas_exchange (dimensionless)"
    legend_algebraic[33] = "x_m1_sec in component gas_exchange (dimensionless)"
    legend_algebraic[34] = "S_m_sec in component gas_exchange (dimensionless)"
    legend_algebraic[35] = "x_f_sec in component gas_exchange (dimensionless)"
    legend_algebraic[36] = "x_f1_sec in component gas_exchange (dimensionless)"
    legend_algebraic[37] = "S_f_sec in component gas_exchange (dimensionless)"
    legend_algebraic[39] = "theta_96_m_sec in component gas_exchange (dimensionless)"
    legend_algebraic[42] = "theta_96_f_sec in component gas_exchange (dimensionless)"
    legend_algebraic[38] = "x1_m_sec in component gas_exchange (dimensionless)"
    legend_algebraic[40] = "theta_m_sec in component gas_exchange (dimensionless)"
    legend_algebraic[41] = "x1_f_sec in component gas_exchange (dimensionless)"
    legend_algebraic[43] = "theta_f_sec in component gas_exchange (dimensionless)"
    legend_algebraic[44] = "R_m_sec in component gas_exchange (dimensionless)"
    legend_algebraic[45] = "R_f_sec in component gas_exchange (dimensionless)"
    legend_constants[31] = "R_p_sec in component gas_exchange (dimensionless)"
    legend_algebraic[46] = "R_t_sec in component gas_exchange (dimensionless)"
    legend_algebraic[47] = "Dp_sec in component gas_exchange (dimensionless)"
    legend_algebraic[48] = "dodt_m_sec in component gas_exchange (dimensionless)"
    legend_algebraic[49] = "dodt_f_sec in component gas_exchange (dimensionless)"
    legend_algebraic[50] = "O2_mnew_sec in component gas_exchange (ml_per_ml)"
    legend_algebraic[51] = "O2_fnew_sec in component gas_exchange (ml_per_ml)"
    legend_algebraic[52] = "f_m_sec in component gas_exchange (dimensionless)"
    legend_algebraic[53] = "f_f_sec in component gas_exchange (dimensionless)"
    legend_algebraic[54] = "upper_m_sec in component gas_exchange (dimensionless)"
    legend_algebraic[55] = "lower_m_sec in component gas_exchange (dimensionless)"
    legend_algebraic[56] = "dfdP_m_sec in component gas_exchange (dimensionless)"
    legend_algebraic[57] = "upper_f_sec in component gas_exchange (dimensionless)"
    legend_algebraic[58] = "lower_f_sec in component gas_exchange (dimensionless)"
    legend_algebraic[59] = "dfdP_f_sec in component gas_exchange (dimensionless)"
    legend_algebraic[60] = "Pend2_m in component gas_exchange (mmHg)"
    legend_algebraic[61] = "Pend2_f in component gas_exchange (mmHg)"
    legend_algebraic[62] = "x_m_trd in component gas_exchange (dimensionless)"
    legend_algebraic[63] = "x_m1_trd in component gas_exchange (dimensionless)"
    legend_algebraic[64] = "S_m_trd in component gas_exchange (dimensionless)"
    legend_algebraic[65] = "x_f_trd in component gas_exchange (dimensionless)"
    legend_algebraic[66] = "x_f1_trd in component gas_exchange (dimensionless)"
    legend_algebraic[67] = "S_f_trd in component gas_exchange (dimensionless)"
    legend_algebraic[69] = "theta_96_m_trd in component gas_exchange (dimensionless)"
    legend_algebraic[72] = "theta_96_f_trd in component gas_exchange (dimensionless)"
    legend_algebraic[68] = "x1_m_trd in component gas_exchange (dimensionless)"
    legend_algebraic[70] = "theta_m_trd in component gas_exchange (dimensionless)"
    legend_algebraic[71] = "x1_f_trd in component gas_exchange (dimensionless)"
    legend_algebraic[73] = "theta_f_trd in component gas_exchange (dimensionless)"
    legend_algebraic[74] = "R_m_trd in component gas_exchange (dimensionless)"
    legend_algebraic[75] = "R_f_trd in component gas_exchange (dimensionless)"
    legend_constants[32] = "R_p_trd in component gas_exchange (dimensionless)"
    legend_algebraic[76] = "R_t_trd in component gas_exchange (dimensionless)"
    legend_algebraic[77] = "Dp_trd in component gas_exchange (dimensionless)"
    legend_algebraic[78] = "dodt_m_trd in component gas_exchange (dimensionless)"
    legend_algebraic[79] = "dodt_f_trd in component gas_exchange (dimensionless)"
    legend_algebraic[80] = "O2_mnew_trd in component gas_exchange (ml_per_ml)"
    legend_algebraic[81] = "O2_fnew_trd in component gas_exchange (ml_per_ml)"
    legend_algebraic[82] = "f_m_trd in component gas_exchange (dimensionless)"
    legend_algebraic[83] = "f_f_trd in component gas_exchange (dimensionless)"
    legend_algebraic[84] = "upper_m_trd in component gas_exchange (dimensionless)"
    legend_algebraic[85] = "lower_m_trd in component gas_exchange (dimensionless)"
    legend_algebraic[86] = "dfdP_m_trd in component gas_exchange (dimensionless)"
    legend_algebraic[87] = "upper_f_trd in component gas_exchange (dimensionless)"
    legend_algebraic[88] = "lower_f_trd in component gas_exchange (dimensionless)"
    legend_algebraic[89] = "dfdP_f_trd in component gas_exchange (dimensionless)"
    legend_algebraic[90] = "Pend3_m in component gas_exchange (mmHg)"
    legend_algebraic[91] = "Pend3_f in component gas_exchange (mmHg)"
    legend_algebraic[92] = "x_m_for in component gas_exchange (dimensionless)"
    legend_algebraic[93] = "x_m1_for in component gas_exchange (dimensionless)"
    legend_algebraic[94] = "S_m_for in component gas_exchange (dimensionless)"
    legend_algebraic[95] = "x_f_for in component gas_exchange (dimensionless)"
    legend_algebraic[96] = "x_f1_for in component gas_exchange (dimensionless)"
    legend_algebraic[97] = "S_f_for in component gas_exchange (dimensionless)"
    legend_algebraic[99] = "theta_96_m_for in component gas_exchange (dimensionless)"
    legend_algebraic[102] = "theta_96_f_for in component gas_exchange (dimensionless)"
    legend_algebraic[98] = "x1_m_for in component gas_exchange (dimensionless)"
    legend_algebraic[100] = "theta_m_for in component gas_exchange (dimensionless)"
    legend_algebraic[101] = "x1_f_for in component gas_exchange (dimensionless)"
    legend_algebraic[103] = "theta_f_for in component gas_exchange (dimensionless)"
    legend_algebraic[104] = "R_m_for in component gas_exchange (dimensionless)"
    legend_algebraic[105] = "R_f_for in component gas_exchange (dimensionless)"
    legend_constants[33] = "R_p_for in component gas_exchange (dimensionless)"
    legend_algebraic[106] = "R_t_for in component gas_exchange (dimensionless)"
    legend_algebraic[107] = "Dp_for in component gas_exchange (dimensionless)"
    legend_algebraic[108] = "dodt_m_for in component gas_exchange (dimensionless)"
    legend_algebraic[110] = "dodt_f_for in component gas_exchange (dimensionless)"
    legend_algebraic[109] = "O2_mnew_final in component gas_exchange (ml_per_ml)"
    legend_algebraic[112] = "O2_fnew_final in component gas_exchange (ml_per_ml)"
    legend_algebraic[111] = "f_m_final in component gas_exchange (dimensionless)"
    legend_algebraic[114] = "f_f_final in component gas_exchange (dimensionless)"
    legend_algebraic[113] = "upper_m_final in component gas_exchange (dimensionless)"
    legend_algebraic[115] = "lower_m_final in component gas_exchange (dimensionless)"
    legend_algebraic[117] = "dfdP_m_final in component gas_exchange (dimensionless)"
    legend_algebraic[116] = "upper_f_final in component gas_exchange (dimensionless)"
    legend_algebraic[118] = "lower_f_final in component gas_exchange (dimensionless)"
    legend_algebraic[119] = "dfdP_f_final in component gas_exchange (dimensionless)"
    legend_algebraic[120] = "Pfinal_m in component gas_exchange (mmHg)"
    legend_algebraic[121] = "Pfinal_f in component gas_exchange (mmHg)"
    legend_rates[0] = "d/dt Pin_m in component gas_exchange (mmHg)"
    legend_rates[1] = "d/dt Pin_f in component gas_exchange (mmHg)"
    return (legend_states, legend_algebraic, legend_voi, legend_constants)

def initConsts():
    constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
    constants[0] = 1.4452
    constants[1] = 1.4452
    constants[2] = 0.456
    constants[3] = 0.3711
    constants[4] = 0.3711
    constants[5] = 7.8322e+03
    constants[6] = 7.4
    constants[7] = 1.302
    constants[8] = 0.464
    constants[9] = 0.395
    constants[10] = 2.2643e+03
    constants[11] = 7.35
    states[0] = 33
    states[1] = 15
    constants[12] = 12.5
    constants[13] = 15.5
    constants[14] = 3.0263e-05
    constants[15] = 164
    constants[16] = 164
    constants[17] = 164
    constants[18] = 1
    constants[19] = 4.89e-7
    constants[20] = 4.89e-7
    constants[21] = 1.2014e-04
    constants[22] = 1.2014e-04
    constants[23] = 6.21
    constants[24] = 3.287
    constants[25] = 0.1117
    constants[26] = 7.05e-3
    constants[27] = 0.8142
    constants[28] = 1.00000/constants[23]
    constants[29] = constants[20]/(constants[22]/60.0000)
    constants[30] = constants[29]/100.000
    constants[31] = 1.00000/constants[23]
    constants[32] = 1.00000/constants[23]
    constants[33] = 1.00000/constants[23]
    return (states, constants)

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

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