# Size of variable arrays: sizeAlgebraic = 93 sizeStates = 33 sizeConstants = 139 from math import * from numpy import * def createLegends(): legend_states = [""] * sizeStates legend_rates = [""] * sizeStates legend_algebraic = [""] * sizeAlgebraic legend_voi = "" legend_constants = [""] * sizeConstants legend_voi = "t in component main (second)" legend_constants[0] = "rho in component main (Js2_per_m5)" legend_constants[1] = "g in component main (m_per_s2)" legend_constants[2] = "theta_deg in component main (dim)" legend_constants[138] = "theta_rad in component main (radian)" legend_algebraic[0] = "z in component main (dim)" legend_constants[3] = "T in component main (second)" legend_algebraic[2] = "mt in component main (dim)" legend_algebraic[1] = "mta in component main (dim)" legend_constants[4] = "delay in component main (dim)" legend_constants[5] = "To_pk in component main (dim)" legend_algebraic[4] = "E_B_LV in component main (kPa_per_L)" legend_constants[6] = "E_B_LVmax in component main (kPa_per_L)" legend_constants[7] = "E_B_LVmin in component main (kPa_per_L)" legend_states[0] = "q_B_LV in component main (litre)" legend_constants[8] = "q_B_LV_us in component main (litre)" legend_algebraic[6] = "u_B_LV in component main (kPa)" legend_algebraic[8] = "E_B_RV in component main (kPa_per_L)" legend_constants[9] = "E_B_RVmax in component main (kPa_per_L)" legend_constants[10] = "E_B_RVmin in component main (kPa_per_L)" legend_states[1] = "q_B_RV in component main (litre)" legend_constants[11] = "q_B_RV_us in component main (litre)" legend_algebraic[10] = "u_B_RV in component main (kPa)" legend_constants[12] = "E_B_LA in component main (kPa_per_L)" legend_states[2] = "q_B_LA in component main (litre)" legend_constants[13] = "q_B_LA_us in component main (litre)" legend_algebraic[12] = "u_B_LA in component main (kPa)" legend_constants[14] = "E_B_RA in component main (kPa_per_L)" legend_states[3] = "q_B_RA in component main (litre)" legend_constants[15] = "q_B_RA_us in component main (litre)" legend_algebraic[14] = "u_B_RA in component main (kPa)" legend_constants[16] = "R_B_AV in component main (kPa_s_per_L)" legend_algebraic[29] = "v_B_AV in component main (L_per_s)" legend_constants[17] = "R_B_PV in component main (kPa_s_per_L)" legend_algebraic[20] = "v_B_PV in component main (L_per_s)" legend_constants[18] = "R_B_TV in component main (kPa_s_per_L)" legend_algebraic[16] = "v_B_TV in component main (L_per_s)" legend_constants[19] = "R_B_MV in component main (kPa_s_per_L)" legend_algebraic[17] = "v_B_MV in component main (L_per_s)" legend_states[4] = "q_B_PA in component main (litre)" legend_constants[20] = "q_B_PA_us in component main (litre)" legend_algebraic[18] = "u_B_PA in component main (kPa)" legend_constants[21] = "E_B_PA in component main (kPa_per_L)" legend_states[5] = "q_B_lung in component main (litre)" legend_constants[22] = "q_B_lung_us in component main (litre)" legend_algebraic[22] = "u_B_lung in component main (kPa)" legend_constants[23] = "E_B_lung in component main (kPa_per_L)" legend_constants[24] = "z_lung in component main (meter)" legend_algebraic[25] = "v_B_lung1 in component main (L_per_s)" legend_algebraic[28] = "v_B_lung2 in component main (L_per_s)" legend_constants[25] = "R_B_lung1 in component main (kPa_s_per_L)" legend_constants[26] = "R_B_lung2 in component main (kPa_s_per_L)" legend_constants[27] = "E_B_pulmVein in component main (kPa_per_L)" legend_states[6] = "q_B_pulmVein in component main (litre)" legend_constants[28] = "q_B_pulmVein_us in component main (litre)" legend_algebraic[27] = "u_B_pulmVein in component main (kPa)" legend_algebraic[30] = "v_B_pulmVein in component main (L_per_s)" legend_constants[29] = "R_B_pulmVein in component main (kPa_s_per_L)" legend_states[7] = "q_B_brain in component main (litre)" legend_constants[30] = "q_B_brain_us in component main (litre)" legend_states[8] = "q_B_brainVein in component main (litre)" legend_constants[31] = "q_B_brainVein_us in component main (litre)" legend_algebraic[19] = "u_B_brain in component main (kPa)" legend_algebraic[21] = "u_B_brainVein in component main (kPa)" legend_constants[32] = "E_B_brain in component main (kPa_per_L)" legend_constants[33] = "E_B_brainVein in component main (kPa_per_L)" legend_algebraic[31] = "v_B_brain1 in component main (L_per_s)" legend_algebraic[23] = "v_B_brain2 in component main (L_per_s)" legend_algebraic[24] = "v_B_brain3 in component main (L_per_s)" legend_constants[34] = "R_B_brain1 in component main (kPa_s_per_L)" legend_constants[35] = "R_B_brain2 in component main (kPa_s_per_L)" legend_constants[36] = "R_B_brain3 in component main (kPa_s_per_L)" legend_states[9] = "q_B_AA in component main (litre)" legend_constants[37] = "q_B_AA_us in component main (litre)" legend_algebraic[26] = "u_B_AA in component main (kPa)" legend_constants[38] = "E_B_AA in component main (kPa_per_L)" legend_states[10] = "q_B_celiac in component main (litre)" legend_constants[39] = "q_B_celiac_us in component main (litre)" legend_algebraic[32] = "u_B_celiac in component main (kPa)" legend_constants[40] = "E_B_celiac in component main (kPa_per_L)" legend_constants[41] = "z_celiac in component main (meter)" legend_algebraic[33] = "v_B_celiac in component main (L_per_s)" legend_constants[42] = "R_B_celiac in component main (kPa_s_per_L)" legend_states[11] = "q_B_supMes in component main (litre)" legend_constants[43] = "q_B_supMes_us in component main (litre)" legend_algebraic[34] = "u_B_supMes in component main (kPa)" legend_constants[44] = "E_B_supMes in component main (kPa_per_L)" legend_constants[45] = "z_supMes in component main (meter)" legend_algebraic[35] = "v_B_supMes in component main (L_per_s)" legend_constants[46] = "R_B_supMes in component main (kPa_s_per_L)" legend_algebraic[46] = "v_B_infMes in component main (L_per_s)" legend_constants[47] = "R_B_infMes in component main (kPa_s_per_L)" legend_states[12] = "q_B_stomach in component main (litre)" legend_constants[48] = "q_B_stomach_us in component main (litre)" legend_algebraic[36] = "u_B_stomach in component main (kPa)" legend_constants[49] = "E_B_stomach in component main (kPa_per_L)" legend_constants[50] = "z_stomach in component main (meter)" legend_algebraic[37] = "v_B_stomach1 in component main (L_per_s)" legend_algebraic[49] = "v_B_stomach2 in component main (L_per_s)" legend_constants[51] = "R_B_stomach1 in component main (kPa_s_per_L)" legend_constants[52] = "R_B_stomach2 in component main (kPa_s_per_L)" legend_states[13] = "q_B_spleen in component main (litre)" legend_constants[53] = "q_B_spleen_us in component main (litre)" legend_algebraic[38] = "u_B_spleen in component main (kPa)" legend_constants[54] = "E_B_spleen in component main (kPa_per_L)" legend_constants[55] = "z_spleen in component main (meter)" legend_algebraic[39] = "v_B_spleen1 in component main (L_per_s)" legend_algebraic[50] = "v_B_spleen2 in component main (L_per_s)" legend_constants[56] = "R_B_spleen1 in component main (kPa_s_per_L)" legend_constants[57] = "R_B_spleen2 in component main (kPa_s_per_L)" legend_states[14] = "q_B_pancreas in component main (litre)" legend_constants[58] = "q_B_pancreas_us in component main (litre)" legend_algebraic[40] = "u_B_pancreas in component main (kPa)" legend_constants[59] = "E_B_pancreas in component main (kPa_per_L)" legend_constants[60] = "z_pancreas in component main (meter)" legend_algebraic[41] = "v_B_pancreas1 in component main (L_per_s)" legend_algebraic[42] = "v_B_pancreas2 in component main (L_per_s)" legend_algebraic[51] = "v_B_pancreas3 in component main (L_per_s)" legend_constants[61] = "R_B_pancreas1 in component main (kPa_s_per_L)" legend_constants[62] = "R_B_pancreas2 in component main (kPa_s_per_L)" legend_constants[63] = "R_B_pancreas3 in component main (kPa_s_per_L)" legend_states[15] = "q_B_intestine in component main (litre)" legend_constants[64] = "q_B_intestine_us in component main (litre)" legend_algebraic[43] = "u_B_intestine in component main (kPa)" legend_constants[65] = "E_B_intestine in component main (kPa_per_L)" legend_constants[66] = "z_ in component main (meter)" legend_algebraic[44] = "v_B_intestine1 in component main (L_per_s)" legend_algebraic[52] = "v_B_intestine2 in component main (L_per_s)" legend_constants[67] = "R_B_intestine1 in component main (kPa_s_per_L)" legend_constants[68] = "R_B_intestine2 in component main (kPa_s_per_L)" legend_states[16] = "q_B_colon in component main (litre)" legend_constants[69] = "q_B_colon_us in component main (litre)" legend_algebraic[45] = "u_B_colon in component main (kPa)" legend_constants[70] = "E_B_colon in component main (kPa_per_L)" legend_constants[71] = "z_colon in component main (meter)" legend_algebraic[47] = "v_B_colon1 in component main (L_per_s)" legend_algebraic[53] = "v_B_colon2 in component main (L_per_s)" legend_constants[72] = "R_B_colon1 in component main (kPa_s_per_L)" legend_constants[73] = "R_B_colon2 in component main (kPa_s_per_L)" legend_states[17] = "q_B_portalVein in component main (litre)" legend_constants[74] = "q_B_portalVein_us in component main (litre)" legend_algebraic[48] = "u_B_portalVein in component main (kPa)" legend_constants[75] = "E_B_portalVein in component main (kPa_per_L)" legend_constants[76] = "z_portalVein in component main (meter)" legend_algebraic[55] = "v_B_portalVein in component main (L_per_s)" legend_constants[77] = "R_B_portalVein in component main (kPa_s_per_L)" legend_states[18] = "q_B_liver in component main (litre)" legend_constants[78] = "q_B_liver_us in component main (litre)" legend_algebraic[54] = "u_B_liver in component main (kPa)" legend_constants[79] = "E_B_liver in component main (kPa_per_L)" legend_constants[80] = "z_liver in component main (meter)" legend_algebraic[56] = "v_B_liver1 in component main (L_per_s)" legend_algebraic[57] = "v_B_liver2 in component main (L_per_s)" legend_constants[81] = "R_B_liver1 in component main (kPa_s_per_L)" legend_constants[82] = "R_B_liver2 in component main (kPa_s_per_L)" legend_states[19] = "q_B_lGlomerulus in component main (litre)" legend_states[20] = "q_B_rGlomerulus in component main (litre)" legend_states[21] = "q_B_lPeritubular in component main (litre)" legend_states[22] = "q_B_rPeritubular in component main (litre)" legend_constants[83] = "q_B_lGlomerulus_us in component main (litre)" legend_constants[84] = "q_B_rGlomerulus_us in component main (litre)" legend_constants[85] = "q_B_lPeritubular_us in component main (litre)" legend_constants[86] = "q_B_rPeritubular_us in component main (litre)" legend_algebraic[58] = "u_B_lGlomerulus in component main (kPa)" legend_algebraic[59] = "u_B_rGlomerulus in component main (kPa)" legend_algebraic[60] = "u_B_lPeritubular in component main (kPa)" legend_algebraic[61] = "u_B_rPeritubular in component main (kPa)" legend_constants[87] = "E_B_lGlomerulus in component main (kPa_per_L)" legend_constants[88] = "E_B_rGlomerulus in component main (kPa_per_L)" legend_constants[89] = "E_B_lPeritubular in component main (kPa_per_L)" legend_constants[90] = "E_B_rPeritubular in component main (kPa_per_L)" legend_constants[91] = "z_kidney in component main (meter)" legend_algebraic[63] = "v_B_lRenal1 in component main (L_per_s)" legend_algebraic[64] = "v_B_rRenal1 in component main (L_per_s)" legend_algebraic[65] = "v_B_lRenal2 in component main (L_per_s)" legend_algebraic[62] = "v_B_rRenal2 in component main (L_per_s)" legend_algebraic[66] = "v_B_lRenal3 in component main (L_per_s)" legend_algebraic[67] = "v_B_rRenal3 in component main (L_per_s)" legend_constants[92] = "R_B_lRenal1 in component main (kPa_s_per_L)" legend_constants[93] = "R_B_rRenal1 in component main (kPa_s_per_L)" legend_constants[94] = "R_B_lRenal2 in component main (kPa_s_per_L)" legend_constants[95] = "R_B_rRenal2 in component main (kPa_s_per_L)" legend_constants[96] = "R_B_lRenal3 in component main (kPa_s_per_L)" legend_constants[97] = "R_B_rRenal3 in component main (kPa_s_per_L)" legend_states[23] = "q_B_musTorso in component main (litre)" legend_states[24] = "q_B_musTorsoVein in component main (litre)" legend_constants[98] = "q_B_musTorso_us in component main (litre)" legend_constants[99] = "q_B_musTorsoVein_us in component main (litre)" legend_algebraic[68] = "u_B_musTorso in component main (kPa)" legend_algebraic[69] = "u_B_musTorsoVein in component main (kPa)" legend_constants[100] = "E_B_musTorso in component main (kPa_per_L)" legend_constants[101] = "E_B_musTorsoVein in component main (kPa_per_L)" legend_constants[102] = "z_musTorso in component main (meter)" legend_algebraic[70] = "v_B_musTorso1 in component main (L_per_s)" legend_algebraic[71] = "v_B_musTorso2 in component main (L_per_s)" legend_algebraic[72] = "v_B_musTorso3 in component main (L_per_s)" legend_constants[103] = "R_B_musTorso1 in component main (kPa_s_per_L)" legend_constants[104] = "R_B_musTorso2 in component main (kPa_s_per_L)" legend_constants[105] = "R_B_musTorso3 in component main (kPa_s_per_L)" legend_states[25] = "q_B_lArm in component main (litre)" legend_states[26] = "q_B_lArmVein in component main (litre)" legend_constants[106] = "q_B_lArm_us in component main (litre)" legend_constants[107] = "q_B_lArmVein_us in component main (litre)" legend_algebraic[73] = "u_B_lArm in component main (kPa)" legend_algebraic[74] = "u_B_lArmVein in component main (kPa)" legend_constants[108] = "E_B_lArm in component main (kPa_per_L)" legend_constants[109] = "E_B_lArmVein in component main (kPa_per_L)" legend_constants[110] = "z_lArm in component main (meter)" legend_algebraic[75] = "v_B_lArm1 in component main (L_per_s)" legend_algebraic[76] = "v_B_lArm2 in component main (L_per_s)" legend_algebraic[77] = "v_B_lArm3 in component main (L_per_s)" legend_constants[111] = "R_B_lArm1 in component main (kPa_s_per_L)" legend_constants[112] = "R_B_lArm2 in component main (kPa_s_per_L)" legend_constants[113] = "R_B_lArm3 in component main (kPa_s_per_L)" legend_states[27] = "q_B_rArm in component main (litre)" legend_states[28] = "q_B_rArmVein in component main (litre)" legend_constants[114] = "q_B_rArm_us in component main (litre)" legend_constants[115] = "q_B_rArmVein_us in component main (litre)" legend_algebraic[78] = "u_B_rArm in component main (kPa)" legend_algebraic[79] = "u_B_rArmVein in component main (kPa)" legend_constants[116] = "E_B_rArm in component main (kPa_per_L)" legend_constants[117] = "E_B_rArmVein in component main (kPa_per_L)" legend_constants[118] = "z_rArm in component main (meter)" legend_algebraic[80] = "v_B_rArm1 in component main (L_per_s)" legend_algebraic[81] = "v_B_rArm2 in component main (L_per_s)" legend_algebraic[82] = "v_B_rArm3 in component main (L_per_s)" legend_constants[119] = "R_B_rArm1 in component main (kPa_s_per_L)" legend_constants[120] = "R_B_rArm2 in component main (kPa_s_per_L)" legend_constants[121] = "R_B_rArm3 in component main (kPa_s_per_L)" legend_states[29] = "q_B_lLeg in component main (litre)" legend_states[30] = "q_B_lLegVein in component main (litre)" legend_constants[122] = "q_B_lLeg_us in component main (litre)" legend_constants[123] = "q_B_lLegVein_us in component main (litre)" legend_algebraic[83] = "u_B_lLeg in component main (kPa)" legend_algebraic[84] = "u_B_lLegVein in component main (kPa)" legend_constants[124] = "E_B_lLeg in component main (kPa_per_L)" legend_constants[125] = "E_B_lLegVein in component main (kPa_per_L)" legend_constants[126] = "z_lLeg in component main (meter)" legend_algebraic[85] = "v_B_lLeg1 in component main (L_per_s)" legend_algebraic[86] = "v_B_lLeg2 in component main (L_per_s)" legend_algebraic[87] = "v_B_lLeg3 in component main (L_per_s)" legend_constants[127] = "R_B_lLeg1 in component main (kPa_s_per_L)" legend_constants[128] = "R_B_lLeg2 in component main (kPa_s_per_L)" legend_constants[129] = "R_B_lLeg3 in component main (kPa_s_per_L)" legend_states[31] = "q_B_rLeg in component main (litre)" legend_states[32] = "q_B_rLegVein in component main (litre)" legend_constants[130] = "q_B_rLeg_us in component main (litre)" legend_constants[131] = "q_B_rLegVein_us in component main (litre)" legend_algebraic[88] = "u_B_rLeg in component main (kPa)" legend_algebraic[89] = "u_B_rLegVein in component main (kPa)" legend_constants[132] = "E_B_rLeg in component main (kPa_per_L)" legend_constants[133] = "E_B_rLegVein in component main (kPa_per_L)" legend_constants[134] = "z_rLeg in component main (meter)" legend_algebraic[90] = "v_B_rLeg1 in component main (L_per_s)" legend_algebraic[91] = "v_B_rLeg2 in component main (L_per_s)" legend_algebraic[92] = "v_B_rLeg3 in component main (L_per_s)" legend_constants[135] = "R_B_rLeg1 in component main (kPa_s_per_L)" legend_constants[136] = "R_B_rLeg2 in component main (kPa_s_per_L)" legend_constants[137] = "R_B_rLeg3 in component main (kPa_s_per_L)" legend_algebraic[3] = "q_B_heartTot in component main (litre)" legend_algebraic[5] = "q_B_lungTot in component main (litre)" legend_algebraic[7] = "q_B_renalTot in component main (litre)" legend_algebraic[9] = "q_B_gutTot in component main (litre)" legend_algebraic[11] = "q_B_brainTot in component main (litre)" legend_algebraic[13] = "q_B_limbTot in component main (litre)" legend_algebraic[15] = "q_Blood in component main (litre)" legend_rates[0] = "d/dt q_B_LV in component main (litre)" legend_rates[1] = "d/dt q_B_RV in component main (litre)" legend_rates[2] = "d/dt q_B_LA in component main (litre)" legend_rates[3] = "d/dt q_B_RA in component main (litre)" legend_rates[4] = "d/dt q_B_PA in component main (litre)" legend_rates[5] = "d/dt q_B_lung in component main (litre)" legend_rates[6] = "d/dt q_B_pulmVein in component main (litre)" legend_rates[7] = "d/dt q_B_brain in component main (litre)" legend_rates[8] = "d/dt q_B_brainVein in component main (litre)" legend_rates[9] = "d/dt q_B_AA in component main (litre)" legend_rates[10] = "d/dt q_B_celiac in component main (litre)" legend_rates[11] = "d/dt q_B_supMes in component main (litre)" legend_rates[12] = "d/dt q_B_stomach in component main (litre)" legend_rates[13] = "d/dt q_B_spleen in component main (litre)" legend_rates[14] = "d/dt q_B_pancreas in component main (litre)" legend_rates[15] = "d/dt q_B_intestine in component main (litre)" legend_rates[16] = "d/dt q_B_colon in component main (litre)" legend_rates[17] = "d/dt q_B_portalVein in component main (litre)" legend_rates[18] = "d/dt q_B_liver in component main (litre)" legend_rates[19] = "d/dt q_B_lGlomerulus in component main (litre)" legend_rates[20] = "d/dt q_B_rGlomerulus in component main (litre)" legend_rates[21] = "d/dt q_B_lPeritubular in component main (litre)" legend_rates[22] = "d/dt q_B_rPeritubular in component main (litre)" legend_rates[23] = "d/dt q_B_musTorso in component main (litre)" legend_rates[24] = "d/dt q_B_musTorsoVein in component main (litre)" legend_rates[25] = "d/dt q_B_lArm in component main (litre)" legend_rates[26] = "d/dt q_B_lArmVein in component main (litre)" legend_rates[27] = "d/dt q_B_rArm in component main (litre)" legend_rates[28] = "d/dt q_B_rArmVein in component main (litre)" legend_rates[29] = "d/dt q_B_lLeg in component main (litre)" legend_rates[30] = "d/dt q_B_lLegVein in component main (litre)" legend_rates[31] = "d/dt q_B_rLeg in component main (litre)" legend_rates[32] = "d/dt q_B_rLegVein in component main (litre)" return (legend_states, legend_algebraic, legend_voi, legend_constants) def initConsts(): constants = [0.0] * sizeConstants; states = [0.0] * sizeStates; constants[0] = 1e3 constants[1] = 9.81 constants[2] = 90 constants[3] = 1 constants[4] = 0.1 constants[5] = 0.1 constants[6] = 320 constants[7] = 6.667 states[0] = 1 constants[8] = 0.005 constants[9] = 70 constants[10] = 2 states[1] = 0.2 constants[11] = 0.01 constants[12] = 10 states[2] = 0.02 constants[13] = 0.004 constants[14] = 10 states[3] = 0.02 constants[15] = 0.004 constants[16] = 1 constants[17] = 1 constants[18] = 1 constants[19] = 1 states[4] = 0.3 constants[20] = 0.2 constants[21] = 200 states[5] = 0.2 constants[22] = 0.1 constants[23] = 100 constants[24] = 0.20 constants[25] = 3 constants[26] = 10 constants[27] = 20 states[6] = 0.02 constants[28] = 0.012 constants[29] = 1 states[7] = 0.1 constants[30] = 0.06 states[8] = 0.1 constants[31] = 0.06 constants[32] = 20 constants[33] = 20 constants[34] = 100 constants[35] = 600 constants[36] = 1 states[9] = 0.2 constants[37] = 0.12 constants[38] = 80 states[10] = 0.1 constants[39] = 0.06 constants[40] = 400 constants[41] = 0.20 constants[42] = 100 states[11] = 0.1 constants[43] = 0.06 constants[44] = 400 constants[45] = 0.20 constants[46] = 100 constants[47] = 100 states[12] = 0.1 constants[48] = 0.06 constants[49] = 100 constants[50] = 0.20 constants[51] = 100 constants[52] = 2400 states[13] = 0.1 constants[53] = 0.06 constants[54] = 100 constants[55] = 0.20 constants[56] = 100 constants[57] = 3000 states[14] = 0.1 constants[58] = 0.06 constants[59] = 100 constants[60] = 0.20 constants[61] = 100 constants[62] = 100 constants[63] = 4000 states[15] = 0.1 constants[64] = 0.06 constants[65] = 100 constants[66] = 0.20 constants[67] = 100 constants[68] = 1200 states[16] = 0.1 constants[69] = 0.06 constants[70] = 100 constants[71] = 0.20 constants[72] = 100 constants[73] = 2000 states[17] = 0.1 constants[74] = 0.06 constants[75] = 50 constants[76] = 0.20 constants[77] = 100 states[18] = 0.2 constants[78] = 0.06 constants[79] = 10 constants[80] = -0.50 constants[81] = 800 constants[82] = 1 states[19] = 0.1 states[20] = 0.1 states[21] = 0.1 states[22] = 0.1 constants[83] = 0.06 constants[84] = 0.06 constants[85] = 0.06 constants[86] = 0.06 constants[87] = 150 constants[88] = 150 constants[89] = 50 constants[90] = 50 constants[91] = -0.25 constants[92] = 100 constants[93] = 100 constants[94] = 2000 constants[95] = 2000 constants[96] = 1 constants[97] = 1 states[23] = 0.2 states[24] = 0.2 constants[98] = 0.12 constants[99] = 0.12 constants[100] = 60 constants[101] = 10 constants[102] = -0.50 constants[103] = 100 constants[104] = 1000 constants[105] = 1 states[25] = 0.1 states[26] = 0.1 constants[106] = 0.06 constants[107] = 0.06 constants[108] = 150 constants[109] = 60 constants[110] = -0.50 constants[111] = 100 constants[112] = 6000 constants[113] = 1 states[27] = 0.1 states[28] = 0.1 constants[114] = 0.06 constants[115] = 0.06 constants[116] = 150 constants[117] = 60 constants[118] = -0.50 constants[119] = 100 constants[120] = 6000 constants[121] = 1 states[29] = 0.1 states[30] = 0.1 constants[122] = 0.06 constants[123] = 0.06 constants[124] = 150 constants[125] = 60 constants[126] = -0.50 constants[127] = 100 constants[128] = 4000 constants[129] = 1 states[31] = 0.1 states[32] = 0.1 constants[130] = 0.06 constants[131] = 0.06 constants[132] = 150 constants[133] = 60 constants[134] = -0.50 constants[135] = 100 constants[136] = 4000 constants[137] = 1 constants[138] = (constants[2]* pi)/180.000 return (states, constants) def computeRates(voi, states, constants): rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic algebraic[2] = voi/constants[3]-floor(voi/constants[3]) algebraic[8] = custom_piecewise([less(algebraic[2] , 0.400000), (((constants[9]-constants[10])*algebraic[2])/constants[5])*exp(0.500000*(1.00000-power(algebraic[2]/constants[5], 2.00000)))+constants[10] , less(algebraic[2] , 0.600000), (((constants[9]-constants[10])*0.400000)/constants[5])*exp(0.500000*(1.00000-power(0.400000/constants[5], 2.00000)))+constants[10] , True, constants[10]]) algebraic[10] = algebraic[8]*(states[1]-constants[11]) algebraic[18] = constants[21]*(states[4]-constants[20]) algebraic[20] = custom_piecewise([greater(algebraic[10] , algebraic[18]), (algebraic[10]-algebraic[18])/constants[17] , True, 0.00000]) algebraic[14] = constants[14]*(states[3]-constants[15]) algebraic[16] = custom_piecewise([greater(algebraic[14] , algebraic[10]), (algebraic[14]-algebraic[10])/constants[18] , True, 0.00000]) rates[1] = algebraic[16]-algebraic[20] algebraic[22] = constants[23]*(states[5]-constants[22])-constants[0]*constants[1]*constants[24]*cos(constants[138]) algebraic[25] = (algebraic[18]-algebraic[22])/constants[25] rates[4] = algebraic[20]-algebraic[25] algebraic[19] = constants[32]*(states[7]-constants[30]) algebraic[21] = constants[33]*(states[8]-constants[31]) algebraic[23] = (algebraic[19]-algebraic[21])/constants[35] algebraic[24] = (algebraic[21]-algebraic[14])/constants[36] rates[8] = algebraic[23]-algebraic[24] algebraic[4] = custom_piecewise([less(algebraic[2] , 0.400000), (((constants[6]-constants[7])*algebraic[2])/constants[5])*exp(0.500000*(1.00000-power(algebraic[2]/constants[5], 2.00000)))+constants[7] , less(algebraic[2] , 0.600000), (((constants[6]-constants[7])*0.400000)/constants[5])*exp(0.500000*(1.00000-power(0.400000/constants[5], 2.00000)))+constants[7] , True, constants[7]]) algebraic[6] = algebraic[4]*(states[0]-constants[8]) algebraic[26] = constants[38]*(states[9]-constants[37]) algebraic[29] = custom_piecewise([greater(algebraic[6] , algebraic[26]), (algebraic[6]-algebraic[26])/constants[16] , True, 0.00000]) algebraic[12] = constants[12]*(states[2]-constants[13]) algebraic[17] = custom_piecewise([greater(algebraic[12] , algebraic[6]), (algebraic[12]-algebraic[6])/constants[19] , True, 0.00000]) rates[0] = algebraic[17]-algebraic[29] algebraic[27] = constants[27]*(states[6]-constants[28]) algebraic[28] = (algebraic[22]-algebraic[27])/constants[26] rates[5] = algebraic[25]-algebraic[28] algebraic[30] = (algebraic[27]-algebraic[12])/constants[29] rates[2] = algebraic[30]-algebraic[17] rates[6] = algebraic[28]-algebraic[30] algebraic[31] = (algebraic[26]-algebraic[19])/constants[34] rates[7] = algebraic[31]-algebraic[23] algebraic[34] = constants[44]*(states[11]-constants[43])-constants[0]*constants[1]*constants[45]*cos(constants[138]) algebraic[35] = (algebraic[26]-algebraic[34])/constants[46] algebraic[40] = constants[59]*(states[14]-constants[58])-constants[0]*constants[1]*constants[60]*cos(constants[138]) algebraic[42] = (algebraic[34]-algebraic[40])/constants[62] algebraic[43] = constants[65]*(states[15]-constants[64]) algebraic[44] = (algebraic[34]-algebraic[43])/constants[67] algebraic[45] = constants[70]*(states[16]-constants[69])-constants[0]*constants[1]*constants[71]*cos(constants[138]) algebraic[47] = (algebraic[34]-algebraic[45])/constants[72] rates[11] = ((algebraic[35]-algebraic[42])-algebraic[44])-algebraic[47] algebraic[32] = constants[40]*(states[10]-constants[39])-constants[0]*constants[1]*constants[41]*cos(constants[138]) algebraic[36] = constants[49]*(states[12]-constants[48])-constants[0]*constants[1]*constants[50]*cos(constants[138]) algebraic[37] = (algebraic[32]-algebraic[36])/constants[51] algebraic[48] = constants[75]*(states[17]-constants[74])-constants[0]*constants[1]*constants[76]*cos(constants[138]) algebraic[49] = (algebraic[36]-algebraic[48])/constants[52] rates[12] = algebraic[37]-algebraic[49] algebraic[38] = constants[70]*(states[13]-constants[53])-constants[0]*constants[1]*constants[55]*cos(constants[138]) algebraic[39] = (algebraic[32]-algebraic[38])/constants[56] algebraic[50] = (algebraic[38]-algebraic[48])/constants[57] rates[13] = algebraic[39]-algebraic[50] algebraic[41] = (algebraic[32]-algebraic[40])/constants[61] algebraic[51] = (algebraic[40]-algebraic[48])/constants[63] rates[14] = (algebraic[41]+algebraic[42])-algebraic[51] algebraic[52] = (algebraic[43]-algebraic[48])/constants[68] rates[15] = algebraic[44]-algebraic[52] algebraic[46] = (algebraic[26]-algebraic[45])/constants[47] algebraic[53] = (algebraic[45]-algebraic[48])/constants[73] rates[16] = (algebraic[47]+algebraic[46])-algebraic[53] algebraic[54] = constants[79]*(states[18]-constants[78])-constants[0]*constants[1]*constants[80]*cos(constants[138]) algebraic[55] = (algebraic[48]-algebraic[54])/constants[77] rates[17] = (algebraic[49]+algebraic[50]+algebraic[51]+algebraic[52]+algebraic[53])-algebraic[55] algebraic[33] = (algebraic[26]-algebraic[32])/constants[42] algebraic[56] = (algebraic[32]-algebraic[54])/constants[81] rates[10] = (((algebraic[33]-algebraic[56])-algebraic[37])-algebraic[39])-algebraic[41] algebraic[57] = (algebraic[54]-algebraic[14])/constants[82] rates[18] = (algebraic[56]+algebraic[55])-algebraic[57] algebraic[58] = constants[87]*(states[19]-constants[83])-constants[0]*constants[1]*constants[91]*cos(constants[138]) algebraic[63] = (algebraic[26]-algebraic[58])/constants[92] algebraic[60] = constants[89]*(states[21]-constants[85])-constants[0]*constants[1]*constants[91]*cos(constants[138]) algebraic[65] = (algebraic[58]-algebraic[60])/constants[94] rates[19] = algebraic[63]-algebraic[65] rates[20] = algebraic[63]-algebraic[65] algebraic[66] = (algebraic[60]-algebraic[14])/constants[96] rates[21] = algebraic[65]-algebraic[66] rates[22] = algebraic[65]-algebraic[66] algebraic[68] = constants[100]*(states[23]-constants[98])-constants[0]*constants[1]*constants[102]*cos(constants[138]) algebraic[70] = (algebraic[26]-algebraic[68])/constants[103] algebraic[69] = constants[101]*(states[24]-constants[99])-constants[0]*constants[1]*constants[102]*cos(constants[138]) algebraic[71] = (algebraic[68]-algebraic[69])/constants[104] rates[23] = algebraic[70]-algebraic[71] algebraic[72] = (algebraic[69]-algebraic[14])/constants[105] rates[24] = algebraic[71]-algebraic[72] algebraic[73] = constants[108]*(states[25]-constants[106])-constants[0]*constants[1]*constants[110]*cos(constants[138]) algebraic[75] = (algebraic[26]-algebraic[73])/constants[111] algebraic[74] = constants[109]*(states[26]-constants[107])-constants[0]*constants[1]*constants[110]*cos(constants[138]) algebraic[76] = (algebraic[73]-algebraic[74])/constants[112] rates[25] = algebraic[75]-algebraic[76] algebraic[77] = (algebraic[74]-algebraic[14])/constants[113] rates[26] = algebraic[76]-algebraic[77] algebraic[78] = constants[116]*(states[27]-constants[114])-constants[0]*constants[1]*constants[118]*cos(constants[138]) algebraic[80] = (algebraic[26]-algebraic[78])/constants[119] algebraic[79] = constants[117]*(states[28]-constants[114])-constants[0]*constants[1]*constants[118]*cos(constants[138]) algebraic[81] = (algebraic[78]-algebraic[79])/constants[120] rates[27] = algebraic[80]-algebraic[81] algebraic[82] = (algebraic[79]-algebraic[14])/constants[121] rates[28] = algebraic[81]-algebraic[82] algebraic[83] = constants[124]*(states[29]-constants[122])-constants[0]*constants[1]*constants[126]*cos(constants[138]) algebraic[85] = (algebraic[26]-algebraic[83])/constants[127] algebraic[84] = constants[125]*(states[30]-constants[122])-constants[0]*constants[1]*constants[126]*cos(constants[138]) algebraic[86] = (algebraic[83]-algebraic[84])/constants[128] rates[29] = algebraic[85]-algebraic[86] algebraic[87] = (algebraic[84]-algebraic[14])/constants[129] rates[30] = algebraic[86]-algebraic[87] algebraic[59] = constants[88]*(states[20]-constants[84])-constants[0]*constants[1]*constants[91]*cos(constants[138]) algebraic[64] = (algebraic[26]-algebraic[59])/constants[93] algebraic[88] = constants[132]*(states[31]-constants[130])-constants[0]*constants[1]*constants[134]*cos(constants[138]) algebraic[90] = (algebraic[26]-algebraic[88])/constants[135] rates[9] = ((((((((((algebraic[29]-algebraic[31])-algebraic[33])-algebraic[35])-algebraic[46])-algebraic[63])-algebraic[64])-algebraic[70])-algebraic[75])-algebraic[80])-algebraic[85])-algebraic[90] algebraic[89] = constants[133]*(states[32]-constants[130])-constants[0]*constants[1]*constants[134]*cos(constants[138]) algebraic[91] = (algebraic[88]-algebraic[89])/constants[136] rates[31] = algebraic[90]-algebraic[91] algebraic[61] = constants[90]*(states[22]-constants[86])-constants[0]*constants[1]*constants[91]*cos(constants[138]) algebraic[67] = (algebraic[61]-algebraic[14])/constants[97] algebraic[92] = (algebraic[89]-algebraic[14])/constants[137] rates[3] = (algebraic[24]+algebraic[57]+algebraic[66]+algebraic[67]+algebraic[72]+algebraic[77]+algebraic[82]+algebraic[87]+algebraic[92])-algebraic[16] rates[32] = algebraic[91]-algebraic[92] return(rates) def computeAlgebraic(constants, states, voi): algebraic = array([[0.0] * len(voi)] * sizeAlgebraic) states = array(states) voi = array(voi) algebraic[2] = voi/constants[3]-floor(voi/constants[3]) algebraic[8] = custom_piecewise([less(algebraic[2] , 0.400000), (((constants[9]-constants[10])*algebraic[2])/constants[5])*exp(0.500000*(1.00000-power(algebraic[2]/constants[5], 2.00000)))+constants[10] , less(algebraic[2] , 0.600000), (((constants[9]-constants[10])*0.400000)/constants[5])*exp(0.500000*(1.00000-power(0.400000/constants[5], 2.00000)))+constants[10] , True, constants[10]]) algebraic[10] = algebraic[8]*(states[1]-constants[11]) algebraic[18] = constants[21]*(states[4]-constants[20]) algebraic[20] = custom_piecewise([greater(algebraic[10] , algebraic[18]), (algebraic[10]-algebraic[18])/constants[17] , True, 0.00000]) algebraic[14] = constants[14]*(states[3]-constants[15]) algebraic[16] = custom_piecewise([greater(algebraic[14] , algebraic[10]), (algebraic[14]-algebraic[10])/constants[18] , True, 0.00000]) algebraic[22] = constants[23]*(states[5]-constants[22])-constants[0]*constants[1]*constants[24]*cos(constants[138]) algebraic[25] = (algebraic[18]-algebraic[22])/constants[25] algebraic[19] = constants[32]*(states[7]-constants[30]) algebraic[21] = constants[33]*(states[8]-constants[31]) algebraic[23] = (algebraic[19]-algebraic[21])/constants[35] algebraic[24] = (algebraic[21]-algebraic[14])/constants[36] algebraic[4] = custom_piecewise([less(algebraic[2] , 0.400000), (((constants[6]-constants[7])*algebraic[2])/constants[5])*exp(0.500000*(1.00000-power(algebraic[2]/constants[5], 2.00000)))+constants[7] , less(algebraic[2] , 0.600000), (((constants[6]-constants[7])*0.400000)/constants[5])*exp(0.500000*(1.00000-power(0.400000/constants[5], 2.00000)))+constants[7] , True, constants[7]]) algebraic[6] = algebraic[4]*(states[0]-constants[8]) algebraic[26] = constants[38]*(states[9]-constants[37]) algebraic[29] = custom_piecewise([greater(algebraic[6] , algebraic[26]), (algebraic[6]-algebraic[26])/constants[16] , True, 0.00000]) algebraic[12] = constants[12]*(states[2]-constants[13]) algebraic[17] = custom_piecewise([greater(algebraic[12] , algebraic[6]), (algebraic[12]-algebraic[6])/constants[19] , True, 0.00000]) algebraic[27] = constants[27]*(states[6]-constants[28]) algebraic[28] = (algebraic[22]-algebraic[27])/constants[26] algebraic[30] = (algebraic[27]-algebraic[12])/constants[29] algebraic[31] = (algebraic[26]-algebraic[19])/constants[34] algebraic[34] = constants[44]*(states[11]-constants[43])-constants[0]*constants[1]*constants[45]*cos(constants[138]) algebraic[35] = (algebraic[26]-algebraic[34])/constants[46] algebraic[40] = constants[59]*(states[14]-constants[58])-constants[0]*constants[1]*constants[60]*cos(constants[138]) algebraic[42] = (algebraic[34]-algebraic[40])/constants[62] algebraic[43] = constants[65]*(states[15]-constants[64]) algebraic[44] = (algebraic[34]-algebraic[43])/constants[67] algebraic[45] = constants[70]*(states[16]-constants[69])-constants[0]*constants[1]*constants[71]*cos(constants[138]) algebraic[47] = (algebraic[34]-algebraic[45])/constants[72] algebraic[32] = constants[40]*(states[10]-constants[39])-constants[0]*constants[1]*constants[41]*cos(constants[138]) algebraic[36] = constants[49]*(states[12]-constants[48])-constants[0]*constants[1]*constants[50]*cos(constants[138]) algebraic[37] = (algebraic[32]-algebraic[36])/constants[51] algebraic[48] = constants[75]*(states[17]-constants[74])-constants[0]*constants[1]*constants[76]*cos(constants[138]) algebraic[49] = (algebraic[36]-algebraic[48])/constants[52] algebraic[38] = constants[70]*(states[13]-constants[53])-constants[0]*constants[1]*constants[55]*cos(constants[138]) algebraic[39] = (algebraic[32]-algebraic[38])/constants[56] algebraic[50] = (algebraic[38]-algebraic[48])/constants[57] algebraic[41] = (algebraic[32]-algebraic[40])/constants[61] algebraic[51] = (algebraic[40]-algebraic[48])/constants[63] algebraic[52] = (algebraic[43]-algebraic[48])/constants[68] algebraic[46] = (algebraic[26]-algebraic[45])/constants[47] algebraic[53] = (algebraic[45]-algebraic[48])/constants[73] algebraic[54] = constants[79]*(states[18]-constants[78])-constants[0]*constants[1]*constants[80]*cos(constants[138]) algebraic[55] = (algebraic[48]-algebraic[54])/constants[77] algebraic[33] = (algebraic[26]-algebraic[32])/constants[42] algebraic[56] = (algebraic[32]-algebraic[54])/constants[81] algebraic[57] = (algebraic[54]-algebraic[14])/constants[82] algebraic[58] = constants[87]*(states[19]-constants[83])-constants[0]*constants[1]*constants[91]*cos(constants[138]) algebraic[63] = (algebraic[26]-algebraic[58])/constants[92] algebraic[60] = constants[89]*(states[21]-constants[85])-constants[0]*constants[1]*constants[91]*cos(constants[138]) algebraic[65] = (algebraic[58]-algebraic[60])/constants[94] algebraic[66] = (algebraic[60]-algebraic[14])/constants[96] algebraic[68] = constants[100]*(states[23]-constants[98])-constants[0]*constants[1]*constants[102]*cos(constants[138]) algebraic[70] = (algebraic[26]-algebraic[68])/constants[103] algebraic[69] = constants[101]*(states[24]-constants[99])-constants[0]*constants[1]*constants[102]*cos(constants[138]) algebraic[71] = (algebraic[68]-algebraic[69])/constants[104] algebraic[72] = (algebraic[69]-algebraic[14])/constants[105] algebraic[73] = constants[108]*(states[25]-constants[106])-constants[0]*constants[1]*constants[110]*cos(constants[138]) algebraic[75] = (algebraic[26]-algebraic[73])/constants[111] algebraic[74] = constants[109]*(states[26]-constants[107])-constants[0]*constants[1]*constants[110]*cos(constants[138]) algebraic[76] = (algebraic[73]-algebraic[74])/constants[112] algebraic[77] = (algebraic[74]-algebraic[14])/constants[113] algebraic[78] = constants[116]*(states[27]-constants[114])-constants[0]*constants[1]*constants[118]*cos(constants[138]) algebraic[80] = (algebraic[26]-algebraic[78])/constants[119] algebraic[79] = constants[117]*(states[28]-constants[114])-constants[0]*constants[1]*constants[118]*cos(constants[138]) algebraic[81] = (algebraic[78]-algebraic[79])/constants[120] algebraic[82] = (algebraic[79]-algebraic[14])/constants[121] algebraic[83] = constants[124]*(states[29]-constants[122])-constants[0]*constants[1]*constants[126]*cos(constants[138]) algebraic[85] = (algebraic[26]-algebraic[83])/constants[127] algebraic[84] = constants[125]*(states[30]-constants[122])-constants[0]*constants[1]*constants[126]*cos(constants[138]) algebraic[86] = (algebraic[83]-algebraic[84])/constants[128] algebraic[87] = (algebraic[84]-algebraic[14])/constants[129] algebraic[59] = constants[88]*(states[20]-constants[84])-constants[0]*constants[1]*constants[91]*cos(constants[138]) algebraic[64] = (algebraic[26]-algebraic[59])/constants[93] algebraic[88] = constants[132]*(states[31]-constants[130])-constants[0]*constants[1]*constants[134]*cos(constants[138]) algebraic[90] = (algebraic[26]-algebraic[88])/constants[135] algebraic[89] = constants[133]*(states[32]-constants[130])-constants[0]*constants[1]*constants[134]*cos(constants[138]) algebraic[91] = (algebraic[88]-algebraic[89])/constants[136] algebraic[61] = constants[90]*(states[22]-constants[86])-constants[0]*constants[1]*constants[91]*cos(constants[138]) algebraic[67] = (algebraic[61]-algebraic[14])/constants[97] algebraic[92] = (algebraic[89]-algebraic[14])/constants[137] algebraic[0] = custom_piecewise([greater(voi , 100.000) & less(voi , 101.000), 1.00000 , True, 0.00000]) algebraic[1] = (voi/constants[3]-constants[4])-floor((voi-constants[4]*constants[3])/constants[3]) algebraic[3] = states[0]+states[1]+states[2]+states[3]+states[9]+states[4] algebraic[5] = states[5]+states[6] algebraic[7] = states[19]+states[20]+states[21]+states[22] algebraic[9] = states[10]+states[11]+states[12]+states[13]+states[14]+states[15]+states[16]+states[17]+states[18]+states[23]+states[24] algebraic[11] = states[7]+states[8] algebraic[13] = states[25]+states[26]+states[27]+states[28]+states[29]+states[30]+states[31]+states[32] algebraic[15] = algebraic[3]+algebraic[5]+algebraic[7]+algebraic[9]+algebraic[11]+algebraic[13] algebraic[62] = (algebraic[59]-algebraic[61])/constants[95] 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)