# Size of variable arrays:
sizeAlgebraic = 53
sizeStates = 7
sizeConstants = 114
from math import *
from numpy import *

def createLegends():
    legend_states = [""] * sizeStates
    legend_rates = [""] * sizeStates
    legend_algebraic = [""] * sizeAlgebraic
    legend_voi = ""
    legend_constants = [""] * sizeConstants
    legend_voi = "time in component environment (day)"
    legend_states[0] = "OBp in component OBp (pM)"
    legend_constants[0] = "OBp_t0 in component OBp (pM)"
    legend_algebraic[49] = "OBp_in in component OBp (flux)"
    legend_algebraic[9] = "OBp_out in component OBp (flux)"
    legend_constants[1] = "OBu_t0 in component model_parameters (pM)"
    legend_constants[77] = "DifferOBpNormal in component OBp (first_order_rate_constant)"
    legend_constants[2] = "D_OBp_t0 in component model_parameters (first_order_rate_constant)"
    legend_constants[3] = "TGFbNormal in component OBp (pM)"
    legend_constants[4] = "KD_TGF_beta_repress in component model_parameters (pM)"
    legend_constants[81] = "ProlifOBpNormal in component OBp (first_order_rate_constant)"
    legend_constants[5] = "frac_prolifOBp_vs_differOBu in component OBp (dimensionless)"
    legend_constants[6] = "OBpNormal in component OBp (pM)"
    legend_constants[7] = "OBp_sat in component OBp (pM)"
    legend_constants[86] = "DifferOBu in component OBp (first_order_rate_constant)"
    legend_algebraic[1] = "ProlifOBp in component OBp (first_order_rate_constant)"
    legend_constants[8] = "D_OBu_t0 in component model_parameters (first_order_rate_constant)"
    legend_constants[9] = "pd_OBp_t0 in component model_parameters (flux)"
    legend_constants[82] = "Pi_TGFbeta_OBu_act_t0 in component TGF_beta (dimensionless)"
    legend_algebraic[8] = "Pi_TGFbeta_OBp_rep in component TGF_beta (dimensionless)"
    legend_algebraic[46] = "Pi_WNT in component Wnt (dimensionless)"
    legend_constants[112] = "Pi_WNT_0 in component Wnt (dimensionless)"
    legend_states[1] = "OBa in component OBa (pM)"
    legend_algebraic[10] = "OBa_in in component OBa (flux)"
    legend_algebraic[2] = "OBa_out in component OBa (flux)"
    legend_constants[10] = "pd_OBa_t0 in component model_parameters (flux)"
    legend_constants[11] = "D_OBa_t0 in component model_parameters (first_order_rate_constant)"
    legend_constants[12] = "A_OBa_t0 in component model_parameters (first_order_rate_constant)"
    legend_states[2] = "OCp in component OCp (pM)"
    legend_algebraic[52] = "OCp_in in component OCp (flux)"
    legend_algebraic[47] = "OCp_out in component OCp (flux)"
    legend_constants[13] = "OCu_t0 in component model_parameters (pM)"
    legend_algebraic[44] = "Pi_RANKL_act_OCp in component RANK_RANKL_OPG (dimensionless)"
    legend_algebraic[50] = "Pi_RANKL_act_OCu in component RANK_RANKL_OPG (dimensionless)"
    legend_constants[14] = "D_OCu_t0 in component model_parameters (first_order_rate_constant)"
    legend_constants[15] = "pd_OCp_t0 in component model_parameters (flux)"
    legend_constants[16] = "D_OCp_t0 in component model_parameters (first_order_rate_constant)"
    legend_states[3] = "OCa in component OCa (pM)"
    legend_algebraic[51] = "OCa_in in component OCa (flux)"
    legend_algebraic[7] = "OCa_out in component OCa (flux)"
    legend_constants[17] = "pd_OCa_t0 in component model_parameters (flux)"
    legend_constants[18] = "A_OCa_t0 in component model_parameters (first_order_rate_constant)"
    legend_algebraic[5] = "Pi_TGFbeta_OCa_act in component TGF_beta (dimensionless)"
    legend_states[4] = "fbm in component fbm (dimensionless)"
    legend_constants[19] = "K_form in component model_parameters (second_order_rate_constant)"
    legend_constants[20] = "K_res in component model_parameters (second_order_rate_constant)"
    legend_algebraic[0] = "dfbmdt in component fbm (first_order_rate_constant)"
    legend_states[5] = "OCY in component OCY (pM)"
    legend_algebraic[3] = "OCY_act in component OCY (pM)"
    legend_constants[79] = "OCY_act_0 in component OCY (pM)"
    legend_constants[21] = "fact_0 in component model_parameters (pM)"
    legend_states[6] = "vm in component vm (dimensionless)"
    legend_constants[22] = "XKAPPA in component model_parameters (first_order_rate_constant)"
    legend_constants[23] = "vmmax in component model_parameters (dimensionless)"
    legend_constants[80] = "PTH_tot in component PTH (pM)"
    legend_constants[83] = "Pi_PTH_act in component PTH (dimensionless)"
    legend_constants[84] = "Pi_PTH_rep in component PTH (dimensionless)"
    legend_constants[24] = "Beta_PTH in component model_parameters (flux)"
    legend_constants[25] = "P_PTH_d in component model_parameters (flux)"
    legend_constants[26] = "Deg_PTH in component model_parameters (first_order_rate_constant)"
    legend_constants[27] = "KD_PTH_act in component model_parameters (pM)"
    legend_constants[28] = "KD_PTH_rep in component model_parameters (pM)"
    legend_algebraic[36] = "RANK in component RANK_RANKL_OPG (pM)"
    legend_algebraic[37] = "RANKL in component RANK_RANKL_OPG (pM)"
    legend_algebraic[38] = "OPG in component RANK_RANKL_OPG (pM)"
    legend_algebraic[39] = "RANKL_tot in component RANK_RANKL_OPG (pM)"
    legend_constants[85] = "RANKL_max in component RANK_RANKL_OPG (pM)"
    legend_algebraic[40] = "P_RANKL1 in component RANK_RANKL_OPG (flux)"
    legend_algebraic[41] = "P_RANKL2 in component RANK_RANKL_OPG (flux)"
    legend_algebraic[42] = "P_RANKL in component RANK_RANKL_OPG (flux)"
    legend_constants[29] = "P_RANKL_d in component model_parameters (flux)"
    legend_algebraic[34] = "Pi_NO_PTH_act_rep in component NO_PTH (dimensionless)"
    legend_constants[30] = "N_RANK_OCp in component model_parameters (dimensionless)"
    legend_constants[31] = "N_RANKL_OBp_max in component model_parameters (dimensionless)"
    legend_constants[32] = "K_RANK_RANKL in component model_parameters (pM)"
    legend_constants[33] = "Beta_OPG in component model_parameters (first_order_rate_constant)"
    legend_constants[34] = "OPG_max in component model_parameters (pM)"
    legend_constants[35] = "Deg_OPG in component model_parameters (first_order_rate_constant)"
    legend_constants[36] = "Deg_OPG_RANKL in component model_parameters (first_order_rate_constant)"
    legend_constants[37] = "K_OPG_RANKL in component model_parameters (pM)"
    legend_constants[38] = "Beta_RANKL_OCY in component model_parameters (first_order_rate_constant)"
    legend_constants[39] = "Beta_RANKL_OBp in component model_parameters (first_order_rate_constant)"
    legend_constants[40] = "Deg_RANKL in component model_parameters (first_order_rate_constant)"
    legend_constants[41] = "Deg_RANK_RANKL in component model_parameters (first_order_rate_constant)"
    legend_constants[42] = "KD_RANKL_act_OCp in component model_parameters (pM)"
    legend_constants[43] = "KD_RANKL_act_OCu in component model_parameters (pM)"
    legend_algebraic[4] = "TGF_beta in component TGF_beta (pM)"
    legend_constants[78] = "TGF_beta_t0 in component TGF_beta (pM)"
    legend_constants[44] = "OCa_t0 in component TGF_beta (pM)"
    legend_constants[45] = "Alpha in component model_parameters (dimensionless)"
    legend_constants[46] = "KD_TGF_beta_activate in component model_parameters (pM)"
    legend_algebraic[6] = "Pi_TGFbeta_OBu_act in component TGF_beta (dimensionless)"
    legend_algebraic[45] = "Scl in component Scl (pM)"
    legend_constants[111] = "Scl_0 in component Scl (pM)"
    legend_algebraic[48] = "P_Scl_b in component Scl (flux)"
    legend_constants[113] = "P_Scl_b_0 in component Scl (flux)"
    legend_algebraic[33] = "A in component Scl (second_order_rate_constant)"
    legend_constants[108] = "A_0 in component Scl (second_order_rate_constant)"
    legend_algebraic[35] = "B in component Scl (first_order_rate_constant)"
    legend_constants[109] = "B_0 in component Scl (first_order_rate_constant)"
    legend_algebraic[43] = "C in component Scl (flux)"
    legend_constants[110] = "C_0 in component Scl (flux)"
    legend_algebraic[31] = "OCYprod in component Scl (flux)"
    legend_constants[107] = "OCYprod_0 in component Scl (flux)"
    legend_constants[89] = "Wnt in component Wnt (pM)"
    legend_constants[47] = "Wnt_0 in component model_parameters (pM)"
    legend_algebraic[29] = "Pi_eps_rep in component SED (dimensionless)"
    legend_constants[106] = "Pi_eps_rep_stst in component SED (dimensionless)"
    legend_constants[87] = "Beta_Scl in component Scl (first_order_rate_constant)"
    legend_constants[48] = "Beta_Scl_0 in component model_parameters (first_order_rate_constant)"
    legend_constants[49] = "KD_SclLRP5 in component model_parameters (pM)"
    legend_constants[50] = "Deg_Scl in component model_parameters (first_order_rate_constant)"
    legend_constants[51] = "Scl_max in component model_parameters (pM)"
    legend_constants[52] = "KD_WntLRP5 in component model_parameters (pM)"
    legend_constants[53] = "Deg_SclLRP5 in component model_parameters (first_order_rate_constant)"
    legend_constants[54] = "LRP5perCell in component model_parameters (dimensionless)"
    legend_algebraic[11] = "LRP5tot in component Scl (pM)"
    legend_constants[88] = "LRP5tot_0 in component Scl (pM)"
    legend_constants[55] = "P_Scl_d in component model_parameters (flux)"
    legend_constants[90] = "PTH_tot in component NO_PTH (pM)"
    legend_algebraic[30] = "NO_tot in component NO_PTH (pM)"
    legend_constants[103] = "NO_eq in component NO_PTH (pM)"
    legend_constants[91] = "Beta_NO in component NO_PTH (first_order_rate_constant)"
    legend_constants[56] = "Beta_NO_0 in component model_parameters (first_order_rate_constant)"
    legend_constants[57] = "P_NO_d in component model_parameters (flux)"
    legend_constants[58] = "Deg_NO in component model_parameters (first_order_rate_constant)"
    legend_algebraic[28] = "Pi_eps_act in component SED (dimensionless)"
    legend_constants[102] = "Pi_eps_act_stst in component SED (dimensionless)"
    legend_constants[59] = "NO_max in component model_parameters (pM)"
    legend_constants[104] = "KD_rep_NO in component NO_PTH (pM)"
    legend_constants[60] = "aa in component NO_PTH (dimensionless)"
    legend_constants[92] = "Pi_PTH_act in component NO_PTH (dimensionless)"
    legend_algebraic[32] = "Pi_NO_rep in component NO_PTH (dimensionless)"
    legend_constants[61] = "lambda_s in component model_parameters (dimensionless)"
    legend_constants[62] = "lambda_c in component model_parameters (dimensionless)"
    legend_algebraic[26] = "SED in component SED (MPa)"
    legend_algebraic[27] = "SED_bm in component SED (MPa)"
    legend_constants[63] = "sig_macro_t0 in component model_parameters (MPa)"
    legend_constants[64] = "de_sig_macro in component model_parameters (MPa)"
    legend_constants[98] = "sig_macro in component SED (MPa)"
    legend_constants[93] = "sig_1 in component SED (MPa)"
    legend_constants[94] = "sig_2 in component SED (MPa)"
    legend_constants[99] = "sig_3 in component SED (MPa)"
    legend_constants[95] = "sig_4 in component SED (MPa)"
    legend_constants[96] = "sig_5 in component SED (MPa)"
    legend_constants[97] = "sig_6 in component SED (MPa)"
    legend_algebraic[20] = "eps_1 in component SED (dimensionless)"
    legend_algebraic[21] = "eps_2 in component SED (dimensionless)"
    legend_algebraic[22] = "eps_3 in component SED (dimensionless)"
    legend_algebraic[23] = "eps_4 in component SED (dimensionless)"
    legend_algebraic[24] = "eps_5 in component SED (dimensionless)"
    legend_algebraic[25] = "eps_6 in component SED (dimensionless)"
    legend_algebraic[12] = "dens_tis in component SED (g_per_cm3)"
    legend_constants[65] = "vo in component model_parameters (dimensionless)"
    legend_constants[66] = "rho_m in component model_parameters (g_per_cm3)"
    legend_constants[67] = "rho_o in component model_parameters (g_per_cm3)"
    legend_algebraic[13] = "fvas in component SED (dimensionless)"
    legend_algebraic[14] = "dens in component SED (g_per_cm3)"
    legend_algebraic[15] = "dens_1 in component SED (dimensionless)"
    legend_algebraic[16] = "E_mod in component SED (MPa)"
    legend_constants[68] = "nu in component model_parameters (dimensionless)"
    legend_algebraic[17] = "S_11 in component SED (per_MPa)"
    legend_algebraic[18] = "S_12 in component SED (per_MPa)"
    legend_algebraic[19] = "S_44 in component SED (per_MPa)"
    legend_constants[69] = "omega in component model_parameters (dimensionless)"
    legend_constants[100] = "valeq in component SED (dimensionless)"
    legend_constants[70] = "taueq in component model_parameters (MPa)"
    legend_constants[71] = "alphaAct in component model_parameters (dimensionless)"
    legend_constants[72] = "alphaRep in component model_parameters (dimensionless)"
    legend_constants[73] = "gammaAct in component model_parameters (dimensionless)"
    legend_constants[74] = "gammaRep in component model_parameters (dimensionless)"
    legend_constants[75] = "rhoAct in component model_parameters (dimensionless)"
    legend_constants[76] = "rhoRep in component model_parameters (dimensionless)"
    legend_constants[101] = "deltaAct in component SED (MPa)"
    legend_constants[105] = "deltaRep in component SED (MPa)"
    legend_rates[0] = "d/dt OBp in component OBp (pM)"
    legend_rates[1] = "d/dt OBa in component OBa (pM)"
    legend_rates[2] = "d/dt OCp in component OCp (pM)"
    legend_rates[3] = "d/dt OCa in component OCa (pM)"
    legend_rates[4] = "d/dt fbm in component fbm (dimensionless)"
    legend_rates[5] = "d/dt OCY in component OCY (pM)"
    legend_rates[6] = "d/dt vm in component vm (dimensionless)"
    return (legend_states, legend_algebraic, legend_voi, legend_constants)

def initConsts():
    constants = [0.0] * sizeConstants; states = [0.0] * sizeStates;
    states[0] = 1.1631869976e-03
    constants[0] = 1.1631869976e-03
    constants[1] = 0.01
    constants[2] = 0.185
    constants[3] = 0.0001
    constants[4] = 0.000175426051821094
    constants[5] = 0.5
    constants[6] = 0.001
    constants[7] = 0.005
    constants[8] = 0.166
    constants[9] = 0
    states[1] = 9.1880026722e-04
    constants[10] = 0
    constants[11] = 0.0212
    constants[12] = 0.1908
    states[2] = 1.3556879825e-03
    constants[13] = 0.01
    constants[14] = 0.0219
    constants[15] = 0
    constants[16] = 0.01958
    states[3] = 1.8376005344e-05
    constants[17] = 0
    constants[18] = 10
    states[4] = 15
    constants[19] = 50
    constants[20] = 2500
    states[5] = 1.2734398187e-02
    constants[21] = 0.000414138
    states[6] = 3.5893417161e-01
    constants[22] = 0.007
    constants[23] = 0.516
    constants[24] = 250
    constants[25] = 0
    constants[26] = 86
    constants[27] = 0.65
    constants[28] = 0.222581427709954
    constants[29] = 0
    constants[30] = 4160
    constants[31] = 2703476
    constants[32] = 10
    constants[33] = 162490033.783568
    constants[34] = 131.428571428571
    constants[35] = 532608.695652174
    constants[36] = 10.132471014805
    constants[37] = 0.0151142857142857
    constants[38] = 5660
    constants[39] = 23600
    constants[40] = 10.132471014805
    constants[41] = 10.132471014805
    constants[42] = 3.34
    constants[43] = 16.7
    constants[44] = 0.0001
    constants[45] = 1
    constants[46] = 0.000563278809675429
    constants[47] = 170
    constants[48] = 24000
    constants[49] = 10
    constants[50] = 1
    constants[51] = 70
    constants[52] = 1000
    constants[53] = 50
    constants[54] = 5
    constants[55] = 0
    constants[56] = 3440
    constants[57] = 0
    constants[58] = 0.0021
    constants[59] = 200000000
    constants[60] = 2
    constants[61] = 0.450454
    constants[62] = 0.900909
    constants[63] = -0.1457
    constants[64] = 0
    constants[65] = 0.428571428571428
    constants[66] = 3.2
    constants[67] = 1.41
    constants[68] = 0.3
    constants[69] = 0.95
    constants[70] = 0.006652
    constants[71] = 1
    constants[72] = 1
    constants[73] = 7
    constants[74] = 8.01559
    constants[75] = 0
    constants[76] = 0
    constants[77] = (constants[2]*1.00000)/(1.00000+constants[3]/constants[4])
    constants[78] = constants[45]*constants[44]
    constants[79] = constants[21]*20.0000
    constants[80] = (constants[24]+constants[25])/constants[26]
    constants[81] = (constants[5]*constants[77])/(1.00000-constants[6]/constants[7])
    constants[82] = constants[78]/(constants[46]+constants[78])
    constants[83] = constants[80]/(constants[80]+constants[27])
    constants[84] = 1.00000/(1.00000+constants[80]/constants[28])
    constants[85] = constants[31]*constants[0]
    constants[86] = (1.00000-constants[5])*constants[8]*constants[82]
    constants[87] = constants[48]
    constants[88] = constants[54]*constants[0]
    constants[89] = constants[47]
    constants[90] = (constants[24]+constants[25])/constants[26]
    constants[91] = constants[56]
    constants[92] = constants[90]/(constants[90]+constants[27])
    constants[93] = 0.00000
    constants[94] = 0.00000
    constants[95] = 0.00000
    constants[96] = 0.00000
    constants[97] = 0.00000
    constants[98] = constants[63]+constants[64]
    constants[99] = constants[98]
    constants[100] = constants[69]
    constants[101] = constants[70]*(power((constants[71]-constants[100])/(constants[100]-constants[75]), 1.00000/constants[73]))
    constants[102] = constants[75]+((constants[71]-constants[75])*(power(constants[70], constants[73])))/(power(constants[101], constants[73])+power(constants[70], constants[73]))
    constants[103] = (constants[56]*constants[102]*constants[79])/(constants[58]+(constants[56]*constants[102]*constants[79])/constants[59])
    constants[104] = constants[103]/constants[60]
    constants[105] = constants[70]*(power((constants[100]-constants[76])/(constants[72]-constants[100]), 1.00000/constants[74]))
    constants[106] = constants[72]-((constants[72]-constants[76])*(power(constants[70], constants[74])))/(power(constants[105], constants[74])+power(constants[70], constants[74]))
    constants[107] = constants[48]*constants[79]*constants[106]
    constants[108] = (1.00000/constants[49])*(constants[50]+constants[107]/constants[51])
    constants[109] = (constants[108]*constants[49]*(1.00000+constants[47]/constants[52])+(constants[53]*constants[88])/constants[49])-(constants[55]+constants[107])/constants[49]
    constants[110] = -(constants[55]+constants[107])*(1.00000+constants[47]/constants[52])
    constants[111] = (-constants[109]+power(power(constants[109], 2.00000)-4.00000*constants[108]*constants[110], 1.0/2))/(2.00000*constants[108])
    constants[112] = constants[47]/(constants[52]*(1.00000+constants[47]/constants[52]+constants[111]/constants[49]))
    constants[113] = constants[107]*(1.00000-constants[111]/constants[51])
    return (states, constants)

def computeRates(voi, states, constants):
    rates = [0.0] * sizeStates; algebraic = [0.0] * sizeAlgebraic
    rates[6] = constants[22]*(constants[23]-states[6])-(states[6]*constants[19]*states[1])/states[4]
    algebraic[0] = constants[19]*states[1]-constants[20]*states[3]
    rates[4] = algebraic[0]
    rates[5] = constants[21]*algebraic[0]
    algebraic[4] = constants[45]*states[3]
    algebraic[8] = 1.00000/(1.00000+algebraic[4]/constants[4])
    algebraic[9] = constants[2]*algebraic[8]*states[0]
    algebraic[10] = algebraic[9]+constants[10]
    algebraic[2] = (constants[11]+constants[12])*states[1]
    rates[1] = algebraic[10]-algebraic[2]
    algebraic[1] = constants[81]*(1.00000-states[0]/constants[7])
    algebraic[3] = constants[21]*states[4]
    algebraic[12] = 1.00000+(constants[66]-1.00000)*states[6]+(constants[67]-1.00000)*constants[65]
    algebraic[14] = (algebraic[12]*states[4])/100.000
    algebraic[15] = algebraic[14]/1.00000
    algebraic[16] = custom_piecewise([less(algebraic[15] , 1.20000), 2014.00*(power(algebraic[15], 2.50000)) , True, 1763.00*(power(algebraic[15], 3.20000))])
    algebraic[17] = 1.00000/algebraic[16]
    algebraic[18] = -constants[68]/algebraic[16]
    algebraic[20] = algebraic[17]*constants[93]+algebraic[18]*constants[94]+algebraic[18]*constants[99]
    algebraic[21] = algebraic[18]*constants[93]+algebraic[17]*constants[94]+algebraic[18]*constants[99]
    algebraic[22] = algebraic[18]*constants[93]+algebraic[18]*constants[94]+algebraic[17]*constants[99]
    algebraic[19] = (2.00000+2.00000*constants[68])/algebraic[16]
    algebraic[23] = algebraic[19]*constants[95]
    algebraic[24] = algebraic[19]*constants[96]
    algebraic[25] = algebraic[19]*constants[97]
    algebraic[26] = 0.500000*(constants[93]*algebraic[20]+constants[94]*algebraic[21]+constants[99]*algebraic[22]+constants[95]*algebraic[23]+constants[96]*algebraic[24]+constants[97]*algebraic[25])
    algebraic[13] = 1.00000-states[4]/100.000
    algebraic[27] = algebraic[26]/(power(1.00000-algebraic[13], 2.00000))
    algebraic[29] = constants[72]-((constants[72]-constants[76])*(power(algebraic[27], constants[74])))/(power(constants[105], constants[74])+power(algebraic[27], constants[74]))
    algebraic[31] = constants[87]*algebraic[3]*algebraic[29]
    algebraic[33] = (1.00000/constants[49])*(constants[50]+algebraic[31]/constants[51])
    algebraic[11] = constants[54]*states[0]
    algebraic[35] = (algebraic[33]*constants[49]*(1.00000+constants[89]/constants[52])+(constants[53]*algebraic[11])/constants[49])-(constants[55]+algebraic[31])/constants[49]
    algebraic[43] = -(constants[55]+algebraic[31])*(1.00000+constants[89]/constants[52])
    algebraic[45] = (-algebraic[35]+power(power(algebraic[35], 2.00000)-4.00000*algebraic[33]*algebraic[43], 1.0/2))/(2.00000*algebraic[33])
    algebraic[46] = constants[89]/(constants[52]*(1.00000+constants[89]/constants[52]+algebraic[45]/constants[49]))
    algebraic[49] = constants[86]*constants[1]+(algebraic[1]*states[0]*algebraic[46])/constants[112]+constants[9]
    rates[0] = algebraic[49]-algebraic[9]
    algebraic[28] = constants[75]+((constants[71]-constants[75])*(power(algebraic[27], constants[73])))/(power(constants[101], constants[73])+power(algebraic[27], constants[73]))
    algebraic[30] = (constants[91]*algebraic[28]*algebraic[3]+constants[57])/(constants[58]+(constants[91]*algebraic[28]*algebraic[3])/constants[59])
    algebraic[32] = 1.00000/(1.00000+algebraic[30]/constants[104])
    algebraic[34] = constants[61]*(constants[92]+algebraic[32])+constants[62]*constants[92]*algebraic[32]
    rootfind_0(voi, constants, rates, states, algebraic)
    algebraic[44] = algebraic[37]/(constants[42]+algebraic[37])
    algebraic[47] = constants[16]*algebraic[44]*states[2]
    algebraic[51] = algebraic[47]+constants[17]
    algebraic[5] = algebraic[4]/(constants[46]+algebraic[4])
    algebraic[7] = constants[18]*algebraic[5]*states[3]
    rates[3] = algebraic[51]-algebraic[7]
    algebraic[50] = algebraic[37]/(constants[43]+algebraic[37])
    algebraic[52] = constants[14]*constants[13]*algebraic[50]+constants[15]
    rates[2] = algebraic[52]-algebraic[47]
    return(rates)

def computeAlgebraic(constants, states, voi):
    algebraic = array([[0.0] * len(voi)] * sizeAlgebraic)
    states = array(states)
    voi = array(voi)
    algebraic[0] = constants[19]*states[1]-constants[20]*states[3]
    algebraic[4] = constants[45]*states[3]
    algebraic[8] = 1.00000/(1.00000+algebraic[4]/constants[4])
    algebraic[9] = constants[2]*algebraic[8]*states[0]
    algebraic[10] = algebraic[9]+constants[10]
    algebraic[2] = (constants[11]+constants[12])*states[1]
    algebraic[1] = constants[81]*(1.00000-states[0]/constants[7])
    algebraic[3] = constants[21]*states[4]
    algebraic[12] = 1.00000+(constants[66]-1.00000)*states[6]+(constants[67]-1.00000)*constants[65]
    algebraic[14] = (algebraic[12]*states[4])/100.000
    algebraic[15] = algebraic[14]/1.00000
    algebraic[16] = custom_piecewise([less(algebraic[15] , 1.20000), 2014.00*(power(algebraic[15], 2.50000)) , True, 1763.00*(power(algebraic[15], 3.20000))])
    algebraic[17] = 1.00000/algebraic[16]
    algebraic[18] = -constants[68]/algebraic[16]
    algebraic[20] = algebraic[17]*constants[93]+algebraic[18]*constants[94]+algebraic[18]*constants[99]
    algebraic[21] = algebraic[18]*constants[93]+algebraic[17]*constants[94]+algebraic[18]*constants[99]
    algebraic[22] = algebraic[18]*constants[93]+algebraic[18]*constants[94]+algebraic[17]*constants[99]
    algebraic[19] = (2.00000+2.00000*constants[68])/algebraic[16]
    algebraic[23] = algebraic[19]*constants[95]
    algebraic[24] = algebraic[19]*constants[96]
    algebraic[25] = algebraic[19]*constants[97]
    algebraic[26] = 0.500000*(constants[93]*algebraic[20]+constants[94]*algebraic[21]+constants[99]*algebraic[22]+constants[95]*algebraic[23]+constants[96]*algebraic[24]+constants[97]*algebraic[25])
    algebraic[13] = 1.00000-states[4]/100.000
    algebraic[27] = algebraic[26]/(power(1.00000-algebraic[13], 2.00000))
    algebraic[29] = constants[72]-((constants[72]-constants[76])*(power(algebraic[27], constants[74])))/(power(constants[105], constants[74])+power(algebraic[27], constants[74]))
    algebraic[31] = constants[87]*algebraic[3]*algebraic[29]
    algebraic[33] = (1.00000/constants[49])*(constants[50]+algebraic[31]/constants[51])
    algebraic[11] = constants[54]*states[0]
    algebraic[35] = (algebraic[33]*constants[49]*(1.00000+constants[89]/constants[52])+(constants[53]*algebraic[11])/constants[49])-(constants[55]+algebraic[31])/constants[49]
    algebraic[43] = -(constants[55]+algebraic[31])*(1.00000+constants[89]/constants[52])
    algebraic[45] = (-algebraic[35]+power(power(algebraic[35], 2.00000)-4.00000*algebraic[33]*algebraic[43], 1.0/2))/(2.00000*algebraic[33])
    algebraic[46] = constants[89]/(constants[52]*(1.00000+constants[89]/constants[52]+algebraic[45]/constants[49]))
    algebraic[49] = constants[86]*constants[1]+(algebraic[1]*states[0]*algebraic[46])/constants[112]+constants[9]
    algebraic[28] = constants[75]+((constants[71]-constants[75])*(power(algebraic[27], constants[73])))/(power(constants[101], constants[73])+power(algebraic[27], constants[73]))
    algebraic[30] = (constants[91]*algebraic[28]*algebraic[3]+constants[57])/(constants[58]+(constants[91]*algebraic[28]*algebraic[3])/constants[59])
    algebraic[32] = 1.00000/(1.00000+algebraic[30]/constants[104])
    algebraic[34] = constants[61]*(constants[92]+algebraic[32])+constants[62]*constants[92]*algebraic[32]
    algebraic[44] = algebraic[37]/(constants[42]+algebraic[37])
    algebraic[47] = constants[16]*algebraic[44]*states[2]
    algebraic[51] = algebraic[47]+constants[17]
    algebraic[5] = algebraic[4]/(constants[46]+algebraic[4])
    algebraic[7] = constants[18]*algebraic[5]*states[3]
    algebraic[50] = algebraic[37]/(constants[43]+algebraic[37])
    algebraic[52] = constants[14]*constants[13]*algebraic[50]+constants[15]
    algebraic[6] = algebraic[4]/(constants[46]+algebraic[4])
    algebraic[48] = algebraic[31]*(1.00000-algebraic[45]/constants[51])
    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(7)*0.1
    if not iterable(voi):
        soln = fsolve(residualSN_0, initialGuess0, args=(algebraic, voi, constants, rates, states), xtol=1E-6)
        initialGuess0 = soln
        algebraic[36] = soln[0]
        algebraic[37] = soln[1]
        algebraic[38] = soln[2]
        algebraic[39] = soln[3]
        algebraic[40] = soln[4]
        algebraic[41] = soln[5]
        algebraic[42] = soln[6]
    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
            algebraic[36][i] = soln[0]
            algebraic[37][i] = soln[1]
            algebraic[38][i] = soln[2]
            algebraic[39][i] = soln[3]
            algebraic[40][i] = soln[4]
            algebraic[41][i] = soln[5]
            algebraic[42][i] = soln[6]

def residualSN_0(algebraicCandidate, algebraic, voi, constants, rates, states):
    resid = array([0.0] * 7)
    algebraic[36] = algebraicCandidate[0]
    algebraic[37] = algebraicCandidate[1]
    algebraic[38] = algebraicCandidate[2]
    algebraic[39] = algebraicCandidate[3]
    algebraic[40] = algebraicCandidate[4]
    algebraic[41] = algebraicCandidate[5]
    algebraic[42] = algebraicCandidate[6]
    resid[0] = (algebraic[36]-(constants[30]*states[2])/(1.00000+(1.00000/constants[32])*algebraic[37]))
    resid[1] = (algebraic[38]-(constants[33]*constants[84]*states[1])/((constants[33]*constants[84]*states[1])/constants[34]+constants[35]+(constants[36]/constants[37])*algebraic[37]))
    resid[2] = (algebraic[39]-algebraic[37]*(1.00000+algebraic[36]/constants[32]+algebraic[38]/constants[37]))
    resid[3] = (algebraic[40]-constants[38]*(1.00000-algebraic[39]/constants[85])*algebraic[3])
    resid[4] = (algebraic[41]-constants[39]*algebraic[34]*(1.00000-algebraic[39]/constants[85])*states[0])
    resid[5] = (algebraic[42]-(algebraic[40]+algebraic[41]))
    resid[6] = (algebraic[37]-(algebraic[42]+constants[29])/(constants[40]+(constants[41]/constants[32])*algebraic[36]+(constants[36]/constants[37])*algebraic[38]))
    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)