Generated Code
The following is python code generated by the CellML API from this CellML file. (Back to language selection)
The raw code is available.
# 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)