# Size of variable arrays: sizeAlgebraic = 5 sizeStates = 48 sizeConstants = 101 from math import * from numpy import * def createLegends(): legend_states = [""] * sizeStates legend_rates = [""] * sizeStates legend_algebraic = [""] * sizeAlgebraic legend_voi = "" legend_constants = [""] * sizeConstants legend_voi = "time in component environment (day)" legend_algebraic[0] = "FAt in component FAt (litre_per_day)" legend_constants[0] = "FA1 in component FAt (litre_per_day)" legend_constants[1] = "tA in component FAt (day)" legend_constants[2] = "FA2 in component FAt (litre_per_day)" legend_algebraic[3] = "FBt in component FBt (litre_per_day)" legend_constants[3] = "FB1 in component FBt (litre_per_day)" legend_constants[4] = "tB in component FBt (day)" legend_constants[5] = "FB2 in component FBt (litre_per_day)" legend_algebraic[4] = "F0t in component F0t (litre_per_day)" legend_states[0] = "c15 in component c15 (micromole_per_litre)" legend_constants[6] = "p0 in component model_parameters (litre_per_day)" legend_constants[7] = "V2 in component model_parameters (litre)" legend_states[1] = "c015 in component c015 (micromole_per_litre)" legend_constants[8] = "p16 in component model_parameters (per_day)" legend_states[2] = "c16 in component c16 (micromole_per_litre)" legend_constants[9] = "p18 in component model_parameters (per_day)" legend_states[3] = "c17 in component c17 (micromole_per_litre)" legend_constants[10] = "p15 in component model_parameters (per_day_per_micromole_litre)" legend_states[4] = "c19 in component c19 (micromole_per_litre)" legend_constants[11] = "p17 in component model_parameters (per_day)" legend_constants[12] = "p19 in component model_parameters (per_day_per_micromole_litre)" legend_constants[13] = "p20 in component model_parameters (micromole_per_litre)" legend_constants[14] = "p29 in component model_parameters (per_day)" legend_constants[15] = "pAA15 in component model_parameters (dimensionless)" legend_algebraic[1] = "pp in component pp (dimensionless)" legend_states[5] = "c016 in component c016 (micromole_per_litre)" legend_constants[16] = "pAA16 in component model_parameters (dimensionless)" legend_states[6] = "c017 in component c017 (micromole_per_litre)" legend_constants[17] = "p22 in component model_parameters (per_day)" legend_states[7] = "c18 in component c18 (micromole_per_litre)" legend_constants[18] = "p21 in component model_parameters (per_day)" legend_constants[19] = "p23 in component model_parameters (per_day)" legend_constants[20] = "p1 in component model_parameters (per_day)" legend_constants[21] = "s1 in component c17 (dimensionless)" legend_states[8] = "c1 in component c1 (micromole_per_litre)" legend_constants[22] = "p2 in component model_parameters (per_day)" legend_constants[23] = "s2 in component c17 (dimensionless)" legend_states[9] = "c2 in component c2 (micromole_per_litre)" legend_constants[24] = "p3 in component model_parameters (per_day)" legend_constants[25] = "s3 in component c17 (dimensionless)" legend_states[10] = "c3 in component c3 (micromole_per_litre)" legend_constants[26] = "p4 in component model_parameters (per_day)" legend_constants[27] = "s4 in component c17 (dimensionless)" legend_states[11] = "c4 in component c4 (micromole_per_litre)" legend_constants[28] = "p5 in component model_parameters (per_day)" legend_constants[29] = "s5 in component c17 (dimensionless)" legend_states[12] = "c5 in component c5 (micromole_per_litre)" legend_constants[30] = "p6 in component model_parameters (per_day)" legend_constants[31] = "s6 in component c17 (dimensionless)" legend_states[13] = "c6 in component c6 (micromole_per_litre)" legend_constants[32] = "p7 in component model_parameters (per_day)" legend_constants[33] = "s7 in component c17 (dimensionless)" legend_states[14] = "c7 in component c7 (micromole_per_litre)" legend_constants[34] = "p8 in component model_parameters (per_day)" legend_constants[35] = "s8 in component c17 (dimensionless)" legend_states[15] = "c8 in component c8 (micromole_per_litre)" legend_constants[36] = "p9 in component model_parameters (per_day)" legend_constants[37] = "s9 in component c17 (dimensionless)" legend_states[16] = "c9 in component c9 (micromole_per_litre)" legend_constants[38] = "p14 in component model_parameters (per_day)" legend_constants[39] = "s14 in component c17 (dimensionless)" legend_states[17] = "c14 in component c14 (micromole_per_litre)" legend_constants[40] = "p25 in component model_parameters (per_day_per_micromole_litre)" legend_constants[41] = "p26 in component model_parameters (per_day)" legend_constants[42] = "p28 in component model_parameters (per_day)" legend_constants[43] = "pAA17 in component model_parameters (dimensionless)" legend_states[18] = "c018 in component c018 (micromole_per_litre)" legend_constants[44] = "pAA18 in component model_parameters (dimensionless)" legend_states[19] = "c019 in component c019 (micromole_per_litre)" legend_algebraic[2] = "gt in component gt (dimensionless)" legend_constants[45] = "p11 in component model_parameters (per_day)" legend_constants[46] = "s11 in component c19 (dimensionless)" legend_states[20] = "c11 in component c11 (micromole_per_litre)" legend_constants[47] = "p12 in component model_parameters (per_day)" legend_constants[48] = "s12 in component c19 (dimensionless)" legend_states[21] = "c12 in component c12 (micromole_per_litre)" legend_constants[49] = "p13 in component model_parameters (per_day)" legend_constants[50] = "s13 in component c19 (dimensionless)" legend_states[22] = "c13 in component c13 (micromole_per_litre)" legend_constants[51] = "p10 in component model_parameters (per_day)" legend_constants[52] = "s10 in component c19 (dimensionless)" legend_states[23] = "c10 in component c10 (micromole_per_litre)" legend_states[24] = "c20 in component c20 (micromole_per_litre)" legend_states[25] = "c020 in component c020 (micromole_per_litre)" legend_states[26] = "c21 in component c21 (micromole_per_litre)" legend_states[27] = "c021 in component c021 (micromole_per_litre)" legend_states[28] = "c22 in component c22 (micromole_per_litre)" legend_states[29] = "c022 in component c022 (micromole_per_litre)" legend_states[30] = "c23 in component c23 (micromole_per_litre)" legend_states[31] = "c023 in component c023 (micromole_per_litre)" legend_constants[53] = "p31 in component model_parameters (per_day)" legend_states[32] = "c24 in component c24 (micromole_per_litre)" legend_states[33] = "c024 in component c024 (micromole_per_litre)" legend_constants[54] = "p30 in component model_parameters (dimensionless)" legend_constants[55] = "p32 in component model_parameters (per_day)" legend_constants[56] = "p24 in component model_parameters (per_day)" legend_constants[57] = "p33 in component pp (micromole_per_litre_per_day)" legend_constants[58] = "t1 in component pp (day)" legend_constants[59] = "p27 in component model_parameters (dimensionless)" legend_states[34] = "c01 in component c01 (micromole_per_litre)" legend_constants[60] = "pAA1 in component model_parameters (dimensionless)" legend_states[35] = "c02 in component c02 (micromole_per_litre)" legend_constants[61] = "pAA2 in component model_parameters (dimensionless)" legend_states[36] = "c03 in component c03 (micromole_per_litre)" legend_constants[62] = "pAA3 in component model_parameters (dimensionless)" legend_states[37] = "c04 in component c04 (micromole_per_litre)" legend_constants[63] = "pAA4 in component model_parameters (dimensionless)" legend_states[38] = "c05 in component c05 (micromole_per_litre)" legend_constants[64] = "pAA5 in component model_parameters (dimensionless)" legend_states[39] = "c06 in component c06 (micromole_per_litre)" legend_constants[65] = "pAA6 in component model_parameters (dimensionless)" legend_states[40] = "c07 in component c07 (micromole_per_litre)" legend_constants[66] = "pAA7 in component model_parameters (dimensionless)" legend_states[41] = "c08 in component c08 (micromole_per_litre)" legend_constants[67] = "pAA8 in component model_parameters (dimensionless)" legend_states[42] = "c09 in component c09 (micromole_per_litre)" legend_constants[68] = "pAA9 in component model_parameters (dimensionless)" legend_states[43] = "c010 in component c010 (micromole_per_litre)" legend_constants[69] = "pAA10 in component model_parameters (dimensionless)" legend_states[44] = "c011 in component c011 (micromole_per_litre)" legend_constants[70] = "pAA11 in component model_parameters (dimensionless)" legend_states[45] = "c012 in component c012 (micromole_per_litre)" legend_constants[71] = "pAA12 in component model_parameters (dimensionless)" legend_states[46] = "c013 in component c013 (micromole_per_litre)" legend_constants[72] = "pAA13 in component model_parameters (dimensionless)" legend_states[47] = "c014 in component c014 (micromole_per_litre)" legend_constants[73] = "pAA14 in component model_parameters (dimensionless)" legend_constants[74] = "V1 in component model_parameters (litre)" legend_constants[75] = "CA1 in component c01 (micromole_per_litre)" legend_constants[76] = "CB1 in component model_parameters (micromole_per_litre)" legend_constants[77] = "CA2 in component c02 (micromole_per_litre)" legend_constants[78] = "CA3 in component c03 (micromole_per_litre)" legend_constants[79] = "CA4 in component c04 (micromole_per_litre)" legend_constants[80] = "CA5 in component c05 (micromole_per_litre)" legend_constants[81] = "CA6 in component c06 (micromole_per_litre)" legend_constants[82] = "CA7 in component c07 (micromole_per_litre)" legend_constants[83] = "CA8 in component c08 (micromole_per_litre)" legend_constants[84] = "CA9 in component c09 (micromole_per_litre)" legend_constants[85] = "CA10 in component c010 (micromole_per_litre)" legend_constants[86] = "CA11 in component c011 (micromole_per_litre)" legend_constants[87] = "CA12 in component c012 (micromole_per_litre)" legend_constants[88] = "CA13 in component c013 (micromole_per_litre)" legend_constants[89] = "CA14 in component c014 (micromole_per_litre)" legend_constants[90] = "CA15 in component c015 (micromole_per_litre)" legend_constants[91] = "CB15 in component model_parameters (micromole_per_litre)" legend_constants[92] = "CA16 in component c016 (micromole_per_litre)" legend_constants[93] = "CA17 in component c017 (micromole_per_litre)" legend_constants[94] = "CA18 in component c018 (micromole_per_litre)" legend_constants[95] = "CA19 in component c019 (micromole_per_litre)" legend_constants[96] = "CA20 in component c020 (micromole_per_litre)" legend_constants[97] = "CA21 in component c021 (micromole_per_litre)" legend_constants[98] = "CA22 in component c022 (micromole_per_litre)" legend_constants[99] = "CA23 in component c023 (micromole_per_litre)" legend_constants[100] = "CA24 in component c024 (micromole_per_litre)" legend_rates[0] = "d/dt c15 in component c15 (micromole_per_litre)" legend_rates[2] = "d/dt c16 in component c16 (micromole_per_litre)" legend_rates[3] = "d/dt c17 in component c17 (micromole_per_litre)" legend_rates[7] = "d/dt c18 in component c18 (micromole_per_litre)" legend_rates[4] = "d/dt c19 in component c19 (micromole_per_litre)" legend_rates[24] = "d/dt c20 in component c20 (micromole_per_litre)" legend_rates[26] = "d/dt c21 in component c21 (micromole_per_litre)" legend_rates[28] = "d/dt c22 in component c22 (micromole_per_litre)" legend_rates[30] = "d/dt c23 in component c23 (micromole_per_litre)" legend_rates[32] = "d/dt c24 in component c24 (micromole_per_litre)" legend_rates[8] = "d/dt c1 in component c1 (micromole_per_litre)" legend_rates[9] = "d/dt c2 in component c2 (micromole_per_litre)" legend_rates[10] = "d/dt c3 in component c3 (micromole_per_litre)" legend_rates[11] = "d/dt c4 in component c4 (micromole_per_litre)" legend_rates[12] = "d/dt c5 in component c5 (micromole_per_litre)" legend_rates[13] = "d/dt c6 in component c6 (micromole_per_litre)" legend_rates[14] = "d/dt c7 in component c7 (micromole_per_litre)" legend_rates[15] = "d/dt c8 in component c8 (micromole_per_litre)" legend_rates[16] = "d/dt c9 in component c9 (micromole_per_litre)" legend_rates[23] = "d/dt c10 in component c10 (micromole_per_litre)" legend_rates[20] = "d/dt c11 in component c11 (micromole_per_litre)" legend_rates[21] = "d/dt c12 in component c12 (micromole_per_litre)" legend_rates[22] = "d/dt c13 in component c13 (micromole_per_litre)" legend_rates[17] = "d/dt c14 in component c14 (micromole_per_litre)" legend_rates[34] = "d/dt c01 in component c01 (micromole_per_litre)" legend_rates[35] = "d/dt c02 in component c02 (micromole_per_litre)" legend_rates[36] = "d/dt c03 in component c03 (micromole_per_litre)" legend_rates[37] = "d/dt c04 in component c04 (micromole_per_litre)" legend_rates[38] = "d/dt c05 in component c05 (micromole_per_litre)" legend_rates[39] = "d/dt c06 in component c06 (micromole_per_litre)" legend_rates[40] = "d/dt c07 in component c07 (micromole_per_litre)" legend_rates[41] = "d/dt c08 in component c08 (micromole_per_litre)" legend_rates[42] = "d/dt c09 in component c09 (micromole_per_litre)" legend_rates[43] = "d/dt c010 in component c010 (micromole_per_litre)" legend_rates[44] = "d/dt c011 in component c011 (micromole_per_litre)" legend_rates[45] = "d/dt c012 in component c012 (micromole_per_litre)" legend_rates[46] = "d/dt c013 in component c013 (micromole_per_litre)" legend_rates[47] = "d/dt c014 in component c014 (micromole_per_litre)" legend_rates[1] = "d/dt c015 in component c015 (micromole_per_litre)" legend_rates[5] = "d/dt c016 in component c016 (micromole_per_litre)" legend_rates[6] = "d/dt c017 in component c017 (micromole_per_litre)" legend_rates[18] = "d/dt c018 in component c018 (micromole_per_litre)" legend_rates[19] = "d/dt c019 in component c019 (micromole_per_litre)" legend_rates[25] = "d/dt c020 in component c020 (micromole_per_litre)" legend_rates[27] = "d/dt c021 in component c021 (micromole_per_litre)" legend_rates[29] = "d/dt c022 in component c022 (micromole_per_litre)" legend_rates[31] = "d/dt c023 in component c023 (micromole_per_litre)" legend_rates[33] = "d/dt c024 in component c024 (micromole_per_litre)" return (legend_states, legend_algebraic, legend_voi, legend_constants) def initConsts(): constants = [0.0] * sizeConstants; states = [0.0] * sizeStates; constants[0] = 3.6 constants[1] = 1 constants[2] = 1.2 constants[3] = 0 constants[4] = 3 constants[5] = 0.024 states[0] = 284 constants[6] = 1191.072 constants[7] = 0.9 states[1] = 284 constants[8] = 3.593 states[2] = 168 constants[9] = 56.73 states[3] = 265 constants[10] = 0.0009 states[4] = 41 constants[11] = 50.08 constants[12] = 225.76 constants[13] = 5 constants[14] = 2.8119 constants[15] = 0.31 states[5] = 168 constants[16] = 0.25 states[6] = 265 constants[17] = 0.0205 states[7] = 687 constants[18] = 0.0033 constants[19] = 40.68 constants[20] = 1.816 constants[21] = 1 states[8] = 2258 constants[22] = 16.29 constants[23] = 3 states[9] = 940 constants[24] = 100 constants[25] = 4 states[10] = 2604 constants[26] = 1.147 constants[27] = 1 states[11] = 987 constants[28] = 15.99 constants[29] = 2 states[12] = 464 constants[30] = 4.148 constants[31] = 1 states[13] = 1300 constants[32] = 3.723 constants[33] = 1 states[14] = 500 constants[34] = 1.573 constants[35] = 1 states[15] = 1825 constants[36] = 1.135 constants[37] = 1 states[16] = 2215 constants[38] = 9.326 constants[39] = 1 states[17] = 859 constants[40] = 0 constants[41] = 0 constants[42] = 38.05 constants[43] = 1.27 states[18] = 687 constants[44] = 0.06 states[19] = 41 constants[45] = 34.53 constants[46] = 1 states[20] = 259 constants[47] = 10.8 constants[48] = 1 states[21] = 861 constants[49] = 10.76 constants[50] = 1 states[22] = 1957 constants[51] = 0.192 constants[52] = 1 states[23] = 327 states[24] = 0 states[25] = 0 states[26] = 0 states[27] = 0 states[28] = 0 states[29] = 0 states[30] = 8159 states[31] = 0 constants[53] = 1.31 states[32] = 8159 states[33] = 0 constants[54] = 0 constants[55] = 478.7 constants[56] = 0.5474 constants[57] = 0 constants[58] = 1 constants[59] = 0.6804 states[34] = 2258 constants[60] = 1.4 states[35] = 940 constants[61] = 0.18 states[36] = 2604 constants[62] = 0 states[37] = 987 constants[63] = 1.19 states[38] = 464 constants[64] = 0.45 states[39] = 1300 constants[65] = 1.05 states[40] = 500 constants[66] = 0.62 states[41] = 1825 constants[67] = 1.15 states[42] = 2215 constants[68] = 1 states[43] = 327 constants[69] = 0.15 states[44] = 259 constants[70] = 0.77 states[45] = 861 constants[71] = 1.55 states[46] = 1957 constants[72] = 0.78 states[47] = 859 constants[73] = 1.81 constants[74] = 0.6 constants[75] = 2258 constants[76] = 0 constants[77] = 940 constants[78] = 2604 constants[79] = 987 constants[80] = 464 constants[81] = 1300 constants[82] = 500 constants[83] = 1825 constants[84] = 2215 constants[85] = 327 constants[86] = 259 constants[87] = 861 constants[88] = 1957 constants[89] = 859 constants[90] = 284 constants[91] = 1500 constants[92] = 168 constants[93] = 265 constants[94] = 687 constants[95] = 41 constants[96] = 0 constants[97] = 0 constants[98] = 0 constants[99] = 0 constants[100] = 0 return (states, constants) def computeRates(voi, states, constants): rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic rates[24] = (constants[6]/constants[7])*(states[25]-states[24])+((constants[12]*states[0])/(constants[13]+states[0]))*states[4] rates[26] = (constants[6]/constants[7])*(states[27]-states[26])-constants[18]*states[26] rates[28] = (constants[6]/constants[7])*(states[29]-states[28])-constants[17]*states[28] rates[30] = (((constants[6]/constants[7])*(states[31]-states[30])+constants[18]*states[26]+constants[17]*states[28])-constants[19]*states[30])+((constants[53]*states[0])/(constants[13]+states[0]))*states[4] rates[32] = ((constants[6]/constants[7])*(states[33]-states[32])+2.00000*constants[54]*constants[19]*states[30]+((constants[55]*states[0])/(constants[13]+states[0]))*states[4])-constants[56]*states[32] algebraic[1] = custom_piecewise([less(voi , constants[58]), constants[57] , True, 0.00000]) rates[0] = ((constants[6]/constants[7])*(states[1]-states[0])+constants[8]*states[2]+constants[11]*states[3]+constants[15]*algebraic[1])-(constants[10]*states[4]+constants[9]+(constants[12]*states[4])/(constants[13]+states[0])+constants[14])*states[0] rates[2] = (((constants[6]/constants[7])*(states[5]-states[2])+constants[10]*states[0]*states[4])-constants[8]*states[2])+constants[16]*algebraic[1] rates[3] = (((constants[6]/constants[7])*(states[6]-states[3])+constants[20]*constants[21]*states[8]+constants[22]*constants[23]*states[9]+constants[24]*constants[25]*states[10]+constants[26]*constants[27]*states[11]+constants[28]*constants[29]*states[12]+constants[30]*constants[31]*states[13]+constants[32]*constants[33]*states[14]+constants[34]*constants[35]*states[15]+constants[36]*constants[37]*states[16]+constants[38]*constants[39]*states[17]+constants[9]*states[0]+constants[41]*states[7])-(constants[11]+constants[40]*states[4]+constants[42])*states[3])+constants[43]*algebraic[1] rates[7] = (((constants[6]/constants[7])*(states[18]-states[7])+constants[40]*states[3]*states[4])-constants[41]*states[7])+constants[44]*algebraic[1] algebraic[2] = custom_piecewise([less(voi , 3.00000), 0.00000 , True, constants[59]]) rates[4] = ((constants[6]/constants[7])*(states[19]-states[4])+constants[45]*constants[46]*states[20]+constants[47]*constants[48]*states[21]+constants[49]*states[22]*constants[50]+constants[51]*states[23]*constants[52]+constants[19]*(1.00000-algebraic[2])*states[3]+constants[8]*states[2]+constants[41]*states[7])-((constants[12]*states[0])/(constants[13]+states[0])+constants[10]*states[0]+constants[40]*states[3])*states[4] rates[8] = ((constants[6]/constants[7])*(states[34]-states[8])-constants[20]*states[8])+constants[60]*algebraic[1] rates[9] = ((constants[6]/constants[7])*(states[35]-states[9])-constants[22]*states[9])+constants[61]*algebraic[1] rates[10] = ((constants[6]/constants[7])*(states[36]-states[10])-constants[24]*states[10])+constants[62]*algebraic[1] rates[11] = ((constants[6]/constants[7])*(states[37]-states[11])-constants[26]*states[11])+constants[63]*algebraic[1] rates[12] = ((constants[6]/constants[7])*(states[38]-states[12])-constants[28]*states[12])+constants[64]*algebraic[1] rates[13] = ((constants[6]/constants[7])*(states[39]-states[13])-constants[30]*states[13])+constants[65]*algebraic[1] rates[14] = ((constants[6]/constants[7])*(states[40]-states[14])-constants[32]*states[14])+constants[66]*algebraic[1] rates[15] = ((constants[6]/constants[7])*(states[41]-states[15])-constants[34]*states[15])+constants[67]*algebraic[1] rates[16] = ((constants[6]/constants[7])*(states[42]-states[16])-constants[36]*states[16])+constants[68]*algebraic[1] rates[23] = ((constants[6]/constants[7])*(states[43]-states[23])-constants[51]*states[23])+constants[69]*algebraic[1] rates[20] = ((constants[6]/constants[7])*(states[44]-states[20])-constants[45]*states[20])+constants[70]*algebraic[1] rates[21] = ((constants[6]/constants[7])*(states[45]-states[21])-constants[47]*states[21])+constants[71]*algebraic[1] rates[22] = ((constants[6]/constants[7])*(states[46]-states[22])-constants[49]*states[22])+constants[72]*algebraic[1] rates[17] = ((constants[6]/constants[7])*(states[47]-states[17])-constants[38]*states[17])+constants[73]*algebraic[1] algebraic[0] = custom_piecewise([greater_equal(voi , constants[1]), constants[2] , True, constants[0]]) algebraic[3] = custom_piecewise([greater_equal(voi , constants[4]), constants[5] , True, constants[3]]) algebraic[4] = algebraic[0]+algebraic[3] rates[34] = (((algebraic[0]/constants[74])*constants[75]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[34])-(constants[6]/constants[74])*(states[34]-states[8]) rates[35] = (((algebraic[0]/constants[74])*constants[77]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[35])-(constants[6]/constants[74])*(states[35]-states[9]) rates[36] = (((algebraic[0]/constants[74])*constants[78]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[36])-(constants[6]/constants[74])*(states[36]-states[10]) rates[37] = (((algebraic[0]/constants[74])*constants[79]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[37])-(constants[6]/constants[74])*(states[37]-states[11]) rates[38] = (((algebraic[0]/constants[74])*constants[80]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[38])-(constants[6]/constants[74])*(states[38]-states[12]) rates[39] = (((algebraic[0]/constants[74])*constants[81]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[39])-(constants[6]/constants[74])*(states[39]-states[13]) rates[40] = (((algebraic[0]/constants[74])*constants[82]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[40])-(constants[6]/constants[74])*(states[40]-states[14]) rates[41] = (((algebraic[0]/constants[74])*constants[83]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[41])-(constants[6]/constants[74])*(states[41]-states[15]) rates[42] = (((algebraic[0]/constants[74])*constants[84]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[42])-(constants[6]/constants[74])*(states[42]-states[16]) rates[43] = (((algebraic[0]/constants[74])*constants[85]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[43])-(constants[6]/constants[74])*(states[43]-states[23]) rates[44] = (((algebraic[0]/constants[74])*constants[86]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[44])-(constants[6]/constants[74])*(states[44]-states[20]) rates[45] = (((algebraic[0]/constants[74])*constants[87]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[45])-(constants[6]/constants[74])*(states[45]-states[21]) rates[46] = (((algebraic[0]/constants[74])*constants[88]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[46])-(constants[6]/constants[74])*(states[46]-states[22]) rates[47] = (((algebraic[0]/constants[74])*constants[89]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[47])-(constants[6]/constants[74])*(states[47]-states[17]) rates[1] = (((algebraic[0]/constants[74])*constants[90]+(algebraic[3]/constants[74])*constants[91])-(algebraic[4]/constants[74])*states[1])-(constants[6]/constants[74])*(states[1]-states[0]) rates[5] = (((algebraic[0]/constants[74])*constants[92]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[5])-(constants[6]/constants[74])*(states[5]-states[2]) rates[6] = (((algebraic[0]/constants[74])*constants[93]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[6])-(constants[6]/constants[74])*(states[6]-states[3]) rates[18] = (((algebraic[0]/constants[74])*constants[94]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[18])-(constants[6]/constants[74])*(states[18]-states[7]) rates[19] = (((algebraic[0]/constants[74])*constants[95]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[19])-(constants[6]/constants[74])*(states[19]-states[4]) rates[25] = (((algebraic[0]/constants[74])*constants[96]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[25])-(constants[6]/constants[74])*(states[25]-states[24]) rates[27] = (((algebraic[0]/constants[74])*constants[97]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[27])-(constants[6]/constants[74])*(states[27]-states[26]) rates[29] = (((algebraic[0]/constants[74])*constants[98]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[29])-(constants[6]/constants[74])*(states[29]-states[28]) rates[31] = (((algebraic[0]/constants[74])*constants[99]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[31])-(constants[6]/constants[74])*(states[31]-states[30]) rates[33] = (((algebraic[0]/constants[74])*constants[100]+(algebraic[3]/constants[74])*constants[76])-(algebraic[4]/constants[74])*states[33])-(constants[6]/constants[74])*(states[33]-states[32]) return(rates) def computeAlgebraic(constants, states, voi): algebraic = array([[0.0] * len(voi)] * sizeAlgebraic) states = array(states) voi = array(voi) algebraic[1] = custom_piecewise([less(voi , constants[58]), constants[57] , True, 0.00000]) algebraic[2] = custom_piecewise([less(voi , 3.00000), 0.00000 , True, constants[59]]) algebraic[0] = custom_piecewise([greater_equal(voi , constants[1]), constants[2] , True, constants[0]]) algebraic[3] = custom_piecewise([greater_equal(voi , constants[4]), constants[5] , True, constants[3]]) algebraic[4] = algebraic[0]+algebraic[3] 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)