# Size of variable arrays:
sizeAlgebraic = 119
sizeStates = 31
sizeConstants = 130
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 (second)"
    legend_constants[0] = "L_totmax in component b1_AR_Gs_parameters (uM)"
    legend_constants[1] = "sum_b1_AR in component b1_AR_Gs_parameters (uM)"
    legend_constants[2] = "Gs_tot in component b1_AR_Gs_parameters (uM)"
    legend_constants[3] = "Kl in component b1_AR_Gs_parameters (uM)"
    legend_constants[4] = "Kr in component b1_AR_Gs_parameters (uM)"
    legend_constants[5] = "Kc in component b1_AR_Gs_parameters (uM)"
    legend_constants[6] = "k_bar_kp in component b1_AR_Gs_parameters (per_sec)"
    legend_constants[7] = "k_bar_km in component b1_AR_Gs_parameters (per_sec)"
    legend_constants[8] = "k_p_kap in component b1_AR_Gs_parameters (per_uM_per_sec)"
    legend_constants[9] = "k_p_kam in component b1_AR_Gs_parameters (per_sec)"
    legend_constants[10] = "k_g_act in component b1_AR_Gs_parameters (per_sec)"
    legend_constants[11] = "k_hyd in component b1_AR_Gs_parameters (per_sec)"
    legend_constants[12] = "k_reassoc in component b1_AR_Gs_parameters (per_uM_per_sec)"
    legend_constants[13] = "AC_tot in component cAMP_parameters (uM)"
    legend_constants[14] = "ATP in component cAMP_parameters (uM)"
    legend_constants[15] = "PDE_tot in component cAMP_parameters (uM)"
    legend_constants[16] = "IBMX_tot in component cAMP_parameters (uM)"
    legend_constants[17] = "Fsk_tot in component cAMP_parameters (uM)"
    legend_constants[18] = "k_ac_basal in component cAMP_parameters (per_sec)"
    legend_constants[19] = "k_ac_gsa in component cAMP_parameters (per_sec)"
    legend_constants[20] = "k_ac_fsk in component cAMP_parameters (per_sec)"
    legend_constants[21] = "k_pde in component cAMP_parameters (per_sec)"
    legend_constants[22] = "Km_basal in component cAMP_parameters (uM)"
    legend_constants[23] = "Km_gsa in component cAMP_parameters (uM)"
    legend_constants[24] = "Km_fsk in component cAMP_parameters (uM)"
    legend_constants[25] = "Km_pde in component cAMP_parameters (uM)"
    legend_constants[26] = "K_gsa in component cAMP_parameters (uM)"
    legend_constants[27] = "K_fsk in component cAMP_parameters (uM)"
    legend_constants[28] = "Ki_ibmx in component cAMP_parameters (uM)"
    legend_constants[29] = "PKAI_tot in component PKA_parameters (uM)"
    legend_constants[30] = "PKAII_tot in component PKA_parameters (uM)"
    legend_constants[31] = "PKI_tot in component PKA_parameters (uM)"
    legend_constants[32] = "K_a in component PKA_parameters (uM)"
    legend_constants[33] = "K_b in component PKA_parameters (uM)"
    legend_constants[34] = "K_d in component PKA_parameters (uM)"
    legend_constants[35] = "Ki_pki in component PKA_parameters (uM)"
    legend_constants[36] = "PLB_tot in component PLB_parameters (uM)"
    legend_constants[37] = "PP1_tot in component PLB_parameters (uM)"
    legend_constants[38] = "Inhib1_tot in component PLB_parameters (uM)"
    legend_constants[39] = "k_pka_plb in component PLB_parameters (per_sec)"
    legend_constants[40] = "Km_pka_plb in component PLB_parameters (uM)"
    legend_constants[41] = "k_pp1_plb in component PLB_parameters (per_sec)"
    legend_constants[42] = "Km_pp1_plb in component PLB_parameters (uM)"
    legend_constants[43] = "k_pka_i1 in component PLB_parameters (per_sec)"
    legend_constants[44] = "Km_pka_i1 in component PLB_parameters (uM)"
    legend_constants[45] = "Vmax_pp2a_i1 in component PLB_parameters (uM_per_sec)"
    legend_constants[46] = "Km_pp2a_i1 in component PLB_parameters (uM)"
    legend_constants[47] = "Ki_inhib1 in component PLB_parameters (uM)"
    legend_constants[48] = "epsilon in component LCC_parameters (dimensionless)"
    legend_constants[49] = "LCC_tot in component LCC_parameters (uM)"
    legend_constants[50] = "PP1_lcc_tot in component LCC_parameters (uM)"
    legend_constants[51] = "PP2A_lcc_tot in component LCC_parameters (uM)"
    legend_constants[52] = "k_pka_lcc in component LCC_parameters (per_sec)"
    legend_constants[53] = "Km_pka_lcc in component LCC_parameters (uM)"
    legend_constants[54] = "k_pp1_lcc in component LCC_parameters (per_sec)"
    legend_constants[55] = "Km_pp1_lcc in component LCC_parameters (uM)"
    legend_constants[56] = "k_pp2a_lcc in component LCC_parameters (per_sec)"
    legend_constants[57] = "Km_pp2a_lcc in component LCC_parameters (uM)"
    legend_constants[58] = "V_myo in component EC_Coupling_Parameters (uL)"
    legend_constants[59] = "V_nsr in component EC_Coupling_Parameters (uL)"
    legend_constants[60] = "V_jsr in component EC_Coupling_Parameters (uL)"
    legend_constants[61] = "A_Cap in component EC_Coupling_Parameters (cm2)"
    legend_constants[62] = "Temp in component EC_Coupling_Parameters (kelvin)"
    legend_constants[63] = "Na_ext in component EC_Coupling_Parameters (mM)"
    legend_constants[64] = "K_ext in component EC_Coupling_Parameters (mM)"
    legend_constants[65] = "Ca_ext in component EC_Coupling_Parameters (mM)"
    legend_constants[66] = "G_Na in component EC_Coupling_Parameters (mS_per_uF)"
    legend_constants[67] = "G_to in component EC_Coupling_Parameters (mS_per_uF)"
    legend_constants[68] = "G_ss in component EC_Coupling_Parameters (mS_per_uF)"
    legend_constants[69] = "G_Ki_bar in component EC_Coupling_Parameters (mS_per_uF)"
    legend_constants[70] = "G_Kp in component EC_Coupling_Parameters (mS_per_uF)"
    legend_constants[71] = "f in component EC_Coupling_Parameters (per_sec)"
    legend_constants[72] = "g in component EC_Coupling_Parameters (per_sec)"
    legend_constants[73] = "gamma_o in component EC_Coupling_Parameters (per_mM_per_sec)"
    legend_constants[74] = "omega in component EC_Coupling_Parameters (per_sec)"
    legend_constants[75] = "p_Ca in component EC_Coupling_Parameters (cm_per_sec)"
    legend_constants[76] = "p_K in component EC_Coupling_Parameters (cm_per_sec)"
    legend_constants[77] = "N_lcc in component EC_Coupling_Parameters (dimensionless)"
    legend_constants[78] = "I_Ca05 in component EC_Coupling_Parameters (uA_per_uF)"
    legend_constants[79] = "k_NaCa in component EC_Coupling_Parameters (uA_per_uF)"
    legend_constants[80] = "Km_Na in component EC_Coupling_Parameters (mM)"
    legend_constants[81] = "Km_Ca in component EC_Coupling_Parameters (mM)"
    legend_constants[82] = "k_sat in component EC_Coupling_Parameters (dimensionless)"
    legend_constants[83] = "eta in component EC_Coupling_Parameters (dimensionless)"
    legend_constants[84] = "I_bar_NaK in component EC_Coupling_Parameters (uA_per_uF)"
    legend_constants[85] = "Km_Nai in component EC_Coupling_Parameters (mM)"
    legend_constants[86] = "Km_Ko in component EC_Coupling_Parameters (mM)"
    legend_constants[87] = "I_bar_PCa in component EC_Coupling_Parameters (uA_per_uF)"
    legend_constants[88] = "Km_PCa in component EC_Coupling_Parameters (mM)"
    legend_constants[89] = "G_CaB in component EC_Coupling_Parameters (uA_per_uF)"
    legend_constants[90] = "G_NaB in component EC_Coupling_Parameters (uA_per_uF)"
    legend_constants[91] = "Pns in component EC_Coupling_Parameters (dimensionless)"
    legend_constants[92] = "Km_NS in component EC_Coupling_Parameters (mM)"
    legend_constants[93] = "I_up_bar in component EC_Coupling_Parameters (mM_per_sec)"
    legend_constants[94] = "Km_up0 in component EC_Coupling_Parameters (mM)"
    legend_constants[95] = "NSR_bar in component EC_Coupling_Parameters (mM)"
    legend_constants[96] = "tau_on in component EC_Coupling_Parameters (second)"
    legend_constants[97] = "tau_off in component EC_Coupling_Parameters (second)"
    legend_constants[98] = "G_max_rel in component EC_Coupling_Parameters (mM_per_sec)"
    legend_constants[99] = "d_Cai_th in component EC_Coupling_Parameters (mM)"
    legend_constants[100] = "Km_rel in component EC_Coupling_Parameters (mM)"
    legend_constants[101] = "CSQN_th in component EC_Coupling_Parameters (mM)"
    legend_constants[102] = "CSQN_bar in component EC_Coupling_Parameters (mM)"
    legend_constants[103] = "Km_CSQN in component EC_Coupling_Parameters (mM)"
    legend_constants[104] = "tau_tr in component EC_Coupling_Parameters (second)"
    legend_constants[105] = "TRPN_bar in component EC_Coupling_Parameters (mM)"
    legend_constants[106] = "CMDN_bar in component EC_Coupling_Parameters (mM)"
    legend_constants[107] = "INDO_bar in component EC_Coupling_Parameters (mM)"
    legend_constants[108] = "Km_TRPN in component EC_Coupling_Parameters (mM)"
    legend_constants[109] = "Km_CMDN in component EC_Coupling_Parameters (mM)"
    legend_constants[110] = "Km_INDO in component EC_Coupling_Parameters (mM)"
    legend_algebraic[62] = "LR in component b1_AR_module (uM)"
    legend_algebraic[63] = "LRG in component b1_AR_module (uM)"
    legend_algebraic[64] = "RG in component b1_AR_module (uM)"
    legend_algebraic[80] = "BARK_DESENS in component b1_AR_module (uM_per_sec)"
    legend_algebraic[11] = "BARK_RESENS in component b1_AR_module (uM_per_sec)"
    legend_algebraic[88] = "PKA_DESENS in component b1_AR_module (uM_per_sec)"
    legend_algebraic[23] = "PKA_RESENS in component b1_AR_module (uM_per_sec)"
    legend_algebraic[81] = "G_ACT in component b1_AR_module (uM_per_sec)"
    legend_algebraic[25] = "HYD in component b1_AR_module (uM_per_sec)"
    legend_algebraic[27] = "REASSOC in component b1_AR_module (per_sec)"
    legend_algebraic[65] = "L in component b1_AR_module (uM)"
    legend_algebraic[66] = "R in component b1_AR_module (uM)"
    legend_algebraic[67] = "Gs in component b1_AR_module (uM)"
    legend_states[0] = "b1_AR_tot in component b1_AR_module (uM)"
    legend_states[1] = "b1_AR_d in component b1_AR_module (uM)"
    legend_states[2] = "b1_AR_p in component b1_AR_module (uM)"
    legend_states[3] = "Gs_agtp_tot in component b1_AR_module (uM)"
    legend_states[4] = "Gs_agdp in component b1_AR_module (uM)"
    legend_states[5] = "Gs_bg in component b1_AR_module (uM)"
    legend_algebraic[68] = "PKAC_I in component PKA_module (uM)"
    legend_algebraic[69] = "cAMP in component PKA_module (uM)"
    legend_algebraic[42] = "Gsa_GTP in component cAMP_module (uM)"
    legend_algebraic[50] = "Fsk in component cAMP_module (uM)"
    legend_algebraic[43] = "AC in component cAMP_module (uM)"
    legend_constants[127] = "PDE in component cAMP_module (uM)"
    legend_constants[128] = "IBMX in component cAMP_module (uM)"
    legend_states[6] = "cAMP_tot in component cAMP_module (uM)"
    legend_algebraic[44] = "Gsa_GTP_AC in component cAMP_module (uM)"
    legend_algebraic[51] = "Fsk_AC in component cAMP_module (uM)"
    legend_algebraic[48] = "AC_ACT_GSA in component cAMP_module (uM)"
    legend_algebraic[46] = "AC_ACT_BASAL in component cAMP_module (uM)"
    legend_algebraic[53] = "AC_ACT_FSK in component cAMP_module (uM)"
    legend_algebraic[82] = "PDE_ACT in component cAMP_module (uM)"
    legend_constants[129] = "PDE_IBMX in component cAMP_module (uM)"
    legend_algebraic[70] = "PKI in component PKA_module (uM)"
    legend_algebraic[71] = "A2RC_I in component PKA_module (uM)"
    legend_algebraic[72] = "A2R_I in component PKA_module (uM)"
    legend_algebraic[73] = "A2RC_II in component PKA_module (uM)"
    legend_algebraic[74] = "A2R_II in component PKA_module (uM)"
    legend_algebraic[75] = "ARC_I in component PKA_module (uM)"
    legend_algebraic[76] = "ARC_II in component PKA_module (uM)"
    legend_algebraic[77] = "PKA_temp in component PKA_module (uM)"
    legend_algebraic[78] = "PKAC_II in component PKA_module (uM)"
    legend_algebraic[29] = "PLB in component PLB_module (uM)"
    legend_algebraic[83] = "PLB_PHOSPH in component PLB_module (uM_per_sec)"
    legend_algebraic[60] = "PLB_DEPHOSPH in component PLB_module (uM_per_sec)"
    legend_algebraic[30] = "Inhib1 in component PLB_module (uM)"
    legend_algebraic[55] = "Inhib1p_PP1 in component PLB_module (uM)"
    legend_algebraic[84] = "Inhib1_PHOSPH in component PLB_module (uM_per_sec)"
    legend_algebraic[32] = "Inhib1_DEPHOSPH in component PLB_module (uM_per_sec)"
    legend_states[7] = "PLB_p in component PLB_module (uM)"
    legend_states[8] = "Inhib1p_tot in component PLB_module (uM)"
    legend_algebraic[56] = "Inhib1p in component PLB_module (uM)"
    legend_algebraic[57] = "PP1 in component PLB_module (uM)"
    legend_algebraic[0] = "frac_PLB_p in component PLB_module (dimensionless)"
    legend_algebraic[31] = "frac_PLB in component PLB_module (dimensionless)"
    legend_constants[119] = "frac_PLB_o in component PLB_module (dimensionless)"
    legend_algebraic[34] = "LCCa in component LCC_module (uM)"
    legend_algebraic[85] = "LCCa_PHOSPH in component LCC_module (uM_per_sec)"
    legend_algebraic[36] = "LCCa_DEPHOSPH in component LCC_module (uM_per_sec)"
    legend_algebraic[38] = "LCCb in component LCC_module (uM)"
    legend_algebraic[86] = "LCCb_PHOSPH in component LCC_module (uM_per_sec)"
    legend_algebraic[40] = "LCCb_DEPHOSPH in component LCC_module (uM_per_sec)"
    legend_states[9] = "LCCa_p in component LCC_module (uM)"
    legend_states[10] = "LCCb_p in component LCC_module (uM)"
    legend_algebraic[1] = "frac_LCCa_p in component LCC_module (dimensionless)"
    legend_constants[120] = "frac_LCCa_po in component LCC_module (dimensionless)"
    legend_algebraic[33] = "frac_LCCb_p in component LCC_module (dimensionless)"
    legend_constants[121] = "frac_LCCb_po in component LCC_module (dimensionless)"
    legend_algebraic[35] = "E_Na in component Nernst_Potentials (mV)"
    legend_algebraic[37] = "E_K in component Nernst_Potentials (mV)"
    legend_algebraic[39] = "E_Ca in component Nernst_Potentials (mV)"
    legend_constants[122] = "E_Cl in component Nernst_Potentials (mV)"
    legend_constants[111] = "R in component Nernst_Potentials (joules_per_mole_kelvin)"
    legend_constants[112] = "Frdy in component Nernst_Potentials (coulombs_per_mole)"
    legend_constants[123] = "FoRT in component Nernst_Potentials (per_mV)"
    legend_constants[113] = "z_Na in component Nernst_Potentials (dimensionless)"
    legend_constants[114] = "z_K in component Nernst_Potentials (dimensionless)"
    legend_constants[115] = "z_Ca in component Nernst_Potentials (dimensionless)"
    legend_states[11] = "Na_i in component Ion_Concentrations_and_Membrane_Potential (mM)"
    legend_states[12] = "K_i in component Ion_Concentrations_and_Membrane_Potential (mM)"
    legend_states[13] = "Ca_i in component Ion_Concentrations_and_Membrane_Potential (mM)"
    legend_algebraic[2] = "am in component Fast_Na_Current (per_sec)"
    legend_algebraic[13] = "bm in component Fast_Na_Current (per_sec)"
    legend_algebraic[3] = "ah in component Fast_Na_Current (per_sec)"
    legend_algebraic[4] = "aj in component Fast_Na_Current (per_sec)"
    legend_algebraic[14] = "bh in component Fast_Na_Current (per_sec)"
    legend_algebraic[15] = "bj in component Fast_Na_Current (per_sec)"
    legend_states[14] = "m in component Fast_Na_Current (dimensionless)"
    legend_states[15] = "h in component Fast_Na_Current (dimensionless)"
    legend_states[16] = "j in component Fast_Na_Current (dimensionless)"
    legend_algebraic[41] = "I_Na in component Fast_Na_Current (uA_per_uF)"
    legend_states[17] = "V_m in component Ion_Concentrations_and_Membrane_Potential (mV)"
    legend_algebraic[5] = "a_lcc in component L_Type_Calcium_Current (per_sec)"
    legend_algebraic[16] = "b_lcc in component L_Type_Calcium_Current (per_sec)"
    legend_algebraic[18] = "f_lcc in component L_Type_Calcium_Current (dimensionless)"
    legend_algebraic[6] = "y_lcc_inf in component L_Type_Calcium_Current (dimensionless)"
    legend_algebraic[17] = "tau_y_lcc in component L_Type_Calcium_Current (second)"
    legend_algebraic[24] = "gamma in component L_Type_Calcium_Current (per_mM_per_sec)"
    legend_algebraic[26] = "v_gamma in component L_Type_Calcium_Current (per_mM_per_sec)"
    legend_algebraic[28] = "v_omega in component L_Type_Calcium_Current (dimensionless)"
    legend_states[18] = "v in component L_Type_Calcium_Current (dimensionless)"
    legend_states[19] = "w in component L_Type_Calcium_Current (dimensionless)"
    legend_states[20] = "x in component L_Type_Calcium_Current (dimensionless)"
    legend_states[21] = "y in component L_Type_Calcium_Current (dimensionless)"
    legend_states[22] = "z in component L_Type_Calcium_Current (dimensionless)"
    legend_algebraic[45] = "i_bar_Ca in component L_Type_Calcium_Current (dimensionless)"
    legend_algebraic[47] = "i_bar_K in component L_Type_Calcium_Current (dimensionless)"
    legend_algebraic[49] = "f_avail in component L_Type_Calcium_Current (dimensionless)"
    legend_algebraic[52] = "I_Ca in component L_Type_Calcium_Current (uA_per_uF)"
    legend_algebraic[54] = "I_CaK in component L_Type_Calcium_Current (uA_per_uF)"
    legend_algebraic[58] = "I_Ca_tot in component L_Type_Calcium_Current (uA_per_uF)"
    legend_algebraic[7] = "r_toss in component Transient_Outward_K_Current (dimensionless)"
    legend_algebraic[8] = "s_toss in component Transient_Outward_K_Current (dimensionless)"
    legend_algebraic[19] = "tau_r_to in component Transient_Outward_K_Current (second)"
    legend_algebraic[20] = "tau_s_to in component Transient_Outward_K_Current (second)"
    legend_algebraic[21] = "tau_ss_to in component Transient_Outward_K_Current (second)"
    legend_states[23] = "r_to in component Transient_Outward_K_Current (dimensionless)"
    legend_states[24] = "s_to in component Transient_Outward_K_Current (dimensionless)"
    legend_states[25] = "ss_to in component Transient_Outward_K_Current (dimensionless)"
    legend_algebraic[59] = "I_to in component Transient_Outward_K_Current (uA_per_uF)"
    legend_algebraic[9] = "r_ss_inf in component Steady_State_K_Current (dimensionless)"
    legend_algebraic[22] = "tau_r_ss in component Steady_State_K_Current (second)"
    legend_algebraic[10] = "s_ss_inf in component Steady_State_K_Current (dimensionless)"
    legend_constants[124] = "tau_s_ss in component Steady_State_K_Current (second)"
    legend_states[26] = "r_ss in component Steady_State_K_Current (dimensionless)"
    legend_states[27] = "s_ss in component Steady_State_K_Current (dimensionless)"
    legend_algebraic[61] = "I_ss in component Steady_State_K_Current (uA_per_uF)"
    legend_algebraic[79] = "a_Ki in component Time_Independent_K_Current (dimensionless)"
    legend_algebraic[87] = "b_Ki in component Time_Independent_K_Current (dimensionless)"
    legend_algebraic[89] = "Ki_ss in component Time_Independent_K_Current (dimensionless)"
    legend_algebraic[90] = "I_Ki in component Time_Independent_K_Current (uA_per_uF)"
    legend_algebraic[91] = "Kp in component Plateau_K_Current (dimensionless)"
    legend_algebraic[92] = "I_Kp in component Plateau_K_Current (uA_per_uF)"
    legend_algebraic[93] = "s4 in component Na_Ca_Exchanger_Current (dimensionless)"
    legend_algebraic[94] = "s5 in component Na_Ca_Exchanger_Current (dimensionless)"
    legend_algebraic[95] = "I_NCX in component Na_Ca_Exchanger_Current (uA_per_uF)"
    legend_constants[125] = "sigma in component Na_K_Pump_Current (dimensionless)"
    legend_algebraic[96] = "f_NaK in component Na_K_Pump_Current (dimensionless)"
    legend_algebraic[97] = "I_NaK in component Na_K_Pump_Current (uA_per_uF)"
    legend_algebraic[98] = "I_PCa in component Sarcolemmal_Ca_Pump_Current (uA_per_uF)"
    legend_algebraic[99] = "I_CaB in component Ca_Background_Current (uA_per_uF)"
    legend_algebraic[100] = "I_NaB in component Na_Background_Current (uA_per_uF)"
    legend_algebraic[101] = "I_Na_tot in component Total_Membrane_Currents (uA_per_uF)"
    legend_algebraic[102] = "I_K_tot in component Total_Membrane_Currents (uA_per_uF)"
    legend_algebraic[103] = "I_Ca_tot in component Total_Membrane_Currents (uA_per_uF)"
    legend_algebraic[104] = "t_rel in component Calcium_Induced_Calcium_Release (second)"
    legend_algebraic[105] = "ryr_on in component Calcium_Induced_Calcium_Release (dimensionless)"
    legend_algebraic[106] = "ryr_off in component Calcium_Induced_Calcium_Release (dimensionless)"
    legend_algebraic[107] = "g_rel in component Calcium_Induced_Calcium_Release (per_sec)"
    legend_algebraic[108] = "I_rel in component Calcium_Induced_Calcium_Release (mM_per_sec)"
    legend_states[28] = "Ca_jsr in component Other_SR_Fluxes_and_Concentrations (mM)"
    legend_states[29] = "trel in component Ion_Concentrations_and_Membrane_Potential (second)"
    legend_algebraic[109] = "Km_up in component Other_SR_Fluxes_and_Concentrations (mM)"
    legend_algebraic[110] = "I_up in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)"
    legend_algebraic[111] = "I_leak in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)"
    legend_algebraic[112] = "I_tr in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)"
    legend_algebraic[114] = "B_jsr in component Other_SR_Fluxes_and_Concentrations (dimensionless)"
    legend_states[30] = "Ca_nsr in component Other_SR_Fluxes_and_Concentrations (mM)"
    legend_algebraic[116] = "SR_content in component Other_SR_Fluxes_and_Concentrations (uM)"
    legend_algebraic[113] = "b_trpn in component Cytoplasmic_Calcium_Buffering (dimensionless)"
    legend_algebraic[115] = "b_cmdn in component Cytoplasmic_Calcium_Buffering (dimensionless)"
    legend_algebraic[117] = "b_indo in component Cytoplasmic_Calcium_Buffering (dimensionless)"
    legend_algebraic[118] = "B_myo in component Cytoplasmic_Calcium_Buffering (dimensionless)"
    legend_constants[126] = "I_app in component Ion_Concentrations_and_Membrane_Potential (uA_per_uF)"
    legend_constants[116] = "V_hold in component Ion_Concentrations_and_Membrane_Potential (mV)"
    legend_constants[117] = "V_test in component Ion_Concentrations_and_Membrane_Potential (mV)"
    legend_algebraic[12] = "V_clamp in component Ion_Concentrations_and_Membrane_Potential (mV)"
    legend_constants[118] = "R_clamp in component Ion_Concentrations_and_Membrane_Potential (dimensionless)"
    legend_rates[0] = "d/dt b1_AR_tot in component b1_AR_module (uM)"
    legend_rates[1] = "d/dt b1_AR_d in component b1_AR_module (uM)"
    legend_rates[2] = "d/dt b1_AR_p in component b1_AR_module (uM)"
    legend_rates[3] = "d/dt Gs_agtp_tot in component b1_AR_module (uM)"
    legend_rates[4] = "d/dt Gs_agdp in component b1_AR_module (uM)"
    legend_rates[5] = "d/dt Gs_bg in component b1_AR_module (uM)"
    legend_rates[6] = "d/dt cAMP_tot in component cAMP_module (uM)"
    legend_rates[7] = "d/dt PLB_p in component PLB_module (uM)"
    legend_rates[8] = "d/dt Inhib1p_tot in component PLB_module (uM)"
    legend_rates[9] = "d/dt LCCa_p in component LCC_module (uM)"
    legend_rates[10] = "d/dt LCCb_p in component LCC_module (uM)"
    legend_rates[14] = "d/dt m in component Fast_Na_Current (dimensionless)"
    legend_rates[15] = "d/dt h in component Fast_Na_Current (dimensionless)"
    legend_rates[16] = "d/dt j in component Fast_Na_Current (dimensionless)"
    legend_rates[18] = "d/dt v in component L_Type_Calcium_Current (dimensionless)"
    legend_rates[19] = "d/dt w in component L_Type_Calcium_Current (dimensionless)"
    legend_rates[20] = "d/dt x in component L_Type_Calcium_Current (dimensionless)"
    legend_rates[21] = "d/dt y in component L_Type_Calcium_Current (dimensionless)"
    legend_rates[22] = "d/dt z in component L_Type_Calcium_Current (dimensionless)"
    legend_rates[23] = "d/dt r_to in component Transient_Outward_K_Current (dimensionless)"
    legend_rates[24] = "d/dt s_to in component Transient_Outward_K_Current (dimensionless)"
    legend_rates[25] = "d/dt ss_to in component Transient_Outward_K_Current (dimensionless)"
    legend_rates[26] = "d/dt r_ss in component Steady_State_K_Current (dimensionless)"
    legend_rates[27] = "d/dt s_ss in component Steady_State_K_Current (dimensionless)"
    legend_rates[30] = "d/dt Ca_nsr in component Other_SR_Fluxes_and_Concentrations (mM)"
    legend_rates[28] = "d/dt Ca_jsr in component Other_SR_Fluxes_and_Concentrations (mM)"
    legend_rates[11] = "d/dt Na_i in component Ion_Concentrations_and_Membrane_Potential (mM)"
    legend_rates[12] = "d/dt K_i in component Ion_Concentrations_and_Membrane_Potential (mM)"
    legend_rates[13] = "d/dt Ca_i in component Ion_Concentrations_and_Membrane_Potential (mM)"
    legend_rates[17] = "d/dt V_m in component Ion_Concentrations_and_Membrane_Potential (mV)"
    legend_rates[29] = "d/dt trel in component Ion_Concentrations_and_Membrane_Potential (second)"
    return (legend_states, legend_algebraic, legend_voi, legend_constants)

def initConsts():
    constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
    constants[0] = 1
    constants[1] = 0.0132
    constants[2] = 3.83
    constants[3] = 0.285
    constants[4] = 0.062
    constants[5] = 33
    constants[6] = 1.1e-3
    constants[7] = 2.2e-3
    constants[8] = 3.6e-3
    constants[9] = 2.2e-3
    constants[10] = 16
    constants[11] = 0.8
    constants[12] = 1.21e3
    constants[13] = 49.7e-3
    constants[14] = 5e3
    constants[15] = 38.9e-3
    constants[16] = 0
    constants[17] = 0
    constants[18] = 0.2
    constants[19] = 8.5
    constants[20] = 7.3
    constants[21] = 5
    constants[22] = 1.03e3
    constants[23] = 315
    constants[24] = 860
    constants[25] = 1.3
    constants[26] = 0.4
    constants[27] = 44
    constants[28] = 30
    constants[29] = 0.59
    constants[30] = 0.025
    constants[31] = 0.18
    constants[32] = 9.14
    constants[33] = 1.64
    constants[34] = 4.375
    constants[35] = 0.2e-3
    constants[36] = 106
    constants[37] = 0.89
    constants[38] = 0.3
    constants[39] = 54
    constants[40] = 21
    constants[41] = 8.5
    constants[42] = 7
    constants[43] = 60
    constants[44] = 1
    constants[45] = 14
    constants[46] = 1
    constants[47] = 1e-3
    constants[48] = 10
    constants[49] = 25e-3
    constants[50] = 25e-3
    constants[51] = 25e-3
    constants[52] = 54
    constants[53] = 21
    constants[54] = 8.52
    constants[55] = 3
    constants[56] = 10.1
    constants[57] = 3
    constants[58] = 20.8e-6
    constants[59] = 9.88e-7
    constants[60] = 9.3e-8
    constants[61] = 1.534e-4
    constants[62] = 310
    constants[63] = 140
    constants[64] = 5.4
    constants[65] = 1.8
    constants[66] = 8
    constants[67] = 0.35
    constants[68] = 0.07
    constants[69] = 0.24
    constants[70] = 0.008
    constants[71] = 300
    constants[72] = 2e3
    constants[73] = 5187.5
    constants[74] = 10
    constants[75] = 1.7469e-8
    constants[76] = 3.234e-11
    constants[77] = 3e5
    constants[78] = -0.458
    constants[79] = 1483
    constants[80] = 87.5
    constants[81] = 1.38
    constants[82] = 0.1
    constants[83] = 0.35
    constants[84] = 1.1
    constants[85] = 10
    constants[86] = 1.5
    constants[87] = 1.15
    constants[88] = 0.5e-3
    constants[89] = 2.8e-3
    constants[90] = 1.18e-3
    constants[91] = 0
    constants[92] = 1.2e-3
    constants[93] = 4.7
    constants[94] = 3e-4
    constants[95] = 15
    constants[96] = 2e-3
    constants[97] = 2e-3
    constants[98] = 60e3
    constants[99] = 0.18e-3
    constants[100] = 0.8e-3
    constants[101] = 8.75
    constants[102] = 15
    constants[103] = 0.8
    constants[104] = 5.7e-4
    constants[105] = 0.07
    constants[106] = 0.05
    constants[107] = 0.07
    constants[108] = 0.5128e-3
    constants[109] = 2.38e-3
    constants[110] = 8.44e-4
    states[0] = 0.01205
    states[1] = 0
    states[2] = 1.154e-3
    states[3] = 0.02505
    states[4] = 6.446e-4
    states[5] = 0.02569
    states[6] = 0.8453
    states[7] = 4.105
    states[8] = 0.0526
    states[9] = 5.103e-3
    states[10] = 5.841e-3
    constants[111] = 8314
    constants[112] = 96485
    constants[113] = 1
    constants[114] = 1
    constants[115] = 2
    states[11] = 16
    states[12] = 145
    states[13] = 1.58e-4
    states[14] = 1.4e-3
    states[15] = 0.99
    states[16] = 0.99
    states[17] = -85.66
    states[18] = 0
    states[19] = 0
    states[20] = 0.13
    states[21] = 0.96
    states[22] = 0.92
    states[23] = 1.4e-3
    states[24] = 1
    states[25] = 0.613
    states[26] = 198e-3
    states[27] = 0.43
    states[28] = 1.92
    states[29] = 0.9
    states[30] = 1.92
    constants[116] = -40
    constants[117] = -10
    constants[118] = 0.02
    constants[119] = 0.961300
    constants[120] = 0.204100
    constants[121] = 0.233600
    constants[122] = -40.0000
    constants[123] = (constants[112]/constants[111])/constants[62]
    constants[124] = 2.10000
    constants[125] = (exp(constants[63]/67.3000)-1.00000)/7.00000
    constants[126] = 0.00000
    rootfind_0(voi, constants, rates, states, algebraic)
    return (states, constants)

def computeRates(voi, states, constants):
    rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
    algebraic[10] = 1.00000/(1.00000+exp((states[17]+87.5000)/10.3000))
    rates[27] = (algebraic[10]-states[27])/constants[124]
    algebraic[2] = (0.320000*(states[17]+47.1300))/(1.00000-exp(-0.100000*(states[17]+47.1300)))
    algebraic[13] = 0.0800000*exp(-states[17]/11.0000)
    rates[14] = 1000.00*(algebraic[2]*(1.00000-states[14])-algebraic[13]*states[14])
    algebraic[3] = custom_piecewise([greater_equal(states[17] , -40.0000), 0.00000 , True, 0.135000*exp((80.0000+states[17])/-6.80000)])
    algebraic[14] = custom_piecewise([greater_equal(states[17] , -40.0000), 1.00000/(0.130000*(1.00000+exp(-(states[17]+10.6600)/11.1000))) , True, 3.56000*exp(0.0790000*states[17])+310000.*exp(0.350000*states[17])])
    rates[15] = 1000.00*(algebraic[3]*(1.00000-states[15])-algebraic[14]*states[15])
    algebraic[4] = custom_piecewise([greater_equal(states[17] , -40.0000), 0.00000 , True, ((-127140.*exp(0.244400*states[17])-3.47400e-05*exp(-0.0439100*states[17]))*(states[17]+37.7800))/(1.00000+exp(0.311000*(states[17]+79.2300)))])
    algebraic[15] = custom_piecewise([greater_equal(states[17] , -40.0000), (0.300000*exp(-2.57500e-07*states[17]))/(1.00000+exp(-0.100000*(states[17]+32.0000))) , True, (0.121200*exp(-0.0105200*states[17]))/(1.00000+exp(-0.137800*(states[17]+40.1400)))])
    rates[16] = 1000.00*(algebraic[4]*(1.00000-states[16])-algebraic[15]*states[16])
    algebraic[5] = 400.000*exp((states[17]+2.00000)/10.0000)
    algebraic[16] = 50.0000*exp((-1.00000*(states[17]+2.00000))/13.0000)
    rates[18] = algebraic[5]*(1.00000-states[18])-algebraic[16]*states[18]
    rates[19] = 2.00000*algebraic[5]*(1.00000-states[19])-(algebraic[16]/2.00000)*states[19]
    algebraic[1] = states[9]/constants[49]
    algebraic[18] = constants[71]*((0.375000*algebraic[1])/constants[120]+0.625000)
    rates[20] = algebraic[18]*(1.00000-states[20])-constants[72]*states[20]
    algebraic[6] = 1.00000/(1.00000+exp((states[17]+55.0000)/7.50000))+0.100000/(1.00000+exp((-states[17]+21.0000)/6.00000))
    algebraic[17] = 0.0200000+0.300000/(1.00000+exp((states[17]+30.0000)/9.50000))
    rates[21] = (algebraic[6]-states[21])/algebraic[17]
    algebraic[7] = 1.00000/(1.00000+exp((states[17]+10.6000)/-11.4200))
    algebraic[19] = 1.00000/(45.1600*exp(0.0357700*(states[17]+50.0000))+98.9000*exp(-0.100000*(states[17]+38.0000)))
    rates[23] = (algebraic[7]-states[23])/algebraic[19]
    algebraic[8] = 1.00000/(1.00000+exp((states[17]+43.5000)/6.88410))
    algebraic[20] = 0.350000*exp(-1.00000*(power((states[17]+70.0000)/15.0000, 2.00000)))+0.0350000
    rates[24] = (algebraic[8]-states[24])/algebraic[20]
    algebraic[21] = 3.70000*exp(-1.00000*(power((states[17]+70.0000)/30.0000, 2.00000)))+0.0350000
    rates[25] = (algebraic[8]-states[25])/algebraic[21]
    algebraic[9] = 1.00000/(1.00000+exp((states[17]+11.5000)/-11.8200))
    algebraic[22] = 10.0000/(45.1600*exp(0.0357700*(states[17]+50.0000))+98.9000*exp(-0.100000*(states[17]+38.0000)))
    rates[26] = (algebraic[9]-states[26])/algebraic[22]
    algebraic[25] = constants[11]*states[3]
    algebraic[27] = constants[12]*states[4]*states[5]
    rates[4] = algebraic[25]-algebraic[27]
    algebraic[24] = constants[73]*states[13]
    algebraic[26] = algebraic[24]*(power(1.00000-states[18], 4.00000)+2.00000*states[18]*(power(1.00000-states[18], 3.00000))+4.00000*(power(states[18], 2.00000))*(power(1.00000-states[18], 2.00000))+8.00000*(power(states[18], 3.00000))*(1.00000-states[18])+16.0000*(power(states[18], 4.00000))*(1.00000-algebraic[18]/constants[72]))
    algebraic[28] = constants[74]*(power(1.00000-states[19], 4.00000)+0.500000*states[19]*(power(1.00000-states[19], 3.00000))+0.250000*(power(states[19], 2.00000))*(power(1.00000-states[19], 2.00000))+0.125000*(power(states[19], 3.00000))*(1.00000-states[19])+0.0625000*(power(states[19], 4.00000)))
    rates[22] = algebraic[28]*(1.00000-states[22])-algebraic[26]*states[22]
    rootfind_1(voi, constants, rates, states, algebraic)
    algebraic[80] = constants[6]*(algebraic[62]+algebraic[63])
    algebraic[11] = constants[7]*states[1]
    rates[1] = algebraic[80]-algebraic[11]
    algebraic[81] = constants[10]*(algebraic[64]+algebraic[63])
    rates[3] = algebraic[81]-algebraic[25]
    rates[5] = algebraic[81]-algebraic[27]
    rootfind_2(voi, constants, rates, states, algebraic)
    algebraic[48] = (constants[19]*algebraic[44]*constants[14])/(constants[23]+constants[14])
    algebraic[46] = (constants[18]*algebraic[43]*constants[14])/(constants[22]+constants[14])
    rootfind_3(voi, constants, rates, states, algebraic)
    algebraic[53] = (constants[20]*algebraic[51]*constants[14])/(constants[24]+constants[14])
    algebraic[82] = (constants[21]*constants[127]*algebraic[69])/(constants[25]+algebraic[69])
    rates[6] = (algebraic[46]+algebraic[48]+algebraic[53])-algebraic[82]
    algebraic[29] = constants[36]-states[7]
    algebraic[83] = (constants[39]*algebraic[68]*algebraic[29])/(constants[40]+algebraic[29])
    rootfind_4(voi, constants, rates, states, algebraic)
    algebraic[60] = (constants[41]*algebraic[57]*states[7])/(constants[42]+states[7])
    rates[7] = algebraic[83]-algebraic[60]
    algebraic[30] = constants[38]-states[8]
    algebraic[84] = (constants[43]*algebraic[68]*algebraic[30])/(constants[44]+algebraic[30])
    algebraic[32] = (constants[45]*states[8])/(constants[46]+states[8])
    rates[8] = algebraic[84]-algebraic[32]
    algebraic[34] = constants[49]-states[9]
    algebraic[85] = (constants[48]*constants[52]*algebraic[78]*algebraic[34])/(constants[53]+constants[48]*algebraic[34])
    algebraic[36] = (constants[48]*constants[56]*constants[51]*states[9])/(constants[57]+constants[48]*states[9])
    rates[9] = algebraic[85]-algebraic[36]
    algebraic[38] = constants[49]-states[10]
    algebraic[86] = (constants[48]*constants[52]*algebraic[78]*algebraic[38])/(constants[53]+constants[48]*algebraic[38])
    algebraic[40] = (constants[48]*constants[54]*constants[50]*states[10])/(constants[55]+constants[48]*states[10])
    rates[10] = algebraic[86]-algebraic[40]
    algebraic[88] = constants[8]*algebraic[68]*states[0]
    algebraic[23] = constants[9]*states[2]
    rates[0] = (algebraic[11]-algebraic[80])+(algebraic[23]-algebraic[88])
    rates[2] = algebraic[88]-algebraic[23]
    algebraic[35] = (1.00000/(constants[123]*constants[113]))*log(constants[63]/states[11])
    algebraic[41] = constants[66]*(power(states[14], 3.00000))*states[15]*states[16]*(states[17]-algebraic[35])
    algebraic[93] = exp(constants[83]*states[17]*constants[123])*(power(states[11], 3.00000))*constants[65]
    algebraic[94] = exp((constants[83]-1.00000)*states[17]*constants[123])*(power(constants[63], 3.00000))*states[13]
    algebraic[95] = (constants[79]/((power(constants[80], 3.00000)+power(constants[63], 3.00000))*(constants[81]+constants[65])*(1.00000+constants[82]*exp((constants[83]-1.00000)*states[17]*constants[123]))))*(algebraic[93]-algebraic[94])
    algebraic[96] = 1.00000/(1.00000+0.124500*exp(-0.100000*states[17]*constants[123])+0.0365000*constants[125]*exp(-states[17]*constants[123]))
    algebraic[97] = (((constants[84]*algebraic[96])/(1.00000+power(constants[85]/states[11], 1.50000)))*constants[64])/(constants[64]+constants[86])
    algebraic[100] = constants[90]*(states[17]-algebraic[35])
    algebraic[101] = algebraic[41]+algebraic[100]+3.00000*algebraic[95]+3.00000*algebraic[97]
    rates[11] = (-1000.00*algebraic[101]*constants[61])/(constants[58]*constants[113]*constants[112])
    algebraic[47] = (constants[76]*states[17]*constants[112]*constants[123]*(states[12]*exp(states[17]*constants[123])-constants[64]))/(exp(states[17]*constants[123])-1.00000)
    algebraic[33] = states[10]/constants[49]
    algebraic[49] = 0.500000*((0.400000*algebraic[33])/constants[121]+0.600000)
    algebraic[45] = (constants[75]*4.00000*states[17]*constants[112]*constants[123]*(0.00100000*exp(2.00000*states[17]*constants[123])-0.341000*constants[65]))/(exp(2.00000*states[17]*constants[123])-1.00000)
    algebraic[52] = algebraic[45]*constants[77]*algebraic[49]*(power(states[18], 4.00000))*states[20]*states[21]*states[22]
    algebraic[54] = (algebraic[47]/(1.00000+algebraic[52]/constants[78]))*constants[77]*algebraic[49]*(power(states[18], 4.00000))*states[20]*states[21]*states[22]
    algebraic[37] = (1.00000/(constants[123]*constants[114]))*log(constants[64]/states[12])
    algebraic[59] = constants[67]*states[23]*(0.886000*states[24]+0.114000*states[25])*(states[17]-algebraic[37])
    algebraic[61] = constants[68]*states[26]*states[27]*(states[17]-algebraic[37])
    algebraic[79] = 1.02000/(1.00000+exp(0.238500*((states[17]-algebraic[37])-59.2150)))
    algebraic[87] = (0.491240*exp(0.0803200*((states[17]+5.47600)-algebraic[37]))+exp(0.0617500*((states[17]-algebraic[37])-594.310)))/(1.00000+exp(-0.514300*((states[17]-algebraic[37])+4.75300)))
    algebraic[89] = algebraic[79]/(algebraic[79]+algebraic[87])
    algebraic[90] = constants[69]*(power(constants[64]/5.40000, 1.0/2))*algebraic[89]*(states[17]-algebraic[37])
    algebraic[91] = 1.00000/(1.00000+exp((7.48800-states[17])/5.98000))
    algebraic[92] = constants[70]*algebraic[91]*(states[17]-algebraic[37])
    algebraic[102] = ((algebraic[59]+algebraic[61]+algebraic[90]+algebraic[92])-2.00000*algebraic[97])+algebraic[54]
    rates[12] = (-1000.00*algebraic[102]*constants[61])/(constants[58]*constants[114]*constants[112])
    algebraic[98] = (constants[87]*states[13])/(constants[88]+states[13])
    algebraic[39] = (1.00000/(constants[123]*constants[115]))*log(constants[65]/states[13])
    algebraic[99] = constants[89]*(states[17]-algebraic[39])
    algebraic[103] = (algebraic[52]+algebraic[99]+algebraic[98])-2.00000*algebraic[95]
    rates[17] = -1000.00*((algebraic[103]+algebraic[102]+algebraic[101])-constants[126])
    rates[29] = custom_piecewise([greater(-1000.00*((algebraic[103]+algebraic[102]+algebraic[101])-constants[126]) , 30000.0), 1.00000-10000.0*states[29] , True, 1.00000])
    algebraic[31] = algebraic[29]/constants[36]
    algebraic[109] = (constants[94]*(1.00000+2.00000*algebraic[31]))/(1.00000+2.00000*constants[119])
    algebraic[110] = (constants[93]*(power(states[13], 2.00000)))/(power(algebraic[109], 2.00000)+power(states[13], 2.00000))
    algebraic[111] = (constants[93]*states[30])/constants[95]
    algebraic[112] = (states[30]-states[28])/constants[104]
    rates[30] = (algebraic[110]-algebraic[111])-(algebraic[112]*constants[60])/constants[59]
    algebraic[104] = states[29]+0.00200000
    algebraic[105] = 1.00000-exp(-algebraic[104]/constants[96])
    algebraic[106] = exp(-algebraic[104]/constants[97])
    algebraic[107] = constants[98]/(1.00000+exp((algebraic[103]+5.00000)/0.900000))
    algebraic[108] = algebraic[107]*algebraic[105]*algebraic[106]*(states[28]-states[13])
    algebraic[114] = 1.00000/(1.00000+(constants[102]*constants[103])/(power(constants[103]+states[28], 2.00000)))
    rates[28] = algebraic[114]*(algebraic[112]-algebraic[108])
    algebraic[113] = (constants[105]*constants[108])/(power(constants[108]+states[13], 2.00000))
    algebraic[115] = (constants[106]*constants[109])/(power(constants[109]+states[13], 2.00000))
    algebraic[117] = (constants[107]*constants[110])/(power(constants[110]+states[13], 2.00000))
    algebraic[118] = 1.00000/(1.00000+algebraic[115]+algebraic[113]+algebraic[113]+algebraic[117])
    rates[13] = -algebraic[118]*(((1000.00*algebraic[103]*constants[61])/(constants[58]*constants[115]*constants[112])+((algebraic[110]-algebraic[111])*constants[59])/constants[58])-(algebraic[108]*constants[60])/constants[58])
    return(rates)

def computeAlgebraic(constants, states, voi):
    algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
    states = array(states)
    voi = array(voi)
    algebraic[10] = 1.00000/(1.00000+exp((states[17]+87.5000)/10.3000))
    algebraic[2] = (0.320000*(states[17]+47.1300))/(1.00000-exp(-0.100000*(states[17]+47.1300)))
    algebraic[13] = 0.0800000*exp(-states[17]/11.0000)
    algebraic[3] = custom_piecewise([greater_equal(states[17] , -40.0000), 0.00000 , True, 0.135000*exp((80.0000+states[17])/-6.80000)])
    algebraic[14] = custom_piecewise([greater_equal(states[17] , -40.0000), 1.00000/(0.130000*(1.00000+exp(-(states[17]+10.6600)/11.1000))) , True, 3.56000*exp(0.0790000*states[17])+310000.*exp(0.350000*states[17])])
    algebraic[4] = custom_piecewise([greater_equal(states[17] , -40.0000), 0.00000 , True, ((-127140.*exp(0.244400*states[17])-3.47400e-05*exp(-0.0439100*states[17]))*(states[17]+37.7800))/(1.00000+exp(0.311000*(states[17]+79.2300)))])
    algebraic[15] = custom_piecewise([greater_equal(states[17] , -40.0000), (0.300000*exp(-2.57500e-07*states[17]))/(1.00000+exp(-0.100000*(states[17]+32.0000))) , True, (0.121200*exp(-0.0105200*states[17]))/(1.00000+exp(-0.137800*(states[17]+40.1400)))])
    algebraic[5] = 400.000*exp((states[17]+2.00000)/10.0000)
    algebraic[16] = 50.0000*exp((-1.00000*(states[17]+2.00000))/13.0000)
    algebraic[1] = states[9]/constants[49]
    algebraic[18] = constants[71]*((0.375000*algebraic[1])/constants[120]+0.625000)
    algebraic[6] = 1.00000/(1.00000+exp((states[17]+55.0000)/7.50000))+0.100000/(1.00000+exp((-states[17]+21.0000)/6.00000))
    algebraic[17] = 0.0200000+0.300000/(1.00000+exp((states[17]+30.0000)/9.50000))
    algebraic[7] = 1.00000/(1.00000+exp((states[17]+10.6000)/-11.4200))
    algebraic[19] = 1.00000/(45.1600*exp(0.0357700*(states[17]+50.0000))+98.9000*exp(-0.100000*(states[17]+38.0000)))
    algebraic[8] = 1.00000/(1.00000+exp((states[17]+43.5000)/6.88410))
    algebraic[20] = 0.350000*exp(-1.00000*(power((states[17]+70.0000)/15.0000, 2.00000)))+0.0350000
    algebraic[21] = 3.70000*exp(-1.00000*(power((states[17]+70.0000)/30.0000, 2.00000)))+0.0350000
    algebraic[9] = 1.00000/(1.00000+exp((states[17]+11.5000)/-11.8200))
    algebraic[22] = 10.0000/(45.1600*exp(0.0357700*(states[17]+50.0000))+98.9000*exp(-0.100000*(states[17]+38.0000)))
    algebraic[25] = constants[11]*states[3]
    algebraic[27] = constants[12]*states[4]*states[5]
    algebraic[24] = constants[73]*states[13]
    algebraic[26] = algebraic[24]*(power(1.00000-states[18], 4.00000)+2.00000*states[18]*(power(1.00000-states[18], 3.00000))+4.00000*(power(states[18], 2.00000))*(power(1.00000-states[18], 2.00000))+8.00000*(power(states[18], 3.00000))*(1.00000-states[18])+16.0000*(power(states[18], 4.00000))*(1.00000-algebraic[18]/constants[72]))
    algebraic[28] = constants[74]*(power(1.00000-states[19], 4.00000)+0.500000*states[19]*(power(1.00000-states[19], 3.00000))+0.250000*(power(states[19], 2.00000))*(power(1.00000-states[19], 2.00000))+0.125000*(power(states[19], 3.00000))*(1.00000-states[19])+0.0625000*(power(states[19], 4.00000)))
    algebraic[80] = constants[6]*(algebraic[62]+algebraic[63])
    algebraic[11] = constants[7]*states[1]
    algebraic[81] = constants[10]*(algebraic[64]+algebraic[63])
    algebraic[48] = (constants[19]*algebraic[44]*constants[14])/(constants[23]+constants[14])
    algebraic[46] = (constants[18]*algebraic[43]*constants[14])/(constants[22]+constants[14])
    algebraic[53] = (constants[20]*algebraic[51]*constants[14])/(constants[24]+constants[14])
    algebraic[82] = (constants[21]*constants[127]*algebraic[69])/(constants[25]+algebraic[69])
    algebraic[29] = constants[36]-states[7]
    algebraic[83] = (constants[39]*algebraic[68]*algebraic[29])/(constants[40]+algebraic[29])
    algebraic[60] = (constants[41]*algebraic[57]*states[7])/(constants[42]+states[7])
    algebraic[30] = constants[38]-states[8]
    algebraic[84] = (constants[43]*algebraic[68]*algebraic[30])/(constants[44]+algebraic[30])
    algebraic[32] = (constants[45]*states[8])/(constants[46]+states[8])
    algebraic[34] = constants[49]-states[9]
    algebraic[85] = (constants[48]*constants[52]*algebraic[78]*algebraic[34])/(constants[53]+constants[48]*algebraic[34])
    algebraic[36] = (constants[48]*constants[56]*constants[51]*states[9])/(constants[57]+constants[48]*states[9])
    algebraic[38] = constants[49]-states[10]
    algebraic[86] = (constants[48]*constants[52]*algebraic[78]*algebraic[38])/(constants[53]+constants[48]*algebraic[38])
    algebraic[40] = (constants[48]*constants[54]*constants[50]*states[10])/(constants[55]+constants[48]*states[10])
    algebraic[88] = constants[8]*algebraic[68]*states[0]
    algebraic[23] = constants[9]*states[2]
    algebraic[35] = (1.00000/(constants[123]*constants[113]))*log(constants[63]/states[11])
    algebraic[41] = constants[66]*(power(states[14], 3.00000))*states[15]*states[16]*(states[17]-algebraic[35])
    algebraic[93] = exp(constants[83]*states[17]*constants[123])*(power(states[11], 3.00000))*constants[65]
    algebraic[94] = exp((constants[83]-1.00000)*states[17]*constants[123])*(power(constants[63], 3.00000))*states[13]
    algebraic[95] = (constants[79]/((power(constants[80], 3.00000)+power(constants[63], 3.00000))*(constants[81]+constants[65])*(1.00000+constants[82]*exp((constants[83]-1.00000)*states[17]*constants[123]))))*(algebraic[93]-algebraic[94])
    algebraic[96] = 1.00000/(1.00000+0.124500*exp(-0.100000*states[17]*constants[123])+0.0365000*constants[125]*exp(-states[17]*constants[123]))
    algebraic[97] = (((constants[84]*algebraic[96])/(1.00000+power(constants[85]/states[11], 1.50000)))*constants[64])/(constants[64]+constants[86])
    algebraic[100] = constants[90]*(states[17]-algebraic[35])
    algebraic[101] = algebraic[41]+algebraic[100]+3.00000*algebraic[95]+3.00000*algebraic[97]
    algebraic[47] = (constants[76]*states[17]*constants[112]*constants[123]*(states[12]*exp(states[17]*constants[123])-constants[64]))/(exp(states[17]*constants[123])-1.00000)
    algebraic[33] = states[10]/constants[49]
    algebraic[49] = 0.500000*((0.400000*algebraic[33])/constants[121]+0.600000)
    algebraic[45] = (constants[75]*4.00000*states[17]*constants[112]*constants[123]*(0.00100000*exp(2.00000*states[17]*constants[123])-0.341000*constants[65]))/(exp(2.00000*states[17]*constants[123])-1.00000)
    algebraic[52] = algebraic[45]*constants[77]*algebraic[49]*(power(states[18], 4.00000))*states[20]*states[21]*states[22]
    algebraic[54] = (algebraic[47]/(1.00000+algebraic[52]/constants[78]))*constants[77]*algebraic[49]*(power(states[18], 4.00000))*states[20]*states[21]*states[22]
    algebraic[37] = (1.00000/(constants[123]*constants[114]))*log(constants[64]/states[12])
    algebraic[59] = constants[67]*states[23]*(0.886000*states[24]+0.114000*states[25])*(states[17]-algebraic[37])
    algebraic[61] = constants[68]*states[26]*states[27]*(states[17]-algebraic[37])
    algebraic[79] = 1.02000/(1.00000+exp(0.238500*((states[17]-algebraic[37])-59.2150)))
    algebraic[87] = (0.491240*exp(0.0803200*((states[17]+5.47600)-algebraic[37]))+exp(0.0617500*((states[17]-algebraic[37])-594.310)))/(1.00000+exp(-0.514300*((states[17]-algebraic[37])+4.75300)))
    algebraic[89] = algebraic[79]/(algebraic[79]+algebraic[87])
    algebraic[90] = constants[69]*(power(constants[64]/5.40000, 1.0/2))*algebraic[89]*(states[17]-algebraic[37])
    algebraic[91] = 1.00000/(1.00000+exp((7.48800-states[17])/5.98000))
    algebraic[92] = constants[70]*algebraic[91]*(states[17]-algebraic[37])
    algebraic[102] = ((algebraic[59]+algebraic[61]+algebraic[90]+algebraic[92])-2.00000*algebraic[97])+algebraic[54]
    algebraic[98] = (constants[87]*states[13])/(constants[88]+states[13])
    algebraic[39] = (1.00000/(constants[123]*constants[115]))*log(constants[65]/states[13])
    algebraic[99] = constants[89]*(states[17]-algebraic[39])
    algebraic[103] = (algebraic[52]+algebraic[99]+algebraic[98])-2.00000*algebraic[95]
    algebraic[31] = algebraic[29]/constants[36]
    algebraic[109] = (constants[94]*(1.00000+2.00000*algebraic[31]))/(1.00000+2.00000*constants[119])
    algebraic[110] = (constants[93]*(power(states[13], 2.00000)))/(power(algebraic[109], 2.00000)+power(states[13], 2.00000))
    algebraic[111] = (constants[93]*states[30])/constants[95]
    algebraic[112] = (states[30]-states[28])/constants[104]
    algebraic[104] = states[29]+0.00200000
    algebraic[105] = 1.00000-exp(-algebraic[104]/constants[96])
    algebraic[106] = exp(-algebraic[104]/constants[97])
    algebraic[107] = constants[98]/(1.00000+exp((algebraic[103]+5.00000)/0.900000))
    algebraic[108] = algebraic[107]*algebraic[105]*algebraic[106]*(states[28]-states[13])
    algebraic[114] = 1.00000/(1.00000+(constants[102]*constants[103])/(power(constants[103]+states[28], 2.00000)))
    algebraic[113] = (constants[105]*constants[108])/(power(constants[108]+states[13], 2.00000))
    algebraic[115] = (constants[106]*constants[109])/(power(constants[109]+states[13], 2.00000))
    algebraic[117] = (constants[107]*constants[110])/(power(constants[110]+states[13], 2.00000))
    algebraic[118] = 1.00000/(1.00000+algebraic[115]+algebraic[113]+algebraic[113]+algebraic[117])
    algebraic[0] = states[7]/constants[36]
    algebraic[12] = custom_piecewise([greater(voi , 59.1000) & less(voi , 59.5000), constants[117] , True, constants[116]])
    algebraic[58] = algebraic[52]+algebraic[54]
    algebraic[116] = 1000.00*(((states[28]+states[28]/algebraic[114])*constants[60])/constants[58]+(states[30]*constants[59])/constants[58])
    return algebraic

initialGuess0 = None
def rootfind_0(voi, constants, rates, states, algebraic):
    """Calculate values of algebraic variables for DAE"""
    from scipy.optimize import fsolve
    global initialGuess0
    if initialGuess0 is None: initialGuess0 = ones(3)*0.1
    if not iterable(voi):
        soln = fsolve(residualSN_0, initialGuess0, args=(algebraic, voi, constants, rates, states), xtol=1E-6)
        initialGuess0 = soln
        constants[127] = soln[0]
        constants[128] = soln[1]
        constants[129] = soln[2]
    else:
        for (i,t) in enumerate(voi):
            soln = fsolve(residualSN_0, initialGuess0, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6)
            initialGuess0 = soln
            constants[127][i] = soln[0]
            constants[128][i] = soln[1]
            constants[129][i] = soln[2]

def residualSN_0(algebraicCandidate, algebraic, voi, constants, rates, states):
    resid = array([0.0] * 3)
    constants[127] = algebraicCandidate[0]
    constants[128] = algebraicCandidate[1]
    constants[129] = algebraicCandidate[2]
    resid[0] = (constants[129]-(constants[127]*constants[128])/constants[28])
    resid[1] = (constants[127]-(constants[15]-constants[129]))
    resid[2] = (constants[128]-(constants[16]-constants[129]))
    return resid

initialGuess1 = None
def rootfind_1(voi, constants, rates, states, algebraic):
    """Calculate values of algebraic variables for DAE"""
    from scipy.optimize import fsolve
    global initialGuess1
    if initialGuess1 is None: initialGuess1 = ones(17)*0.1
    if not iterable(voi):
        soln = fsolve(residualSN_1, initialGuess1, args=(algebraic, voi, constants, rates, states), xtol=1E-6)
        initialGuess1 = soln
        algebraic[62] = soln[0]
        algebraic[63] = soln[1]
        algebraic[64] = soln[2]
        algebraic[65] = soln[3]
        algebraic[66] = soln[4]
        algebraic[67] = soln[5]
        algebraic[68] = soln[6]
        algebraic[69] = soln[7]
        algebraic[70] = soln[8]
        algebraic[71] = soln[9]
        algebraic[72] = soln[10]
        algebraic[73] = soln[11]
        algebraic[74] = soln[12]
        algebraic[75] = soln[13]
        algebraic[76] = soln[14]
        algebraic[77] = soln[15]
        algebraic[78] = soln[16]
    else:
        for (i,t) in enumerate(voi):
            soln = fsolve(residualSN_1, initialGuess1, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6)
            initialGuess1 = soln
            algebraic[62][i] = soln[0]
            algebraic[63][i] = soln[1]
            algebraic[64][i] = soln[2]
            algebraic[65][i] = soln[3]
            algebraic[66][i] = soln[4]
            algebraic[67][i] = soln[5]
            algebraic[68][i] = soln[6]
            algebraic[69][i] = soln[7]
            algebraic[70][i] = soln[8]
            algebraic[71][i] = soln[9]
            algebraic[72][i] = soln[10]
            algebraic[73][i] = soln[11]
            algebraic[74][i] = soln[12]
            algebraic[75][i] = soln[13]
            algebraic[76][i] = soln[14]
            algebraic[77][i] = soln[15]
            algebraic[78][i] = soln[16]

def residualSN_1(algebraicCandidate, algebraic, voi, constants, rates, states):
    resid = array([0.0] * 17)
    algebraic[62] = algebraicCandidate[0]
    algebraic[63] = algebraicCandidate[1]
    algebraic[64] = algebraicCandidate[2]
    algebraic[65] = algebraicCandidate[3]
    algebraic[66] = algebraicCandidate[4]
    algebraic[67] = algebraicCandidate[5]
    algebraic[68] = algebraicCandidate[6]
    algebraic[69] = algebraicCandidate[7]
    algebraic[70] = algebraicCandidate[8]
    algebraic[71] = algebraicCandidate[9]
    algebraic[72] = algebraicCandidate[10]
    algebraic[73] = algebraicCandidate[11]
    algebraic[74] = algebraicCandidate[12]
    algebraic[75] = algebraicCandidate[13]
    algebraic[76] = algebraicCandidate[14]
    algebraic[77] = algebraicCandidate[15]
    algebraic[78] = algebraicCandidate[16]
    resid[0] = (algebraic[62]-(algebraic[65]*algebraic[66])/constants[3])
    resid[1] = (algebraic[63]-(algebraic[62]*algebraic[67])/constants[4])
    resid[2] = (algebraic[64]-(algebraic[66]*algebraic[67])/constants[5])
    resid[3] = (algebraic[65]-((constants[0]-algebraic[62])-algebraic[63]))
    resid[4] = (algebraic[66]-(((states[0]-algebraic[62])-algebraic[63])-algebraic[64]))
    resid[5] = (algebraic[67]-((constants[2]-algebraic[63])-algebraic[64]))
    resid[6] = (algebraic[70]-(constants[31]*constants[35])/(constants[35]+algebraic[68]+algebraic[78]))
    resid[7] = (algebraic[71]-(algebraic[68]/constants[34])*algebraic[68]*(1.00000+algebraic[70]/constants[35]))
    resid[8] = (algebraic[72]-algebraic[68]*(1.00000+algebraic[70]/constants[35]))
    resid[9] = (algebraic[73]-(algebraic[78]/constants[34])*algebraic[78]*(1.00000+algebraic[70]/constants[35]))
    resid[10] = (algebraic[74]-algebraic[78]*(1.00000+algebraic[70]/constants[35]))
    resid[11] = (algebraic[75]-(constants[32]/algebraic[69])*algebraic[71])
    resid[12] = (algebraic[76]-(constants[32]/algebraic[69])*algebraic[73])
    resid[13] = (algebraic[77]-((constants[32]*constants[33])/constants[34]+(constants[32]*algebraic[69])/constants[34]+(power(algebraic[69], 2.00000))/constants[34]))
    resid[14] = (algebraic[69]-((states[6]-(algebraic[75]+2.00000*algebraic[71]+2.00000*algebraic[72]))-(algebraic[76]+2.00000*algebraic[73]+2.00000*algebraic[74])))
    resid[15] = (0.00000-(2.00000*constants[29]*(power(algebraic[69], 2.00000))-algebraic[68]*(1.00000+algebraic[70]/constants[35])*(algebraic[77]*algebraic[68]+power(algebraic[69], 2.00000))))
    resid[16] = (0.00000-(2.00000*constants[30]*(power(algebraic[69], 2.00000))-algebraic[78]*(1.00000+algebraic[70]/constants[35])*(algebraic[77]*algebraic[78]+power(algebraic[69], 2.00000))))
    return resid

initialGuess2 = None
def rootfind_2(voi, constants, rates, states, algebraic):
    """Calculate values of algebraic variables for DAE"""
    from scipy.optimize import fsolve
    global initialGuess2
    if initialGuess2 is None: initialGuess2 = ones(3)*0.1
    if not iterable(voi):
        soln = fsolve(residualSN_2, initialGuess2, args=(algebraic, voi, constants, rates, states), xtol=1E-6)
        initialGuess2 = soln
        algebraic[42] = soln[0]
        algebraic[43] = soln[1]
        algebraic[44] = soln[2]
    else:
        for (i,t) in enumerate(voi):
            soln = fsolve(residualSN_2, initialGuess2, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6)
            initialGuess2 = soln
            algebraic[42][i] = soln[0]
            algebraic[43][i] = soln[1]
            algebraic[44][i] = soln[2]

def residualSN_2(algebraicCandidate, algebraic, voi, constants, rates, states):
    resid = array([0.0] * 3)
    algebraic[42] = algebraicCandidate[0]
    algebraic[43] = algebraicCandidate[1]
    algebraic[44] = algebraicCandidate[2]
    resid[0] = (algebraic[44]-(algebraic[42]*algebraic[43])/constants[26])
    resid[1] = (algebraic[42]-(states[3]-algebraic[44]))
    resid[2] = (algebraic[43]-(constants[13]-algebraic[44]))
    return resid

initialGuess3 = None
def rootfind_3(voi, constants, rates, states, algebraic):
    """Calculate values of algebraic variables for DAE"""
    from scipy.optimize import fsolve
    global initialGuess3
    if initialGuess3 is None: initialGuess3 = ones(2)*0.1
    if not iterable(voi):
        soln = fsolve(residualSN_3, initialGuess3, args=(algebraic, voi, constants, rates, states), xtol=1E-6)
        initialGuess3 = soln
        algebraic[50] = soln[0]
        algebraic[51] = soln[1]
    else:
        for (i,t) in enumerate(voi):
            soln = fsolve(residualSN_3, initialGuess3, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6)
            initialGuess3 = soln
            algebraic[50][i] = soln[0]
            algebraic[51][i] = soln[1]

def residualSN_3(algebraicCandidate, algebraic, voi, constants, rates, states):
    resid = array([0.0] * 2)
    algebraic[50] = algebraicCandidate[0]
    algebraic[51] = algebraicCandidate[1]
    resid[0] = (algebraic[51]-(algebraic[50]*algebraic[43])/constants[27])
    resid[1] = (algebraic[50]-(constants[17]-algebraic[51]))
    return resid

initialGuess4 = None
def rootfind_4(voi, constants, rates, states, algebraic):
    """Calculate values of algebraic variables for DAE"""
    from scipy.optimize import fsolve
    global initialGuess4
    if initialGuess4 is None: initialGuess4 = ones(3)*0.1
    if not iterable(voi):
        soln = fsolve(residualSN_4, initialGuess4, args=(algebraic, voi, constants, rates, states), xtol=1E-6)
        initialGuess4 = soln
        algebraic[55] = soln[0]
        algebraic[56] = soln[1]
        algebraic[57] = soln[2]
    else:
        for (i,t) in enumerate(voi):
            soln = fsolve(residualSN_4, initialGuess4, args=(algebraic[:,i], voi[i], constants, rates[:i], states[:,i]), xtol=1E-6)
            initialGuess4 = soln
            algebraic[55][i] = soln[0]
            algebraic[56][i] = soln[1]
            algebraic[57][i] = soln[2]

def residualSN_4(algebraicCandidate, algebraic, voi, constants, rates, states):
    resid = array([0.0] * 3)
    algebraic[55] = algebraicCandidate[0]
    algebraic[56] = algebraicCandidate[1]
    algebraic[57] = algebraicCandidate[2]
    resid[0] = (algebraic[55]-(algebraic[56]*algebraic[57])/constants[47])
    resid[1] = (algebraic[56]-(states[8]-algebraic[55]))
    resid[2] = (algebraic[57]-(constants[37]-algebraic[55]))
    return resid

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)