# Size of variable arrays: sizeAlgebraic = 86 sizeStates = 44 sizeConstants = 96 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 Circ_Environment (ms)" legend_constants[0] = "HR in component Circ_Environment (ratepm)" legend_constants[91] = "hrf in component Circ_Environment (Hz)" legend_algebraic[0] = "beattime in component LATiming (ms)" legend_constants[1] = "Esys in component LAElastanceFunction (elastance)" legend_constants[2] = "Edia in component LAElastanceFunction (elastance)" legend_constants[3] = "TsK in component LAElastanceFunction (ms)" legend_constants[92] = "Ts in component LAElastanceFunction (ms)" legend_algebraic[1] = "E_LA in component LAElastanceFunction (elastance)" legend_algebraic[82] = "P_PV in component PulmonaryVein (mmHg)" legend_algebraic[9] = "F_LV in component LeftVentricle (flow)" legend_constants[4] = "P_LA_ext in component LeftAtrium (mmHg)" legend_states[0] = "V_LA in component LeftAtrium (mL)" legend_constants[5] = "V_LA_0 in component LeftAtrium (mL)" legend_algebraic[3] = "P_LA in component LeftAtrium (mmHg)" legend_algebraic[83] = "F_LA in component LeftAtrium (flow)" legend_constants[6] = "R_LA in component LeftAtrium (resistance)" legend_algebraic[2] = "conc_X_LA in component LeftAtrium (concentration)" legend_algebraic[81] = "conc_X_PV in component PulmonaryVein (concentration)" legend_states[1] = "sub_X_LA in component LeftAtrium (mL)" legend_constants[7] = "PRint in component LVTiming (ms)" legend_algebraic[4] = "beattime in component LVTiming (ms)" legend_constants[8] = "Esys in component LVElastanceFunction (elastance)" legend_constants[9] = "Edia in component LVElastanceFunction (elastance)" legend_constants[10] = "TsK in component LVElastanceFunction (second)" legend_constants[93] = "Ts in component LVElastanceFunction (second)" legend_algebraic[5] = "E_LV in component LVElastanceFunction (elastance)" legend_algebraic[12] = "F_AO in component Aorta (flow)" legend_algebraic[15] = "F_CO in component CoronaryCirc (flow)" legend_constants[11] = "pLVPext in component LeftVentricle (mmHg)" legend_states[2] = "V_LV in component LeftVentricle (mL)" legend_constants[12] = "V_LV_0 in component LeftVentricle (mL)" legend_algebraic[7] = "P_LV in component LeftVentricle (mmHg)" legend_constants[13] = "R_LV in component LeftVentricle (resistance)" legend_algebraic[6] = "conc_X_LV in component LeftVentricle (concentration)" legend_states[3] = "sub_X_LV in component LeftVentricle (mL)" legend_algebraic[8] = "fiber_stress in component LeftVentricle (mmHg)" legend_algebraic[19] = "F_CR in component CarotidCirc (flow)" legend_algebraic[22] = "F_AR in component Arteries (flow)" legend_constants[14] = "C_AO in component Aorta (capacitance)" legend_constants[15] = "L_AO in component Aorta (inductance)" legend_constants[16] = "R_AO in component Aorta (resistance)" legend_constants[17] = "V_AO_0 in component Aorta (mL)" legend_states[4] = "V_AO in component Aorta (mL)" legend_states[5] = "F_AOalways in component Aorta (flow)" legend_algebraic[11] = "P_AO in component Aorta (mmHg)" legend_algebraic[10] = "conc_X_AO in component Aorta (concentration)" legend_states[6] = "sub_X_AO in component Aorta (mL)" legend_algebraic[62] = "P_RA in component RightAtrium (mmHg)" legend_constants[18] = "R_CO in component CoronaryCirc (resistance)" legend_constants[19] = "R_COV in component CoronaryCirc (resistance)" legend_constants[20] = "V_CO_0 in component CoronaryCirc (mL)" legend_constants[21] = "C_CO in component CoronaryCirc (capacitance)" legend_algebraic[63] = "F_COV in component CoronaryCirc (flow)" legend_states[7] = "V_CO in component CoronaryCirc (mL)" legend_algebraic[14] = "P_CO in component CoronaryCirc (mmHg)" legend_algebraic[13] = "conc_X_CO in component CoronaryCirc (concentration)" legend_states[8] = "sub_X_CO in component CoronaryCirc (mL)" legend_algebraic[55] = "P_VC in component VenaCava (mmHg)" legend_constants[22] = "R_CR in component CarotidCirc (resistance)" legend_constants[23] = "R_CRV in component CarotidCirc (resistance)" legend_constants[24] = "V_CR_0 in component CarotidCirc (mL)" legend_constants[25] = "C_CR in component CarotidCirc (capacitance)" legend_algebraic[56] = "F_CRV in component CarotidCirc (flow)" legend_states[9] = "V_CR in component CarotidCirc (mL)" legend_algebraic[18] = "P_CR in component CarotidCirc (mmHg)" legend_algebraic[17] = "conc_X_CR in component CarotidCirc (concentration)" legend_states[10] = "sub_X_CR in component CarotidCirc (mL)" legend_algebraic[16] = "conc_X_CRtissue in component CarotidCirc (concentration)" legend_constants[26] = "D_CRtissue in component CarotidCirc (flow)" legend_constants[27] = "V_CRtissue in component CarotidCirc (mL)" legend_states[11] = "sub_X_CRtissue in component CarotidCirc (mL)" legend_algebraic[25] = "F_AD in component Adipose (flow)" legend_algebraic[28] = "F_MU in component Muscle (flow)" legend_algebraic[31] = "F_GI in component GutIntestine (flow)" legend_algebraic[35] = "F_LI in component Liver (flow)" legend_algebraic[38] = "F_KI in component Kidney (flow)" legend_algebraic[44] = "F_OT in component OtherTissue (flow)" legend_algebraic[41] = "F_SK in component Skin (flow)" legend_constants[28] = "R_AR in component Arteries (resistance)" legend_constants[29] = "C_AR in component Arteries (capacitance)" legend_constants[30] = "V_AR_0 in component Arteries (mL)" legend_states[12] = "V_AR in component Arteries (mL)" legend_algebraic[21] = "P_AR in component Arteries (mmHg)" legend_algebraic[20] = "conc_X_AR in component Arteries (concentration)" legend_states[13] = "sub_X_AR in component Arteries (mL)" legend_algebraic[46] = "P_VE in component Veins (mmHg)" legend_constants[31] = "C_AD in component Adipose (capacitance)" legend_constants[32] = "R_AD in component Adipose (resistance)" legend_constants[33] = "R_ADV in component Adipose (resistance)" legend_constants[34] = "V_AD_0 in component Adipose (mL)" legend_states[14] = "V_AD in component Adipose (mL)" legend_algebraic[47] = "F_ADV in component Adipose (flow)" legend_algebraic[24] = "P_AD in component Adipose (mmHg)" legend_algebraic[23] = "conc_X_AD in component Adipose (concentration)" legend_states[15] = "sub_X_AD in component Adipose (mL)" legend_constants[35] = "C_MU in component Muscle (capacitance)" legend_constants[36] = "R_MU in component Muscle (resistance)" legend_constants[37] = "R_MUV in component Muscle (resistance)" legend_constants[38] = "V_MU_0 in component Muscle (mL)" legend_states[16] = "V_MU in component Muscle (mL)" legend_algebraic[48] = "F_MUV in component Muscle (flow)" legend_algebraic[27] = "P_MU in component Muscle (mmHg)" legend_algebraic[26] = "conc_X_MU in component Muscle (concentration)" legend_states[17] = "sub_X_MU in component Muscle (mL)" legend_algebraic[33] = "P_LI in component Liver (mmHg)" legend_constants[39] = "C_GI in component GutIntestine (capacitance)" legend_constants[40] = "R_GI in component GutIntestine (resistance)" legend_constants[41] = "R_GIV in component GutIntestine (resistance)" legend_constants[42] = "V_GI_0 in component GutIntestine (mL)" legend_states[18] = "V_GI in component GutIntestine (mL)" legend_algebraic[34] = "F_GIV in component GutIntestine (flow)" legend_algebraic[30] = "P_GI in component GutIntestine (mmHg)" legend_algebraic[29] = "conc_X_GI in component GutIntestine (concentration)" legend_states[19] = "sub_X_GI in component GutIntestine (mL)" legend_constants[43] = "C_LI in component Liver (capacitance)" legend_constants[44] = "R_LI in component Liver (resistance)" legend_constants[45] = "R_LIV in component Liver (resistance)" legend_constants[46] = "V_LI_0 in component Liver (mL)" legend_states[20] = "V_LI in component Liver (mL)" legend_algebraic[49] = "F_LIV in component Liver (flow)" legend_algebraic[32] = "conc_X_LI in component Liver (concentration)" legend_states[21] = "sub_X_LI in component Liver (mL)" legend_constants[47] = "metabolism in component Liver (flow)" legend_constants[48] = "C_KI in component Kidney (capacitance)" legend_constants[49] = "R_KI in component Kidney (resistance)" legend_constants[50] = "V_KI_0 in component Kidney (mL)" legend_constants[51] = "R_KIV in component Kidney (resistance)" legend_states[22] = "V_KI in component Kidney (mL)" legend_algebraic[50] = "F_KIV in component Kidney (flow)" legend_algebraic[37] = "P_KI in component Kidney (mmHg)" legend_algebraic[36] = "conc_X_KI in component Kidney (concentration)" legend_states[23] = "sub_X_KI in component Kidney (mL)" legend_constants[52] = "C_SK in component Skin (capacitance)" legend_constants[53] = "R_SK in component Skin (resistance)" legend_constants[54] = "R_SKV in component Skin (resistance)" legend_constants[55] = "V_SK_0 in component Skin (mL)" legend_states[24] = "V_SK in component Skin (mL)" legend_algebraic[51] = "F_SKV in component Skin (flow)" legend_algebraic[40] = "P_SK in component Skin (mmHg)" legend_algebraic[39] = "conc_X_SK in component Skin (concentration)" legend_states[25] = "sub_X_SK in component Skin (mL)" legend_constants[56] = "C_OT in component OtherTissue (capacitance)" legend_constants[57] = "R_OT in component OtherTissue (resistance)" legend_constants[58] = "R_OTV in component OtherTissue (resistance)" legend_constants[59] = "V_OT_0 in component OtherTissue (mL)" legend_states[26] = "V_OT in component OtherTissue (mL)" legend_algebraic[52] = "F_OTV in component OtherTissue (flow)" legend_algebraic[43] = "P_OT in component OtherTissue (mmHg)" legend_algebraic[42] = "conc_X_OT in component OtherTissue (concentration)" legend_states[27] = "sub_X_OT in component OtherTissue (mL)" legend_constants[60] = "R_VEV in component Veins (resistance)" legend_states[28] = "V_VE in component Veins (mL)" legend_constants[61] = "C_VE in component Veins (capacitance)" legend_constants[62] = "V_VE_0 in component Veins (mL)" legend_algebraic[53] = "F_VE in component Veins (flow)" legend_algebraic[57] = "F_VEV in component Veins (flow)" legend_algebraic[45] = "conc_X_VE in component Veins (concentration)" legend_states[29] = "sub_X_VE in component Veins (mL)" legend_constants[63] = "R_VCV in component VenaCava (resistance)" legend_states[30] = "V_VC in component VenaCava (mL)" legend_constants[64] = "V_VC_0 in component VenaCava (mL)" legend_constants[65] = "C_VC in component VenaCava (capacitance)" legend_algebraic[58] = "F_VC in component VenaCava (flow)" legend_algebraic[64] = "F_VCV in component VenaCava (flow)" legend_algebraic[54] = "conc_X_VC in component VenaCava (concentration)" legend_states[31] = "sub_X_VC in component VenaCava (mL)" legend_algebraic[59] = "beattime in component RATiming (ms)" legend_constants[66] = "Esys in component RAElastanceFunction (elastance)" legend_constants[67] = "Edia in component RAElastanceFunction (elastance)" legend_constants[68] = "TsK in component RAElastanceFunction (second)" legend_constants[94] = "Ts in component RAElastanceFunction (second)" legend_algebraic[60] = "E_RA in component RAElastanceFunction (elastance)" legend_algebraic[70] = "F_RV in component RightVentricle (flow)" legend_constants[69] = "P_RA_ext in component RightAtrium (mmHg)" legend_states[32] = "V_RA in component RightAtrium (mL)" legend_constants[70] = "V_RA_0 in component RightAtrium (mL)" legend_algebraic[65] = "F_RA in component RightAtrium (flow)" legend_algebraic[61] = "conc_X_RA in component RightAtrium (concentration)" legend_states[33] = "sub_X_RA in component RightAtrium (mL)" legend_constants[71] = "PRinterval in component RVTiming (second)" legend_algebraic[66] = "beattime in component RVTiming (ms)" legend_constants[72] = "Esys in component RVElastanceFunction (elastance)" legend_constants[73] = "Edia in component RVElastanceFunction (elastance)" legend_constants[74] = "TsK in component RVElastanceFunction (second)" legend_constants[95] = "Ts in component RVElastanceFunction (second)" legend_algebraic[67] = "E_RV in component RVElastanceFunction (elastance)" legend_algebraic[73] = "F_PA in component PulmonaryArtery (flow)" legend_constants[75] = "P_RV_ext in component RightVentricle (mmHg)" legend_states[34] = "V_RV in component RightVentricle (mL)" legend_constants[76] = "V_RV_0 in component RightVentricle (mL)" legend_algebraic[69] = "P_RV in component RightVentricle (mmHg)" legend_constants[77] = "R_RV in component RightVentricle (resistance)" legend_algebraic[68] = "conc_X_RV in component RightVentricle (concentration)" legend_states[35] = "sub_X_RV in component RightVentricle (mL)" legend_algebraic[76] = "F_Pa in component PulmonaryArteries (flow)" legend_constants[78] = "R_PA in component PulmonaryArtery (resistance)" legend_constants[79] = "C_PA in component PulmonaryArtery (capacitance)" legend_constants[80] = "V_PA_0 in component PulmonaryArtery (mL)" legend_states[36] = "V_PA in component PulmonaryArtery (mL)" legend_algebraic[72] = "P_PA in component PulmonaryArtery (mmHg)" legend_algebraic[71] = "conc_X_PA in component PulmonaryArtery (concentration)" legend_states[37] = "sub_X_PA in component PulmonaryArtery (mL)" legend_algebraic[80] = "F_PC in component PulmonaryCapillaries (flow)" legend_algebraic[79] = "F_SH in component Shunt (flow)" legend_states[38] = "V_Pa in component PulmonaryArteries (mL)" legend_constants[81] = "C_Pa in component PulmonaryArteries (capacitance)" legend_constants[82] = "R_Pa in component PulmonaryArteries (resistance)" legend_constants[83] = "V_Pa_0 in component PulmonaryArteries (mL)" legend_algebraic[75] = "P_Pa in component PulmonaryArteries (mmHg)" legend_algebraic[74] = "conc_X_Pa in component PulmonaryArteries (concentration)" legend_states[39] = "sub_X_Pa in component PulmonaryArteries (mL)" legend_constants[84] = "R_SH in component Shunt (resistance)" legend_algebraic[78] = "P_PC in component PulmonaryCapillaries (mmHg)" legend_constants[85] = "C_PC in component PulmonaryCapillaries (capacitance)" legend_constants[86] = "R_PC in component PulmonaryCapillaries (resistance)" legend_constants[87] = "R_PCV in component PulmonaryCapillaries (resistance)" legend_constants[88] = "V_PC_0 in component PulmonaryCapillaries (mL)" legend_states[40] = "V_PC in component PulmonaryCapillaries (mL)" legend_algebraic[84] = "F_PCV in component PulmonaryCapillaries (flow)" legend_algebraic[77] = "conc_X_PC in component PulmonaryCapillaries (concentration)" legend_states[41] = "sub_X_PC in component PulmonaryCapillaries (mL)" legend_constants[89] = "C_PV in component PulmonaryVein (capacitance)" legend_constants[90] = "V_PV_0 in component PulmonaryVein (mL)" legend_states[42] = "V_PV in component PulmonaryVein (mL)" legend_algebraic[85] = "F_PV in component PulmonaryVein (flow)" legend_states[43] = "sub_X_PV in component PulmonaryVein (mL)" legend_rates[1] = "d/dt sub_X_LA in component LeftAtrium (mL)" legend_rates[0] = "d/dt V_LA in component LeftAtrium (mL)" legend_rates[3] = "d/dt sub_X_LV in component LeftVentricle (mL)" legend_rates[2] = "d/dt V_LV in component LeftVentricle (mL)" legend_rates[6] = "d/dt sub_X_AO in component Aorta (mL)" legend_rates[5] = "d/dt F_AOalways in component Aorta (flow)" legend_rates[4] = "d/dt V_AO in component Aorta (mL)" legend_rates[8] = "d/dt sub_X_CO in component CoronaryCirc (mL)" legend_rates[7] = "d/dt V_CO in component CoronaryCirc (mL)" legend_rates[11] = "d/dt sub_X_CRtissue in component CarotidCirc (mL)" legend_rates[10] = "d/dt sub_X_CR in component CarotidCirc (mL)" legend_rates[9] = "d/dt V_CR in component CarotidCirc (mL)" legend_rates[13] = "d/dt sub_X_AR in component Arteries (mL)" legend_rates[12] = "d/dt V_AR in component Arteries (mL)" legend_rates[15] = "d/dt sub_X_AD in component Adipose (mL)" legend_rates[14] = "d/dt V_AD in component Adipose (mL)" legend_rates[17] = "d/dt sub_X_MU in component Muscle (mL)" legend_rates[16] = "d/dt V_MU in component Muscle (mL)" legend_rates[19] = "d/dt sub_X_GI in component GutIntestine (mL)" legend_rates[18] = "d/dt V_GI in component GutIntestine (mL)" legend_rates[21] = "d/dt sub_X_LI in component Liver (mL)" legend_rates[20] = "d/dt V_LI in component Liver (mL)" legend_rates[23] = "d/dt sub_X_KI in component Kidney (mL)" legend_rates[22] = "d/dt V_KI in component Kidney (mL)" legend_rates[25] = "d/dt sub_X_SK in component Skin (mL)" legend_rates[24] = "d/dt V_SK in component Skin (mL)" legend_rates[27] = "d/dt sub_X_OT in component OtherTissue (mL)" legend_rates[26] = "d/dt V_OT in component OtherTissue (mL)" legend_rates[29] = "d/dt sub_X_VE in component Veins (mL)" legend_rates[28] = "d/dt V_VE in component Veins (mL)" legend_rates[31] = "d/dt sub_X_VC in component VenaCava (mL)" legend_rates[30] = "d/dt V_VC in component VenaCava (mL)" legend_rates[33] = "d/dt sub_X_RA in component RightAtrium (mL)" legend_rates[32] = "d/dt V_RA in component RightAtrium (mL)" legend_rates[35] = "d/dt sub_X_RV in component RightVentricle (mL)" legend_rates[34] = "d/dt V_RV in component RightVentricle (mL)" legend_rates[37] = "d/dt sub_X_PA in component PulmonaryArtery (mL)" legend_rates[36] = "d/dt V_PA in component PulmonaryArtery (mL)" legend_rates[39] = "d/dt sub_X_Pa in component PulmonaryArteries (mL)" legend_rates[38] = "d/dt V_Pa in component PulmonaryArteries (mL)" legend_rates[41] = "d/dt sub_X_PC in component PulmonaryCapillaries (mL)" legend_rates[40] = "d/dt V_PC in component PulmonaryCapillaries (mL)" legend_rates[43] = "d/dt sub_X_PV in component PulmonaryVein (mL)" legend_rates[42] = "d/dt V_PV in component PulmonaryVein (mL)" return (legend_states, legend_algebraic, legend_voi, legend_constants) def initConsts(): constants = [0.0] * sizeConstants; states = [0.0] * sizeStates; constants[0] = 77 constants[1] = 0.145 constants[2] = 0.131 constants[3] = 0.2 constants[4] = 0 states[0] = 144.038898121062 constants[5] = 15 constants[6] = 0.01 states[1] = 0 constants[7] = 0.12 constants[8] = 5.6 constants[9] = 0.19 constants[10] = 0.35 constants[11] = 0 states[2] = 118.75315693072 constants[12] = 30 constants[13] = 0.001 states[3] = 0 constants[14] = 1.2 constants[15] = 5e-5 constants[16] = 0.01 constants[17] = 200 states[4] = 302.714548324312 states[5] = -6898.97886842467 states[6] = 0 constants[18] = 1 constants[19] = 0.2 constants[20] = 10 constants[21] = 0.1 states[7] = 11.4012164759732 states[8] = 0 constants[22] = 5 constants[23] = 0.02 constants[24] = 10 constants[25] = 1 states[9] = 27.2147903053333 states[10] = 0 constants[26] = 1 constants[27] = 1 states[11] = 0 constants[28] = 1.5 constants[29] = 0.5 constants[30] = 370 states[12] = 383.557125288364 states[13] = 0 constants[31] = 1.81 constants[32] = 1 constants[33] = 0.33 constants[34] = 75 states[14] = 111.639088723426 states[15] = 0 constants[35] = 1.81 constants[36] = 1 constants[37] = 0.33 constants[38] = 75 states[16] = 111.639088723426 states[17] = 0 constants[39] = 1.81 constants[40] = 1 constants[41] = 0.33 constants[42] = 75 states[18] = 111.639088723426 states[19] = 0 constants[43] = 1.81 constants[44] = 1 constants[45] = 0.33 constants[46] = 75 states[20] = 111.639088723426 states[21] = 0 constants[47] = 0 constants[48] = 1.81 constants[49] = 1 constants[50] = 75 constants[51] = 0.33 states[22] = 111.639088723426 states[23] = 0 constants[52] = 1.81 constants[53] = 1 constants[54] = 0.33 constants[55] = 75 states[24] = 111.639088723426 states[25] = 0 constants[56] = 1.81 constants[57] = 1 constants[58] = 0.33 constants[59] = 75 states[26] = 111.639088723426 states[27] = 0 constants[60] = 0.0223 states[28] = 833.291684587657 constants[61] = 13.24 constants[62] = 596 states[29] = 5 constants[63] = 0.0267 states[30] = 1849.93788284599 constants[64] = 600 constants[65] = 73.88 states[31] = 0 constants[66] = 0.11 constants[67] = 0.099 constants[68] = 0.2 constants[69] = 0 states[32] = 152.150765432928 constants[70] = 15 states[33] = 0 constants[71] = 0.12 constants[72] = 0.67 constants[73] = 0.104 constants[74] = 0.35 constants[75] = 0 states[34] = 154.337916417774 constants[76] = 25 constants[77] = 0.002 states[35] = 0 constants[78] = 0.0227 constants[79] = 2.222 constants[80] = 50 states[36] = 112.059036699839 states[37] = 0 states[38] = 66.9910537203609 constants[81] = 1.481 constants[82] = 0.053 constants[83] = 30 states[39] = 0 constants[84] = 2 constants[85] = 1.778 constants[86] = 0.0379 constants[87] = 0.0379 constants[88] = 53 states[40] = 91.7784213531534 states[41] = 0 constants[89] = 5 constants[90] = 150 states[42] = 239.25897115602 states[43] = 0 constants[91] = constants[0]/0.0600000 constants[92] = constants[3]*(power(constants[91]/1000.00, 1.0/2)) constants[93] = constants[10]*(power(constants[91]/1000.00, 1.0/2)) constants[94] = constants[68]*(power(constants[91]/1000.00, 1.0/2)) constants[95] = constants[74]*(power(constants[91]/1000.00, 1.0/2)) return (states, constants) def computeRates(voi, states, constants): rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic algebraic[4] = (voi-floor(voi/constants[91])*constants[91])/1000.00-constants[7] algebraic[5] = custom_piecewise([greater_equal(algebraic[4]*0.00100000 , 0.00000) & less_equal(algebraic[4]*0.00100000 , constants[93]), constants[9]+((constants[8]-constants[9])*(1.00000-cos(( pi*(algebraic[4]*0.00100000))/constants[93])))/2.00000 , less(algebraic[4]*0.00100000 , 1.50000*constants[93]) & greater_equal(algebraic[4]*0.00100000 , constants[93]), constants[9]+((constants[8]-constants[9])*(1.00000+cos((2.00000* pi*(algebraic[4]*0.00100000-constants[93]))/constants[93])))/2.00000 , True, constants[9]]) algebraic[7] = algebraic[5]*(states[2]-constants[12])+constants[11] algebraic[11] = (1.00000/constants[14])*(states[4]-constants[17]) rates[5] = ((1.00000/constants[15])*(algebraic[7]-algebraic[11])-(constants[16]/constants[15])*states[5])/1000.00 algebraic[0] = (voi-floor(voi/constants[91])*constants[91])/1000.00 algebraic[1] = custom_piecewise([greater_equal(algebraic[0] , 0.00000) & less_equal(algebraic[0] , constants[92]), constants[2]+((constants[1]-constants[2])*(1.00000-cos(( pi*algebraic[0])/constants[92])))/2.00000 , less(algebraic[0] , 1.50000*constants[92]) & greater_equal(algebraic[0] , constants[92]), constants[2]+((constants[1]-constants[2])*(1.00000+cos((2.00000* pi*(algebraic[0]-constants[92]))/constants[92])))/2.00000 , True, constants[2]]) algebraic[3] = algebraic[1]*(states[0]-constants[5])+constants[4] algebraic[9] = custom_piecewise([greater(algebraic[3] , algebraic[7]), (algebraic[3]-algebraic[7])/constants[13] , True, 0.00000]) algebraic[2] = states[1]/states[0] algebraic[12] = custom_piecewise([greater(algebraic[7] , algebraic[11]), states[5] , True, 0.00000]) algebraic[14] = (1.00000/constants[21])*(states[7]-constants[20]) algebraic[15] = custom_piecewise([greater(algebraic[7] , algebraic[14]), (algebraic[7]-algebraic[14])/constants[18] , True, 0.00000]) algebraic[6] = states[3]/states[2] rates[3] = (algebraic[2]*algebraic[9]-algebraic[6]*(algebraic[12]+algebraic[15]))/1000.00 rates[2] = ((algebraic[9]-algebraic[12])-algebraic[15])/1000.00 algebraic[17] = states[10]/states[9] algebraic[16] = states[11]/constants[27] rates[11] = ((algebraic[17]-algebraic[16])*constants[26])/1000.00 algebraic[18] = (1.00000/constants[25])*(states[9]-constants[24]) algebraic[19] = (algebraic[11]-algebraic[18])/constants[22] algebraic[21] = (1.00000/constants[29])*(states[12]-constants[30]) algebraic[22] = (algebraic[11]-algebraic[21])/constants[28] algebraic[10] = states[6]/states[4] rates[6] = (algebraic[6]*algebraic[12]-algebraic[10]*(algebraic[19]+algebraic[22]))/1000.00 rates[4] = ((algebraic[12]-algebraic[22])-algebraic[19])/1000.00 algebraic[30] = (1.00000/constants[39])*(states[18]-constants[42]) algebraic[31] = (algebraic[21]-algebraic[30])/constants[40] algebraic[20] = states[13]/states[12] algebraic[33] = (1.00000/constants[43])*(states[20]-constants[46]) algebraic[34] = (algebraic[30]-algebraic[33])/constants[41] algebraic[29] = states[19]/states[18] rates[19] = ((algebraic[20]*algebraic[31]-algebraic[29]*algebraic[34])/1000.00)*0.00100000 rates[18] = ((algebraic[31]-algebraic[34])/1000.00)*0.00100000 algebraic[24] = (1.00000/constants[31])*(states[14]-constants[34]) algebraic[25] = (algebraic[21]-algebraic[24])/constants[32] algebraic[27] = (1.00000/constants[35])*(states[16]-constants[38]) algebraic[28] = (algebraic[21]-algebraic[27])/constants[36] algebraic[35] = (algebraic[21]-algebraic[33])/constants[44] algebraic[37] = (1.00000/constants[48])*(states[22]-constants[50]) algebraic[38] = (algebraic[21]-algebraic[37])/constants[49] algebraic[43] = (1.00000/constants[56])*(states[26]-constants[59]) algebraic[44] = (algebraic[21]-algebraic[43])/constants[57] algebraic[40] = (1.00000/constants[52])*(states[24]-constants[55]) algebraic[41] = (algebraic[21]-algebraic[40])/constants[53] rates[13] = (algebraic[10]*algebraic[22]-algebraic[20]*(algebraic[41]+algebraic[25]+algebraic[28]+algebraic[31]+algebraic[35]+algebraic[38]+algebraic[44]))/1000.00 rates[12] = (((((((algebraic[22]-algebraic[25])-algebraic[28])-algebraic[31])-algebraic[35])-algebraic[38])-algebraic[44])-algebraic[41])/1000.00 algebraic[46] = (1.00000/constants[61])*(states[28]-constants[62]) algebraic[47] = (algebraic[24]-algebraic[46])/constants[33] algebraic[23] = states[15]/states[14] rates[15] = ((algebraic[20]*algebraic[25]-algebraic[23]*algebraic[47])/1000.00)*0.00100000 rates[14] = ((algebraic[25]-algebraic[47])/1000.00)*0.00100000 algebraic[48] = (algebraic[27]-algebraic[46])/constants[37] algebraic[26] = states[17]/states[16] rates[17] = ((algebraic[20]*algebraic[28]-algebraic[26]*algebraic[48])/1000.00)*0.00100000 rates[16] = ((algebraic[28]-algebraic[48])/1000.00)*0.00100000 algebraic[49] = (algebraic[33]-algebraic[46])/constants[45] algebraic[32] = states[21]/states[20] rates[21] = ((((algebraic[29]*algebraic[34]+algebraic[20]*algebraic[35])-algebraic[32]*algebraic[49])+algebraic[32]*constants[47])/1000.00)*0.00100000 rates[20] = (((algebraic[34]+algebraic[35])-algebraic[49])/1000.00)*0.00100000 algebraic[50] = (algebraic[37]-algebraic[46])/constants[51] algebraic[36] = states[23]/states[22] rates[23] = ((algebraic[20]*algebraic[38]-algebraic[36]*algebraic[50])/1000.00)*0.00100000 rates[22] = ((algebraic[38]-algebraic[50])/1000.00)*0.00100000 algebraic[51] = (algebraic[40]-algebraic[46])/constants[54] algebraic[39] = states[25]/states[24] rates[25] = ((algebraic[20]*algebraic[41]-algebraic[39]*algebraic[51])/1000.00)*0.00100000 rates[24] = ((algebraic[41]-algebraic[51])/1000.00)*0.00100000 algebraic[52] = (algebraic[43]-algebraic[46])/constants[58] algebraic[42] = states[27]/states[26] rates[27] = ((algebraic[20]*algebraic[44]-algebraic[42]*algebraic[52])/1000.00)*0.00100000 rates[26] = ((algebraic[44]-algebraic[52])/1000.00)*0.00100000 algebraic[55] = (1.00000/constants[65])*(states[30]-constants[64]) algebraic[56] = (algebraic[18]-algebraic[55])/constants[23] rates[10] = ((algebraic[10]*algebraic[19]-algebraic[17]*algebraic[56])-(algebraic[17]-algebraic[16])*constants[26])/1000.00 rates[9] = (algebraic[19]-algebraic[56])/1000.00 algebraic[57] = (algebraic[46]-algebraic[55])/constants[60] algebraic[45] = states[29]/states[28] rates[29] = ((algebraic[39]*algebraic[51]+algebraic[23]*algebraic[47]+algebraic[26]*algebraic[48]+algebraic[32]*algebraic[49]+algebraic[36]*algebraic[50]+algebraic[42]*algebraic[52])-algebraic[45]*algebraic[57])/1000.00 algebraic[53] = algebraic[51]+algebraic[47]+algebraic[48]+algebraic[49]+algebraic[50]+algebraic[52] rates[28] = (algebraic[53]-algebraic[57])/1000.00 algebraic[59] = (voi-floor(voi/constants[91])*constants[91])/1000.00 algebraic[60] = custom_piecewise([greater_equal(algebraic[59]*0.00100000 , 0.00000) & less_equal(algebraic[59]*0.00100000 , constants[94]), constants[67]+((constants[66]-constants[67])*(1.00000-cos(( pi*(algebraic[59]*0.00100000))/constants[94])))/2.00000 , less(algebraic[59]*0.00100000 , 1.50000*constants[94]) & greater_equal(algebraic[59]*0.00100000 , constants[94]), constants[67]+((constants[66]-constants[67])*(1.00000+cos((2.00000* pi*(algebraic[59]*0.00100000-constants[94]))/constants[94])))/2.00000 , True, constants[67]]) algebraic[62] = algebraic[60]*(states[32]-constants[70])+constants[69] algebraic[63] = (algebraic[14]-algebraic[62])/constants[19] algebraic[13] = states[8]/states[7] rates[8] = (algebraic[6]*algebraic[15]-algebraic[13]*algebraic[63])/1000.00 rates[7] = (algebraic[15]-algebraic[63])/1000.00 algebraic[64] = (algebraic[55]-algebraic[62])/constants[63] algebraic[54] = states[31]/states[30] rates[31] = ((algebraic[45]*algebraic[57]+algebraic[17]*algebraic[56])-algebraic[54]*algebraic[64])/1000.00 algebraic[58] = algebraic[56]+algebraic[57] rates[30] = (algebraic[58]-algebraic[64])/1000.00 algebraic[66] = (voi-floor(voi/constants[91])*constants[91])/1000.00-constants[71] algebraic[67] = custom_piecewise([greater_equal(algebraic[66]*0.00100000 , 0.00000) & less_equal(algebraic[66]*0.00100000 , constants[95]), constants[73]+((constants[72]-constants[73])*(1.00000-cos(( pi*(algebraic[66]*0.00100000))/constants[95])))/2.00000 , less(algebraic[66]*0.00100000 , 1.50000*constants[95]) & greater_equal(algebraic[66]*0.00100000 , constants[95]), constants[73]+((constants[72]-constants[73])*(1.00000+cos((2.00000* pi*(algebraic[66]*0.00100000-constants[95]))/constants[95])))/2.00000 , True, constants[73]]) algebraic[69] = algebraic[67]*(states[34]-constants[76])+constants[75] algebraic[70] = custom_piecewise([greater(algebraic[62] , algebraic[69]), (algebraic[62]-algebraic[69])/constants[77] , True, 0.00000]) algebraic[61] = states[33]/states[32] rates[33] = ((algebraic[54]*algebraic[64]+algebraic[13]*algebraic[63])-algebraic[61]*algebraic[70])/1000.00 algebraic[65] = algebraic[64]+algebraic[63] rates[32] = (algebraic[65]-algebraic[70])/1000.00 algebraic[72] = (1.00000/constants[79])*(states[36]-constants[80]) algebraic[73] = custom_piecewise([greater(algebraic[69] , algebraic[72]), (algebraic[69]-algebraic[72])/constants[78] , True, 0.00000]) algebraic[68] = states[35]/states[34] rates[35] = (algebraic[61]*algebraic[70]-algebraic[68]*algebraic[73])/1000.00 rates[34] = (algebraic[70]-algebraic[73])/1000.00 algebraic[75] = (1.00000/constants[81])*(states[38]-constants[83]) algebraic[76] = (algebraic[72]-algebraic[75])/constants[82] algebraic[71] = states[37]/states[36] rates[37] = (algebraic[68]*algebraic[73]-algebraic[71]*algebraic[76])/1000.00 rates[36] = (algebraic[73]-algebraic[76])/1000.00 algebraic[78] = (1.00000/constants[85])*(states[40]-constants[88]) algebraic[80] = (algebraic[75]-algebraic[78])/constants[86] algebraic[79] = (algebraic[75]-algebraic[78])/constants[84] algebraic[74] = states[39]/states[38] rates[39] = (algebraic[71]*algebraic[76]-algebraic[74]*(algebraic[80]+algebraic[79]))/1000.00 rates[38] = ((algebraic[76]-algebraic[80])-algebraic[79])/1000.00 algebraic[82] = (1.00000/constants[89])*(states[42]-constants[90]) algebraic[83] = (algebraic[82]-algebraic[3])/constants[6] algebraic[81] = states[43]/states[42] rates[1] = (algebraic[81]*algebraic[83]-algebraic[2]*algebraic[9])/1000.00 rates[0] = (algebraic[83]-algebraic[9])/1000.00 algebraic[84] = (algebraic[78]-algebraic[82])/constants[87] algebraic[77] = states[41]/states[40] rates[41] = (algebraic[74]*algebraic[80]-algebraic[77]*algebraic[84])/1000.00 rates[40] = (algebraic[80]-algebraic[84])/1000.00 rates[43] = ((algebraic[77]*algebraic[84]+algebraic[74]*algebraic[79])-algebraic[81]*algebraic[83])/1000.00 algebraic[85] = algebraic[84]+algebraic[79] rates[42] = (algebraic[85]-algebraic[83])/1000.00 return(rates) def computeAlgebraic(constants, states, voi): algebraic = array([[0.0] * len(voi)] * sizeAlgebraic) states = array(states) voi = array(voi) algebraic[4] = (voi-floor(voi/constants[91])*constants[91])/1000.00-constants[7] algebraic[5] = custom_piecewise([greater_equal(algebraic[4]*0.00100000 , 0.00000) & less_equal(algebraic[4]*0.00100000 , constants[93]), constants[9]+((constants[8]-constants[9])*(1.00000-cos(( pi*(algebraic[4]*0.00100000))/constants[93])))/2.00000 , less(algebraic[4]*0.00100000 , 1.50000*constants[93]) & greater_equal(algebraic[4]*0.00100000 , constants[93]), constants[9]+((constants[8]-constants[9])*(1.00000+cos((2.00000* pi*(algebraic[4]*0.00100000-constants[93]))/constants[93])))/2.00000 , True, constants[9]]) algebraic[7] = algebraic[5]*(states[2]-constants[12])+constants[11] algebraic[11] = (1.00000/constants[14])*(states[4]-constants[17]) algebraic[0] = (voi-floor(voi/constants[91])*constants[91])/1000.00 algebraic[1] = custom_piecewise([greater_equal(algebraic[0] , 0.00000) & less_equal(algebraic[0] , constants[92]), constants[2]+((constants[1]-constants[2])*(1.00000-cos(( pi*algebraic[0])/constants[92])))/2.00000 , less(algebraic[0] , 1.50000*constants[92]) & greater_equal(algebraic[0] , constants[92]), constants[2]+((constants[1]-constants[2])*(1.00000+cos((2.00000* pi*(algebraic[0]-constants[92]))/constants[92])))/2.00000 , True, constants[2]]) algebraic[3] = algebraic[1]*(states[0]-constants[5])+constants[4] algebraic[9] = custom_piecewise([greater(algebraic[3] , algebraic[7]), (algebraic[3]-algebraic[7])/constants[13] , True, 0.00000]) algebraic[2] = states[1]/states[0] algebraic[12] = custom_piecewise([greater(algebraic[7] , algebraic[11]), states[5] , True, 0.00000]) algebraic[14] = (1.00000/constants[21])*(states[7]-constants[20]) algebraic[15] = custom_piecewise([greater(algebraic[7] , algebraic[14]), (algebraic[7]-algebraic[14])/constants[18] , True, 0.00000]) algebraic[6] = states[3]/states[2] algebraic[17] = states[10]/states[9] algebraic[16] = states[11]/constants[27] algebraic[18] = (1.00000/constants[25])*(states[9]-constants[24]) algebraic[19] = (algebraic[11]-algebraic[18])/constants[22] algebraic[21] = (1.00000/constants[29])*(states[12]-constants[30]) algebraic[22] = (algebraic[11]-algebraic[21])/constants[28] algebraic[10] = states[6]/states[4] algebraic[30] = (1.00000/constants[39])*(states[18]-constants[42]) algebraic[31] = (algebraic[21]-algebraic[30])/constants[40] algebraic[20] = states[13]/states[12] algebraic[33] = (1.00000/constants[43])*(states[20]-constants[46]) algebraic[34] = (algebraic[30]-algebraic[33])/constants[41] algebraic[29] = states[19]/states[18] algebraic[24] = (1.00000/constants[31])*(states[14]-constants[34]) algebraic[25] = (algebraic[21]-algebraic[24])/constants[32] algebraic[27] = (1.00000/constants[35])*(states[16]-constants[38]) algebraic[28] = (algebraic[21]-algebraic[27])/constants[36] algebraic[35] = (algebraic[21]-algebraic[33])/constants[44] algebraic[37] = (1.00000/constants[48])*(states[22]-constants[50]) algebraic[38] = (algebraic[21]-algebraic[37])/constants[49] algebraic[43] = (1.00000/constants[56])*(states[26]-constants[59]) algebraic[44] = (algebraic[21]-algebraic[43])/constants[57] algebraic[40] = (1.00000/constants[52])*(states[24]-constants[55]) algebraic[41] = (algebraic[21]-algebraic[40])/constants[53] algebraic[46] = (1.00000/constants[61])*(states[28]-constants[62]) algebraic[47] = (algebraic[24]-algebraic[46])/constants[33] algebraic[23] = states[15]/states[14] algebraic[48] = (algebraic[27]-algebraic[46])/constants[37] algebraic[26] = states[17]/states[16] algebraic[49] = (algebraic[33]-algebraic[46])/constants[45] algebraic[32] = states[21]/states[20] algebraic[50] = (algebraic[37]-algebraic[46])/constants[51] algebraic[36] = states[23]/states[22] algebraic[51] = (algebraic[40]-algebraic[46])/constants[54] algebraic[39] = states[25]/states[24] algebraic[52] = (algebraic[43]-algebraic[46])/constants[58] algebraic[42] = states[27]/states[26] algebraic[55] = (1.00000/constants[65])*(states[30]-constants[64]) algebraic[56] = (algebraic[18]-algebraic[55])/constants[23] algebraic[57] = (algebraic[46]-algebraic[55])/constants[60] algebraic[45] = states[29]/states[28] algebraic[53] = algebraic[51]+algebraic[47]+algebraic[48]+algebraic[49]+algebraic[50]+algebraic[52] algebraic[59] = (voi-floor(voi/constants[91])*constants[91])/1000.00 algebraic[60] = custom_piecewise([greater_equal(algebraic[59]*0.00100000 , 0.00000) & less_equal(algebraic[59]*0.00100000 , constants[94]), constants[67]+((constants[66]-constants[67])*(1.00000-cos(( pi*(algebraic[59]*0.00100000))/constants[94])))/2.00000 , less(algebraic[59]*0.00100000 , 1.50000*constants[94]) & greater_equal(algebraic[59]*0.00100000 , constants[94]), constants[67]+((constants[66]-constants[67])*(1.00000+cos((2.00000* pi*(algebraic[59]*0.00100000-constants[94]))/constants[94])))/2.00000 , True, constants[67]]) algebraic[62] = algebraic[60]*(states[32]-constants[70])+constants[69] algebraic[63] = (algebraic[14]-algebraic[62])/constants[19] algebraic[13] = states[8]/states[7] algebraic[64] = (algebraic[55]-algebraic[62])/constants[63] algebraic[54] = states[31]/states[30] algebraic[58] = algebraic[56]+algebraic[57] algebraic[66] = (voi-floor(voi/constants[91])*constants[91])/1000.00-constants[71] algebraic[67] = custom_piecewise([greater_equal(algebraic[66]*0.00100000 , 0.00000) & less_equal(algebraic[66]*0.00100000 , constants[95]), constants[73]+((constants[72]-constants[73])*(1.00000-cos(( pi*(algebraic[66]*0.00100000))/constants[95])))/2.00000 , less(algebraic[66]*0.00100000 , 1.50000*constants[95]) & greater_equal(algebraic[66]*0.00100000 , constants[95]), constants[73]+((constants[72]-constants[73])*(1.00000+cos((2.00000* pi*(algebraic[66]*0.00100000-constants[95]))/constants[95])))/2.00000 , True, constants[73]]) algebraic[69] = algebraic[67]*(states[34]-constants[76])+constants[75] algebraic[70] = custom_piecewise([greater(algebraic[62] , algebraic[69]), (algebraic[62]-algebraic[69])/constants[77] , True, 0.00000]) algebraic[61] = states[33]/states[32] algebraic[65] = algebraic[64]+algebraic[63] algebraic[72] = (1.00000/constants[79])*(states[36]-constants[80]) algebraic[73] = custom_piecewise([greater(algebraic[69] , algebraic[72]), (algebraic[69]-algebraic[72])/constants[78] , True, 0.00000]) algebraic[68] = states[35]/states[34] algebraic[75] = (1.00000/constants[81])*(states[38]-constants[83]) algebraic[76] = (algebraic[72]-algebraic[75])/constants[82] algebraic[71] = states[37]/states[36] algebraic[78] = (1.00000/constants[85])*(states[40]-constants[88]) algebraic[80] = (algebraic[75]-algebraic[78])/constants[86] algebraic[79] = (algebraic[75]-algebraic[78])/constants[84] algebraic[74] = states[39]/states[38] algebraic[82] = (1.00000/constants[89])*(states[42]-constants[90]) algebraic[83] = (algebraic[82]-algebraic[3])/constants[6] algebraic[81] = states[43]/states[42] algebraic[84] = (algebraic[78]-algebraic[82])/constants[87] algebraic[77] = states[41]/states[40] algebraic[85] = algebraic[84]+algebraic[79] algebraic[8] = algebraic[7]*(1.00000+(3.00000*states[2])/200.000) 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)