Generated Code
The following is matlab code generated by the CellML API from this CellML file. (Back to language selection)
The raw code is available.
function [VOI, STATES, ALGEBRAIC, CONSTANTS] = mainFunction()
% This is the "main function". In Matlab, things work best if you rename this function to match the filename.
[VOI, STATES, ALGEBRAIC, CONSTANTS] = solveModel();
end
function [algebraicVariableCount] = getAlgebraicVariableCount()
% Used later when setting a global variable with the number of algebraic variables.
% Note: This is not the "main method".
algebraicVariableCount =121;
end
% There are a total of 31 entries in each of the rate and state variable arrays.
% There are a total of 130 entries in the constant variable array.
%
function [VOI, STATES, ALGEBRAIC, CONSTANTS] = solveModel()
% Create ALGEBRAIC of correct size
global algebraicVariableCount; algebraicVariableCount = getAlgebraicVariableCount();
% Initialise constants and state variables
[INIT_STATES, CONSTANTS] = initConsts;
% Set timespan to solve over
tspan = [0, 10];
% Set numerical accuracy options for ODE solver
options = odeset('RelTol', 1e-06, 'AbsTol', 1e-06, 'MaxStep', 1);
% Solve model with ODE solver
[VOI, STATES] = ode15s(@(VOI, STATES)computeRates(VOI, STATES, CONSTANTS), tspan, INIT_STATES, options);
% Compute algebraic variables
[RATES, ALGEBRAIC] = computeRates(VOI, STATES, CONSTANTS);
ALGEBRAIC = computeAlgebraic(ALGEBRAIC, CONSTANTS, STATES, VOI);
% Plot state variables against variable of integration
[LEGEND_STATES, LEGEND_ALGEBRAIC, LEGEND_VOI, LEGEND_CONSTANTS] = createLegends();
figure();
plot(VOI, STATES);
xlabel(LEGEND_VOI);
l = legend(LEGEND_STATES);
set(l,'Interpreter','none');
end
function [LEGEND_STATES, LEGEND_ALGEBRAIC, LEGEND_VOI, LEGEND_CONSTANTS] = createLegends()
LEGEND_STATES = ''; LEGEND_ALGEBRAIC = ''; LEGEND_VOI = ''; LEGEND_CONSTANTS = '';
LEGEND_VOI = strpad('time in component environment (second)');
LEGEND_CONSTANTS(:,1) = strpad('protocol in component environment (dimensionless)');
LEGEND_CONSTANTS(:,2) = strpad('L_totmax in component b1_AR_Gs_parameters (uM)');
LEGEND_CONSTANTS(:,3) = strpad('sum_b1_AR in component b1_AR_Gs_parameters (uM)');
LEGEND_CONSTANTS(:,4) = strpad('Gs_tot in component b1_AR_Gs_parameters (uM)');
LEGEND_CONSTANTS(:,5) = strpad('Kl in component b1_AR_Gs_parameters (uM)');
LEGEND_CONSTANTS(:,6) = strpad('Kr in component b1_AR_Gs_parameters (uM)');
LEGEND_CONSTANTS(:,7) = strpad('Kc in component b1_AR_Gs_parameters (uM)');
LEGEND_CONSTANTS(:,8) = strpad('k_bar_kp in component b1_AR_Gs_parameters (per_sec)');
LEGEND_CONSTANTS(:,9) = strpad('k_bar_km in component b1_AR_Gs_parameters (per_sec)');
LEGEND_CONSTANTS(:,10) = strpad('k_p_kap in component b1_AR_Gs_parameters (per_uM_per_sec)');
LEGEND_CONSTANTS(:,11) = strpad('k_p_kam in component b1_AR_Gs_parameters (per_sec)');
LEGEND_CONSTANTS(:,12) = strpad('k_g_act in component b1_AR_Gs_parameters (per_sec)');
LEGEND_CONSTANTS(:,13) = strpad('k_hyd in component b1_AR_Gs_parameters (per_sec)');
LEGEND_CONSTANTS(:,14) = strpad('k_reassoc in component b1_AR_Gs_parameters (per_uM_per_sec)');
LEGEND_CONSTANTS(:,15) = strpad('AC_tot in component cAMP_parameters (uM)');
LEGEND_CONSTANTS(:,16) = strpad('ATP in component cAMP_parameters (uM)');
LEGEND_CONSTANTS(:,17) = strpad('PDE_tot in component cAMP_parameters (uM)');
LEGEND_CONSTANTS(:,18) = strpad('IBMX_tot in component cAMP_parameters (uM)');
LEGEND_CONSTANTS(:,19) = strpad('Fsk_tot in component cAMP_parameters (uM)');
LEGEND_CONSTANTS(:,20) = strpad('k_ac_basal in component cAMP_parameters (per_sec)');
LEGEND_CONSTANTS(:,21) = strpad('k_ac_gsa in component cAMP_parameters (per_sec)');
LEGEND_CONSTANTS(:,22) = strpad('k_ac_fsk in component cAMP_parameters (per_sec)');
LEGEND_CONSTANTS(:,23) = strpad('k_pde in component cAMP_parameters (per_sec)');
LEGEND_CONSTANTS(:,24) = strpad('Km_basal in component cAMP_parameters (uM)');
LEGEND_CONSTANTS(:,25) = strpad('Km_gsa in component cAMP_parameters (uM)');
LEGEND_CONSTANTS(:,26) = strpad('Km_fsk in component cAMP_parameters (uM)');
LEGEND_CONSTANTS(:,27) = strpad('Km_pde in component cAMP_parameters (uM)');
LEGEND_CONSTANTS(:,28) = strpad('K_gsa in component cAMP_parameters (uM)');
LEGEND_CONSTANTS(:,29) = strpad('K_fsk in component cAMP_parameters (uM)');
LEGEND_CONSTANTS(:,30) = strpad('Ki_ibmx in component cAMP_parameters (uM)');
LEGEND_CONSTANTS(:,31) = strpad('PKAI_tot in component PKA_parameters (uM)');
LEGEND_CONSTANTS(:,32) = strpad('PKAII_tot in component PKA_parameters (uM)');
LEGEND_CONSTANTS(:,33) = strpad('PKI_tot in component PKA_parameters (uM)');
LEGEND_CONSTANTS(:,34) = strpad('K_a in component PKA_parameters (uM)');
LEGEND_CONSTANTS(:,35) = strpad('K_b in component PKA_parameters (uM)');
LEGEND_CONSTANTS(:,36) = strpad('K_d in component PKA_parameters (uM)');
LEGEND_CONSTANTS(:,37) = strpad('Ki_pki in component PKA_parameters (uM)');
LEGEND_CONSTANTS(:,38) = strpad('PLB_tot in component PLB_parameters (uM)');
LEGEND_CONSTANTS(:,39) = strpad('PP1_tot in component PLB_parameters (uM)');
LEGEND_CONSTANTS(:,40) = strpad('Inhib1_tot in component PLB_parameters (uM)');
LEGEND_CONSTANTS(:,41) = strpad('k_pka_plb in component PLB_parameters (per_sec)');
LEGEND_CONSTANTS(:,42) = strpad('Km_pka_plb in component PLB_parameters (uM)');
LEGEND_CONSTANTS(:,43) = strpad('k_pp1_plb in component PLB_parameters (per_sec)');
LEGEND_CONSTANTS(:,44) = strpad('Km_pp1_plb in component PLB_parameters (uM)');
LEGEND_CONSTANTS(:,45) = strpad('k_pka_i1 in component PLB_parameters (per_sec)');
LEGEND_CONSTANTS(:,46) = strpad('Km_pka_i1 in component PLB_parameters (uM)');
LEGEND_CONSTANTS(:,47) = strpad('Vmax_pp2a_i1 in component PLB_parameters (uM_per_sec)');
LEGEND_CONSTANTS(:,48) = strpad('Km_pp2a_i1 in component PLB_parameters (uM)');
LEGEND_CONSTANTS(:,49) = strpad('Ki_inhib1 in component PLB_parameters (uM)');
LEGEND_CONSTANTS(:,50) = strpad('epsilon in component LCC_parameters (dimensionless)');
LEGEND_CONSTANTS(:,51) = strpad('LCC_tot in component LCC_parameters (uM)');
LEGEND_CONSTANTS(:,52) = strpad('PP1_lcc_tot in component LCC_parameters (uM)');
LEGEND_CONSTANTS(:,53) = strpad('PP2A_lcc_tot in component LCC_parameters (uM)');
LEGEND_CONSTANTS(:,54) = strpad('k_pka_lcc in component LCC_parameters (per_sec)');
LEGEND_CONSTANTS(:,55) = strpad('Km_pka_lcc in component LCC_parameters (uM)');
LEGEND_CONSTANTS(:,56) = strpad('k_pp1_lcc in component LCC_parameters (per_sec)');
LEGEND_CONSTANTS(:,57) = strpad('Km_pp1_lcc in component LCC_parameters (uM)');
LEGEND_CONSTANTS(:,58) = strpad('k_pp2a_lcc in component LCC_parameters (per_sec)');
LEGEND_CONSTANTS(:,59) = strpad('Km_pp2a_lcc in component LCC_parameters (uM)');
LEGEND_CONSTANTS(:,60) = strpad('V_myo in component EC_Coupling_Parameters (uL)');
LEGEND_CONSTANTS(:,61) = strpad('V_nsr in component EC_Coupling_Parameters (uL)');
LEGEND_CONSTANTS(:,62) = strpad('V_jsr in component EC_Coupling_Parameters (uL)');
LEGEND_CONSTANTS(:,63) = strpad('A_Cap in component EC_Coupling_Parameters (cm2)');
LEGEND_CONSTANTS(:,64) = strpad('Temp in component EC_Coupling_Parameters (kelvin)');
LEGEND_CONSTANTS(:,65) = strpad('Na_ext in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,66) = strpad('K_ext in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,67) = strpad('Ca_ext in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,68) = strpad('G_Na in component EC_Coupling_Parameters (mS_per_uF)');
LEGEND_CONSTANTS(:,69) = strpad('G_to in component EC_Coupling_Parameters (mS_per_uF)');
LEGEND_CONSTANTS(:,70) = strpad('G_ss in component EC_Coupling_Parameters (mS_per_uF)');
LEGEND_CONSTANTS(:,71) = strpad('G_Ki_bar in component EC_Coupling_Parameters (mS_per_uF)');
LEGEND_CONSTANTS(:,72) = strpad('G_Kp in component EC_Coupling_Parameters (mS_per_uF)');
LEGEND_CONSTANTS(:,73) = strpad('f in component EC_Coupling_Parameters (per_sec)');
LEGEND_CONSTANTS(:,74) = strpad('g in component EC_Coupling_Parameters (per_sec)');
LEGEND_CONSTANTS(:,75) = strpad('gamma_o in component EC_Coupling_Parameters (per_mM_per_sec)');
LEGEND_CONSTANTS(:,76) = strpad('omega in component EC_Coupling_Parameters (per_sec)');
LEGEND_CONSTANTS(:,77) = strpad('p_Ca in component EC_Coupling_Parameters (cm_per_sec)');
LEGEND_CONSTANTS(:,78) = strpad('p_K in component EC_Coupling_Parameters (cm_per_sec)');
LEGEND_CONSTANTS(:,79) = strpad('N_lcc in component EC_Coupling_Parameters (dimensionless)');
LEGEND_CONSTANTS(:,80) = strpad('I_Ca05 in component EC_Coupling_Parameters (uA_per_uF)');
LEGEND_CONSTANTS(:,81) = strpad('k_NaCa in component EC_Coupling_Parameters (uA_per_uF)');
LEGEND_CONSTANTS(:,82) = strpad('Km_Na in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,83) = strpad('Km_Ca in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,84) = strpad('k_sat in component EC_Coupling_Parameters (dimensionless)');
LEGEND_CONSTANTS(:,85) = strpad('eta in component EC_Coupling_Parameters (dimensionless)');
LEGEND_CONSTANTS(:,86) = strpad('I_bar_NaK in component EC_Coupling_Parameters (uA_per_uF)');
LEGEND_CONSTANTS(:,87) = strpad('Km_Nai in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,88) = strpad('Km_Ko in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,89) = strpad('I_bar_PCa in component EC_Coupling_Parameters (uA_per_uF)');
LEGEND_CONSTANTS(:,90) = strpad('Km_PCa in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,91) = strpad('G_CaB in component EC_Coupling_Parameters (uA_per_uF)');
LEGEND_CONSTANTS(:,92) = strpad('G_NaB in component EC_Coupling_Parameters (uA_per_uF)');
LEGEND_CONSTANTS(:,93) = strpad('Pns in component EC_Coupling_Parameters (dimensionless)');
LEGEND_CONSTANTS(:,94) = strpad('Km_NS in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,95) = strpad('I_up_bar in component EC_Coupling_Parameters (mM_per_sec)');
LEGEND_CONSTANTS(:,96) = strpad('Km_up0 in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,97) = strpad('NSR_bar in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,98) = strpad('tau_on in component EC_Coupling_Parameters (second)');
LEGEND_CONSTANTS(:,99) = strpad('tau_off in component EC_Coupling_Parameters (second)');
LEGEND_CONSTANTS(:,100) = strpad('G_max_rel in component EC_Coupling_Parameters (mM_per_sec)');
LEGEND_CONSTANTS(:,101) = strpad('d_Cai_th in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,102) = strpad('Km_rel in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,103) = strpad('CSQN_th in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,104) = strpad('CSQN_bar in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,105) = strpad('Km_CSQN in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,106) = strpad('tau_tr in component EC_Coupling_Parameters (second)');
LEGEND_CONSTANTS(:,107) = strpad('TRPN_bar in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,108) = strpad('CMDN_bar in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,109) = strpad('INDO_bar in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,110) = strpad('Km_TRPN in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,111) = strpad('Km_CMDN in component EC_Coupling_Parameters (mM)');
LEGEND_CONSTANTS(:,112) = strpad('Km_INDO in component EC_Coupling_Parameters (mM)');
LEGEND_ALGEBRAIC(:,62) = strpad('LR in component b1_AR_module (uM)');
LEGEND_ALGEBRAIC(:,63) = strpad('LRG in component b1_AR_module (uM)');
LEGEND_ALGEBRAIC(:,64) = strpad('RG in component b1_AR_module (uM)');
LEGEND_ALGEBRAIC(:,80) = strpad('BARK_DESENS in component b1_AR_module (uM_per_sec)');
LEGEND_ALGEBRAIC(:,12) = strpad('BARK_RESENS in component b1_AR_module (uM_per_sec)');
LEGEND_ALGEBRAIC(:,88) = strpad('PKA_DESENS in component b1_AR_module (uM_per_sec)');
LEGEND_ALGEBRAIC(:,23) = strpad('PKA_RESENS in component b1_AR_module (uM_per_sec)');
LEGEND_ALGEBRAIC(:,81) = strpad('G_ACT in component b1_AR_module (uM_per_sec)');
LEGEND_ALGEBRAIC(:,25) = strpad('HYD in component b1_AR_module (uM_per_sec)');
LEGEND_ALGEBRAIC(:,27) = strpad('REASSOC in component b1_AR_module (per_sec)');
LEGEND_ALGEBRAIC(:,65) = strpad('L in component b1_AR_module (uM)');
LEGEND_ALGEBRAIC(:,66) = strpad('R in component b1_AR_module (uM)');
LEGEND_ALGEBRAIC(:,67) = strpad('Gs in component b1_AR_module (uM)');
LEGEND_STATES(:,1) = strpad('b1_AR_tot in component b1_AR_module (uM)');
LEGEND_STATES(:,2) = strpad('b1_AR_d in component b1_AR_module (uM)');
LEGEND_STATES(:,3) = strpad('b1_AR_p in component b1_AR_module (uM)');
LEGEND_STATES(:,4) = strpad('Gs_agtp_tot in component b1_AR_module (uM)');
LEGEND_STATES(:,5) = strpad('Gs_agdp in component b1_AR_module (uM)');
LEGEND_STATES(:,6) = strpad('Gs_bg in component b1_AR_module (uM)');
LEGEND_ALGEBRAIC(:,68) = strpad('PKAC_I in component PKA_module (uM)');
LEGEND_ALGEBRAIC(:,69) = strpad('cAMP in component PKA_module (uM)');
LEGEND_ALGEBRAIC(:,42) = strpad('Gsa_GTP in component cAMP_module (uM)');
LEGEND_ALGEBRAIC(:,50) = strpad('Fsk in component cAMP_module (uM)');
LEGEND_ALGEBRAIC(:,43) = strpad('AC in component cAMP_module (uM)');
LEGEND_CONSTANTS(:,128) = strpad('PDE in component cAMP_module (uM)');
LEGEND_CONSTANTS(:,129) = strpad('IBMX in component cAMP_module (uM)');
LEGEND_STATES(:,7) = strpad('cAMP_tot in component cAMP_module (uM)');
LEGEND_ALGEBRAIC(:,44) = strpad('Gsa_GTP_AC in component cAMP_module (uM)');
LEGEND_ALGEBRAIC(:,51) = strpad('Fsk_AC in component cAMP_module (uM)');
LEGEND_ALGEBRAIC(:,48) = strpad('AC_ACT_GSA in component cAMP_module (uM)');
LEGEND_ALGEBRAIC(:,46) = strpad('AC_ACT_BASAL in component cAMP_module (uM)');
LEGEND_ALGEBRAIC(:,53) = strpad('AC_ACT_FSK in component cAMP_module (uM)');
LEGEND_ALGEBRAIC(:,82) = strpad('PDE_ACT in component cAMP_module (uM)');
LEGEND_CONSTANTS(:,130) = strpad('PDE_IBMX in component cAMP_module (uM)');
LEGEND_ALGEBRAIC(:,70) = strpad('PKI in component PKA_module (uM)');
LEGEND_ALGEBRAIC(:,71) = strpad('A2RC_I in component PKA_module (uM)');
LEGEND_ALGEBRAIC(:,72) = strpad('A2R_I in component PKA_module (uM)');
LEGEND_ALGEBRAIC(:,73) = strpad('A2RC_II in component PKA_module (uM)');
LEGEND_ALGEBRAIC(:,74) = strpad('A2R_II in component PKA_module (uM)');
LEGEND_ALGEBRAIC(:,75) = strpad('ARC_I in component PKA_module (uM)');
LEGEND_ALGEBRAIC(:,76) = strpad('ARC_II in component PKA_module (uM)');
LEGEND_ALGEBRAIC(:,77) = strpad('PKA_temp in component PKA_module (uM)');
LEGEND_ALGEBRAIC(:,78) = strpad('PKAC_II in component PKA_module (uM)');
LEGEND_ALGEBRAIC(:,29) = strpad('PLB in component PLB_module (uM)');
LEGEND_ALGEBRAIC(:,83) = strpad('PLB_PHOSPH in component PLB_module (uM_per_sec)');
LEGEND_ALGEBRAIC(:,60) = strpad('PLB_DEPHOSPH in component PLB_module (uM_per_sec)');
LEGEND_ALGEBRAIC(:,30) = strpad('Inhib1 in component PLB_module (uM)');
LEGEND_ALGEBRAIC(:,55) = strpad('Inhib1p_PP1 in component PLB_module (uM)');
LEGEND_ALGEBRAIC(:,84) = strpad('Inhib1_PHOSPH in component PLB_module (uM_per_sec)');
LEGEND_ALGEBRAIC(:,32) = strpad('Inhib1_DEPHOSPH in component PLB_module (uM_per_sec)');
LEGEND_STATES(:,8) = strpad('PLB_p in component PLB_module (uM)');
LEGEND_STATES(:,9) = strpad('Inhib1p_tot in component PLB_module (uM)');
LEGEND_ALGEBRAIC(:,56) = strpad('Inhib1p in component PLB_module (uM)');
LEGEND_ALGEBRAIC(:,57) = strpad('PP1 in component PLB_module (uM)');
LEGEND_ALGEBRAIC(:,1) = strpad('frac_PLB_p in component PLB_module (dimensionless)');
LEGEND_ALGEBRAIC(:,31) = strpad('frac_PLB in component PLB_module (dimensionless)');
LEGEND_CONSTANTS(:,121) = strpad('frac_PLB_o in component PLB_module (dimensionless)');
LEGEND_ALGEBRAIC(:,34) = strpad('LCCa in component LCC_module (uM)');
LEGEND_ALGEBRAIC(:,85) = strpad('LCCa_PHOSPH in component LCC_module (uM_per_sec)');
LEGEND_ALGEBRAIC(:,36) = strpad('LCCa_DEPHOSPH in component LCC_module (uM_per_sec)');
LEGEND_ALGEBRAIC(:,38) = strpad('LCCb in component LCC_module (uM)');
LEGEND_ALGEBRAIC(:,86) = strpad('LCCb_PHOSPH in component LCC_module (uM_per_sec)');
LEGEND_ALGEBRAIC(:,40) = strpad('LCCb_DEPHOSPH in component LCC_module (uM_per_sec)');
LEGEND_STATES(:,10) = strpad('LCCa_p in component LCC_module (uM)');
LEGEND_STATES(:,11) = strpad('LCCb_p in component LCC_module (uM)');
LEGEND_ALGEBRAIC(:,2) = strpad('frac_LCCa_p in component LCC_module (dimensionless)');
LEGEND_CONSTANTS(:,122) = strpad('frac_LCCa_po in component LCC_module (dimensionless)');
LEGEND_ALGEBRAIC(:,33) = strpad('frac_LCCb_p in component LCC_module (dimensionless)');
LEGEND_CONSTANTS(:,123) = strpad('frac_LCCb_po in component LCC_module (dimensionless)');
LEGEND_ALGEBRAIC(:,35) = strpad('E_Na in component Nernst_Potentials (mV)');
LEGEND_ALGEBRAIC(:,37) = strpad('E_K in component Nernst_Potentials (mV)');
LEGEND_ALGEBRAIC(:,39) = strpad('E_Ca in component Nernst_Potentials (mV)');
LEGEND_CONSTANTS(:,124) = strpad('E_Cl in component Nernst_Potentials (mV)');
LEGEND_CONSTANTS(:,113) = strpad('R in component Nernst_Potentials (joules_per_mole_kelvin)');
LEGEND_CONSTANTS(:,114) = strpad('Frdy in component Nernst_Potentials (coulombs_per_mole)');
LEGEND_CONSTANTS(:,125) = strpad('FoRT in component Nernst_Potentials (per_mV)');
LEGEND_CONSTANTS(:,115) = strpad('z_Na in component Nernst_Potentials (dimensionless)');
LEGEND_CONSTANTS(:,116) = strpad('z_K in component Nernst_Potentials (dimensionless)');
LEGEND_CONSTANTS(:,117) = strpad('z_Ca in component Nernst_Potentials (dimensionless)');
LEGEND_STATES(:,12) = strpad('Na_i in component Ion_Concentrations_and_Membrane_Potential (mM)');
LEGEND_STATES(:,13) = strpad('K_i in component Ion_Concentrations_and_Membrane_Potential (mM)');
LEGEND_STATES(:,14) = strpad('Ca_i in component Ion_Concentrations_and_Membrane_Potential (mM)');
LEGEND_ALGEBRAIC(:,3) = strpad('am in component Fast_Na_Current (per_sec)');
LEGEND_ALGEBRAIC(:,13) = strpad('bm in component Fast_Na_Current (per_sec)');
LEGEND_ALGEBRAIC(:,4) = strpad('ah in component Fast_Na_Current (per_sec)');
LEGEND_ALGEBRAIC(:,5) = strpad('aj in component Fast_Na_Current (per_sec)');
LEGEND_ALGEBRAIC(:,14) = strpad('bh in component Fast_Na_Current (per_sec)');
LEGEND_ALGEBRAIC(:,15) = strpad('bj in component Fast_Na_Current (per_sec)');
LEGEND_STATES(:,15) = strpad('m in component Fast_Na_Current (dimensionless)');
LEGEND_STATES(:,16) = strpad('h in component Fast_Na_Current (dimensionless)');
LEGEND_STATES(:,17) = strpad('j in component Fast_Na_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,41) = strpad('I_Na in component Fast_Na_Current (uA_per_uF)');
LEGEND_STATES(:,18) = strpad('V_m in component Ion_Concentrations_and_Membrane_Potential (mV)');
LEGEND_ALGEBRAIC(:,6) = strpad('a_lcc in component L_Type_Calcium_Current (per_sec)');
LEGEND_ALGEBRAIC(:,16) = strpad('b_lcc in component L_Type_Calcium_Current (per_sec)');
LEGEND_ALGEBRAIC(:,18) = strpad('f_lcc in component L_Type_Calcium_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,7) = strpad('y_lcc_inf in component L_Type_Calcium_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,17) = strpad('tau_y_lcc in component L_Type_Calcium_Current (second)');
LEGEND_ALGEBRAIC(:,24) = strpad('gamma in component L_Type_Calcium_Current (per_mM_per_sec)');
LEGEND_ALGEBRAIC(:,26) = strpad('v_gamma in component L_Type_Calcium_Current (per_mM_per_sec)');
LEGEND_ALGEBRAIC(:,28) = strpad('v_omega in component L_Type_Calcium_Current (dimensionless)');
LEGEND_STATES(:,19) = strpad('v in component L_Type_Calcium_Current (dimensionless)');
LEGEND_STATES(:,20) = strpad('w in component L_Type_Calcium_Current (dimensionless)');
LEGEND_STATES(:,21) = strpad('x in component L_Type_Calcium_Current (dimensionless)');
LEGEND_STATES(:,22) = strpad('y in component L_Type_Calcium_Current (dimensionless)');
LEGEND_STATES(:,23) = strpad('z in component L_Type_Calcium_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,45) = strpad('i_bar_Ca in component L_Type_Calcium_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,47) = strpad('i_bar_K in component L_Type_Calcium_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,49) = strpad('f_avail in component L_Type_Calcium_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,52) = strpad('I_Ca in component L_Type_Calcium_Current (uA_per_uF)');
LEGEND_ALGEBRAIC(:,54) = strpad('I_CaK in component L_Type_Calcium_Current (uA_per_uF)');
LEGEND_ALGEBRAIC(:,58) = strpad('I_Ca_tot in component L_Type_Calcium_Current (uA_per_uF)');
LEGEND_ALGEBRAIC(:,8) = strpad('r_toss in component Transient_Outward_K_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,9) = strpad('s_toss in component Transient_Outward_K_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,19) = strpad('tau_r_to in component Transient_Outward_K_Current (second)');
LEGEND_ALGEBRAIC(:,20) = strpad('tau_s_to in component Transient_Outward_K_Current (second)');
LEGEND_ALGEBRAIC(:,21) = strpad('tau_ss_to in component Transient_Outward_K_Current (second)');
LEGEND_STATES(:,24) = strpad('r_to in component Transient_Outward_K_Current (dimensionless)');
LEGEND_STATES(:,25) = strpad('s_to in component Transient_Outward_K_Current (dimensionless)');
LEGEND_STATES(:,26) = strpad('ss_to in component Transient_Outward_K_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,59) = strpad('I_to in component Transient_Outward_K_Current (uA_per_uF)');
LEGEND_ALGEBRAIC(:,10) = strpad('r_ss_inf in component Steady_State_K_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,22) = strpad('tau_r_ss in component Steady_State_K_Current (second)');
LEGEND_ALGEBRAIC(:,11) = strpad('s_ss_inf in component Steady_State_K_Current (dimensionless)');
LEGEND_CONSTANTS(:,126) = strpad('tau_s_ss in component Steady_State_K_Current (second)');
LEGEND_STATES(:,27) = strpad('r_ss in component Steady_State_K_Current (dimensionless)');
LEGEND_STATES(:,28) = strpad('s_ss in component Steady_State_K_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,61) = strpad('I_ss in component Steady_State_K_Current (uA_per_uF)');
LEGEND_ALGEBRAIC(:,79) = strpad('a_Ki in component Time_Independent_K_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,87) = strpad('b_Ki in component Time_Independent_K_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,89) = strpad('Ki_ss in component Time_Independent_K_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,90) = strpad('I_Ki in component Time_Independent_K_Current (uA_per_uF)');
LEGEND_ALGEBRAIC(:,91) = strpad('Kp in component Plateau_K_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,92) = strpad('I_Kp in component Plateau_K_Current (uA_per_uF)');
LEGEND_ALGEBRAIC(:,93) = strpad('s4 in component Na_Ca_Exchanger_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,94) = strpad('s5 in component Na_Ca_Exchanger_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,95) = strpad('I_NCX in component Na_Ca_Exchanger_Current (uA_per_uF)');
LEGEND_CONSTANTS(:,127) = strpad('sigma in component Na_K_Pump_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,96) = strpad('f_NaK in component Na_K_Pump_Current (dimensionless)');
LEGEND_ALGEBRAIC(:,97) = strpad('I_NaK in component Na_K_Pump_Current (uA_per_uF)');
LEGEND_ALGEBRAIC(:,98) = strpad('I_PCa in component Sarcolemmal_Ca_Pump_Current (uA_per_uF)');
LEGEND_ALGEBRAIC(:,99) = strpad('I_CaB in component Ca_Background_Current (uA_per_uF)');
LEGEND_ALGEBRAIC(:,100) = strpad('I_NaB in component Na_Background_Current (uA_per_uF)');
LEGEND_ALGEBRAIC(:,101) = strpad('I_Na_tot in component Total_Membrane_Currents (uA_per_uF)');
LEGEND_ALGEBRAIC(:,102) = strpad('I_K_tot in component Total_Membrane_Currents (uA_per_uF)');
LEGEND_ALGEBRAIC(:,103) = strpad('I_Ca_tot in component Total_Membrane_Currents (uA_per_uF)');
LEGEND_ALGEBRAIC(:,104) = strpad('t_rel in component Calcium_Induced_Calcium_Release (second)');
LEGEND_ALGEBRAIC(:,106) = strpad('ryr_on in component Calcium_Induced_Calcium_Release (dimensionless)');
LEGEND_ALGEBRAIC(:,108) = strpad('ryr_off in component Calcium_Induced_Calcium_Release (dimensionless)');
LEGEND_ALGEBRAIC(:,110) = strpad('g_rel in component Calcium_Induced_Calcium_Release (per_sec)');
LEGEND_ALGEBRAIC(:,111) = strpad('I_rel in component Calcium_Induced_Calcium_Release (mM_per_sec)');
LEGEND_STATES(:,29) = strpad('Ca_jsr in component Other_SR_Fluxes_and_Concentrations (mM)');
LEGEND_STATES(:,30) = strpad('trel in component Ion_Concentrations_and_Membrane_Potential (second)');
LEGEND_ALGEBRAIC(:,112) = strpad('Km_up in component Other_SR_Fluxes_and_Concentrations (mM)');
LEGEND_ALGEBRAIC(:,113) = strpad('I_up in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)');
LEGEND_ALGEBRAIC(:,114) = strpad('I_leak in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)');
LEGEND_ALGEBRAIC(:,115) = strpad('I_tr in component Other_SR_Fluxes_and_Concentrations (mM_per_sec)');
LEGEND_ALGEBRAIC(:,117) = strpad('B_jsr in component Other_SR_Fluxes_and_Concentrations (dimensionless)');
LEGEND_STATES(:,31) = strpad('Ca_nsr in component Other_SR_Fluxes_and_Concentrations (mM)');
LEGEND_ALGEBRAIC(:,119) = strpad('SR_content in component Other_SR_Fluxes_and_Concentrations (uM)');
LEGEND_ALGEBRAIC(:,116) = strpad('b_trpn in component Cytoplasmic_Calcium_Buffering (dimensionless)');
LEGEND_ALGEBRAIC(:,118) = strpad('b_cmdn in component Cytoplasmic_Calcium_Buffering (dimensionless)');
LEGEND_ALGEBRAIC(:,120) = strpad('b_indo in component Cytoplasmic_Calcium_Buffering (dimensionless)');
LEGEND_ALGEBRAIC(:,121) = strpad('B_myo in component Cytoplasmic_Calcium_Buffering (dimensionless)');
LEGEND_ALGEBRAIC(:,109) = strpad('I_app in component Ion_Concentrations_and_Membrane_Potential (uA_per_uF)');
LEGEND_CONSTANTS(:,118) = strpad('V_hold in component Ion_Concentrations_and_Membrane_Potential (mV)');
LEGEND_CONSTANTS(:,119) = strpad('V_test in component Ion_Concentrations_and_Membrane_Potential (mV)');
LEGEND_ALGEBRAIC(:,105) = strpad('V_clamp in component Ion_Concentrations_and_Membrane_Potential (mV)');
LEGEND_CONSTANTS(:,120) = strpad('R_clamp in component Ion_Concentrations_and_Membrane_Potential (dimensionless)');
LEGEND_ALGEBRAIC(:,107) = strpad('I_pace in component Ion_Concentrations_and_Membrane_Potential (uA_per_uF)');
LEGEND_RATES(:,1) = strpad('d/dt b1_AR_tot in component b1_AR_module (uM)');
LEGEND_RATES(:,2) = strpad('d/dt b1_AR_d in component b1_AR_module (uM)');
LEGEND_RATES(:,3) = strpad('d/dt b1_AR_p in component b1_AR_module (uM)');
LEGEND_RATES(:,4) = strpad('d/dt Gs_agtp_tot in component b1_AR_module (uM)');
LEGEND_RATES(:,5) = strpad('d/dt Gs_agdp in component b1_AR_module (uM)');
LEGEND_RATES(:,6) = strpad('d/dt Gs_bg in component b1_AR_module (uM)');
LEGEND_RATES(:,7) = strpad('d/dt cAMP_tot in component cAMP_module (uM)');
LEGEND_RATES(:,8) = strpad('d/dt PLB_p in component PLB_module (uM)');
LEGEND_RATES(:,9) = strpad('d/dt Inhib1p_tot in component PLB_module (uM)');
LEGEND_RATES(:,10) = strpad('d/dt LCCa_p in component LCC_module (uM)');
LEGEND_RATES(:,11) = strpad('d/dt LCCb_p in component LCC_module (uM)');
LEGEND_RATES(:,15) = strpad('d/dt m in component Fast_Na_Current (dimensionless)');
LEGEND_RATES(:,16) = strpad('d/dt h in component Fast_Na_Current (dimensionless)');
LEGEND_RATES(:,17) = strpad('d/dt j in component Fast_Na_Current (dimensionless)');
LEGEND_RATES(:,19) = strpad('d/dt v in component L_Type_Calcium_Current (dimensionless)');
LEGEND_RATES(:,20) = strpad('d/dt w in component L_Type_Calcium_Current (dimensionless)');
LEGEND_RATES(:,21) = strpad('d/dt x in component L_Type_Calcium_Current (dimensionless)');
LEGEND_RATES(:,22) = strpad('d/dt y in component L_Type_Calcium_Current (dimensionless)');
LEGEND_RATES(:,23) = strpad('d/dt z in component L_Type_Calcium_Current (dimensionless)');
LEGEND_RATES(:,24) = strpad('d/dt r_to in component Transient_Outward_K_Current (dimensionless)');
LEGEND_RATES(:,25) = strpad('d/dt s_to in component Transient_Outward_K_Current (dimensionless)');
LEGEND_RATES(:,26) = strpad('d/dt ss_to in component Transient_Outward_K_Current (dimensionless)');
LEGEND_RATES(:,27) = strpad('d/dt r_ss in component Steady_State_K_Current (dimensionless)');
LEGEND_RATES(:,28) = strpad('d/dt s_ss in component Steady_State_K_Current (dimensionless)');
LEGEND_RATES(:,31) = strpad('d/dt Ca_nsr in component Other_SR_Fluxes_and_Concentrations (mM)');
LEGEND_RATES(:,29) = strpad('d/dt Ca_jsr in component Other_SR_Fluxes_and_Concentrations (mM)');
LEGEND_RATES(:,12) = strpad('d/dt Na_i in component Ion_Concentrations_and_Membrane_Potential (mM)');
LEGEND_RATES(:,13) = strpad('d/dt K_i in component Ion_Concentrations_and_Membrane_Potential (mM)');
LEGEND_RATES(:,14) = strpad('d/dt Ca_i in component Ion_Concentrations_and_Membrane_Potential (mM)');
LEGEND_RATES(:,18) = strpad('d/dt V_m in component Ion_Concentrations_and_Membrane_Potential (mV)');
LEGEND_RATES(:,30) = strpad('d/dt trel in component Ion_Concentrations_and_Membrane_Potential (second)');
LEGEND_STATES = LEGEND_STATES';
LEGEND_ALGEBRAIC = LEGEND_ALGEBRAIC';
LEGEND_RATES = LEGEND_RATES';
LEGEND_CONSTANTS = LEGEND_CONSTANTS';
end
function [STATES, CONSTANTS] = initConsts()
VOI = 0; CONSTANTS = []; STATES = []; ALGEBRAIC = [];
CONSTANTS(:,1) = 0;
CONSTANTS(:,2) = 1;
CONSTANTS(:,3) = 0.0132;
CONSTANTS(:,4) = 3.83;
CONSTANTS(:,5) = 0.285;
CONSTANTS(:,6) = 0.062;
CONSTANTS(:,7) = 33;
CONSTANTS(:,8) = 1.1e-3;
CONSTANTS(:,9) = 2.2e-3;
CONSTANTS(:,10) = 3.6e-3;
CONSTANTS(:,11) = 2.2e-3;
CONSTANTS(:,12) = 16;
CONSTANTS(:,13) = 0.8;
CONSTANTS(:,14) = 1.21e3;
CONSTANTS(:,15) = 49.7e-3;
CONSTANTS(:,16) = 5e3;
CONSTANTS(:,17) = 38.9e-3;
CONSTANTS(:,18) = 0;
CONSTANTS(:,19) = 0;
CONSTANTS(:,20) = 0.2;
CONSTANTS(:,21) = 8.5;
CONSTANTS(:,22) = 7.3;
CONSTANTS(:,23) = 5;
CONSTANTS(:,24) = 1.03e3;
CONSTANTS(:,25) = 315;
CONSTANTS(:,26) = 860;
CONSTANTS(:,27) = 1.3;
CONSTANTS(:,28) = 0.4;
CONSTANTS(:,29) = 44;
CONSTANTS(:,30) = 30;
CONSTANTS(:,31) = 0.59;
CONSTANTS(:,32) = 0.025;
CONSTANTS(:,33) = 0.18;
CONSTANTS(:,34) = 9.14;
CONSTANTS(:,35) = 1.64;
CONSTANTS(:,36) = 4.375;
CONSTANTS(:,37) = 0.2e-3;
CONSTANTS(:,38) = 106;
CONSTANTS(:,39) = 0.89;
CONSTANTS(:,40) = 0.3;
CONSTANTS(:,41) = 54;
CONSTANTS(:,42) = 21;
CONSTANTS(:,43) = 8.5;
CONSTANTS(:,44) = 7;
CONSTANTS(:,45) = 60;
CONSTANTS(:,46) = 1;
CONSTANTS(:,47) = 14;
CONSTANTS(:,48) = 1;
CONSTANTS(:,49) = 1e-3;
CONSTANTS(:,50) = 10;
CONSTANTS(:,51) = 25e-3;
CONSTANTS(:,52) = 25e-3;
CONSTANTS(:,53) = 25e-3;
CONSTANTS(:,54) = 54;
CONSTANTS(:,55) = 21;
CONSTANTS(:,56) = 8.52;
CONSTANTS(:,57) = 3;
CONSTANTS(:,58) = 10.1;
CONSTANTS(:,59) = 3;
CONSTANTS(:,60) = 20.8e-6;
CONSTANTS(:,61) = 9.88e-7;
CONSTANTS(:,62) = 9.3e-8;
CONSTANTS(:,63) = 1.534e-4;
CONSTANTS(:,64) = 310;
CONSTANTS(:,65) = 140;
CONSTANTS(:,66) = 5.4;
CONSTANTS(:,67) = 1.8;
CONSTANTS(:,68) = 8;
CONSTANTS(:,69) = 0.35;
CONSTANTS(:,70) = 0.07;
CONSTANTS(:,71) = 0.24;
CONSTANTS(:,72) = 0.008;
CONSTANTS(:,73) = 300;
CONSTANTS(:,74) = 2e3;
CONSTANTS(:,75) = 5187.5;
CONSTANTS(:,76) = 10;
CONSTANTS(:,77) = 1.7469e-8;
CONSTANTS(:,78) = 3.234e-11;
CONSTANTS(:,79) = 3e5;
CONSTANTS(:,80) = -0.458;
CONSTANTS(:,81) = 1483;
CONSTANTS(:,82) = 87.5;
CONSTANTS(:,83) = 1.38;
CONSTANTS(:,84) = 0.1;
CONSTANTS(:,85) = 0.35;
CONSTANTS(:,86) = 1.1;
CONSTANTS(:,87) = 10;
CONSTANTS(:,88) = 1.5;
CONSTANTS(:,89) = 1.15;
CONSTANTS(:,90) = 0.5e-3;
CONSTANTS(:,91) = 2.8e-3;
CONSTANTS(:,92) = 1.18e-3;
CONSTANTS(:,93) = 0;
CONSTANTS(:,94) = 1.2e-3;
CONSTANTS(:,95) = 4.7;
CONSTANTS(:,96) = 3e-4;
CONSTANTS(:,97) = 15;
CONSTANTS(:,98) = 2e-3;
CONSTANTS(:,99) = 2e-3;
CONSTANTS(:,100) = 60e3;
CONSTANTS(:,101) = 0.18e-3;
CONSTANTS(:,102) = 0.8e-3;
CONSTANTS(:,103) = 8.75;
CONSTANTS(:,104) = 15;
CONSTANTS(:,105) = 0.8;
CONSTANTS(:,106) = 5.7e-4;
CONSTANTS(:,107) = 0.07;
CONSTANTS(:,108) = 0.05;
CONSTANTS(:,109) = 0.07;
CONSTANTS(:,110) = 0.5128e-3;
CONSTANTS(:,111) = 2.38e-3;
CONSTANTS(:,112) = 8.44e-4;
STATES(:,1) = 0.01205;
STATES(:,2) = 0;
STATES(:,3) = 1.154e-3;
STATES(:,4) = 0.02505;
STATES(:,5) = 6.446e-4;
STATES(:,6) = 0.02569;
STATES(:,7) = 0.8453;
STATES(:,8) = 4.105;
STATES(:,9) = 0.0526;
STATES(:,10) = 5.103e-3;
STATES(:,11) = 5.841e-3;
CONSTANTS(:,113) = 8314;
CONSTANTS(:,114) = 96485;
CONSTANTS(:,115) = 1;
CONSTANTS(:,116) = 1;
CONSTANTS(:,117) = 2;
STATES(:,12) = 16;
STATES(:,13) = 145;
STATES(:,14) = 1.58e-4;
STATES(:,15) = 1.4e-3;
STATES(:,16) = 0.99;
STATES(:,17) = 0.99;
STATES(:,18) = -85.66;
STATES(:,19) = 0;
STATES(:,20) = 0;
STATES(:,21) = 0.13;
STATES(:,22) = 0.96;
STATES(:,23) = 0.92;
STATES(:,24) = 1.4e-3;
STATES(:,25) = 1;
STATES(:,26) = 0.613;
STATES(:,27) = 198e-3;
STATES(:,28) = 0.43;
STATES(:,29) = 1.92;
STATES(:,30) = 0.9;
STATES(:,31) = 1.92;
CONSTANTS(:,118) = -40;
CONSTANTS(:,119) = -10;
CONSTANTS(:,120) = 0.02;
CONSTANTS(:,121) = 0.961300;
CONSTANTS(:,122) = 0.204100;
CONSTANTS(:,123) = 0.233600;
CONSTANTS(:,124) = - 40.0000;
CONSTANTS(:,125) = (CONSTANTS(:,114)./CONSTANTS(:,113))./CONSTANTS(:,64);
CONSTANTS(:,126) = 2.10000;
CONSTANTS(:,127) = (exp(CONSTANTS(:,65)./67.3000) - 1.00000)./7.00000;
[CONSTANTS, STATES, ALGEBRAIC] = rootfind_0(VOI, CONSTANTS, STATES, ALGEBRAIC);
if (isempty(STATES)), warning('Initial values for states not set');, end
end
function [RATES, ALGEBRAIC] = computeRates(VOI, STATES, CONSTANTS)
global algebraicVariableCount;
statesSize = size(STATES);
statesColumnCount = statesSize(2);
if ( statesColumnCount == 1)
STATES = STATES';
ALGEBRAIC = zeros(1, algebraicVariableCount);
utilOnes = 1;
else
statesRowCount = statesSize(1);
ALGEBRAIC = zeros(statesRowCount, algebraicVariableCount);
RATES = zeros(statesRowCount, statesColumnCount);
utilOnes = ones(statesRowCount, 1);
end
ALGEBRAIC(:,11) = 1.00000./(1.00000+exp((STATES(:,18)+87.5000)./10.3000));
RATES(:,28) = (ALGEBRAIC(:,11) - STATES(:,28))./CONSTANTS(:,126);
ALGEBRAIC(:,3) = ( 0.320000.*(STATES(:,18)+47.1300))./(1.00000 - exp( - 0.100000.*(STATES(:,18)+47.1300)));
ALGEBRAIC(:,13) = 0.0800000.*exp( - STATES(:,18)./11.0000);
RATES(:,15) = 1000.00.*( ALGEBRAIC(:,3).*(1.00000 - STATES(:,15)) - ALGEBRAIC(:,13).*STATES(:,15));
ALGEBRAIC(:,4) = piecewise({STATES(:,18)>= - 40.0000, 0.00000 }, 0.135000.*exp((80.0000+STATES(:,18))./ - 6.80000));
ALGEBRAIC(:,14) = piecewise({STATES(:,18)>= - 40.0000, 1.00000./( 0.130000.*(1.00000+exp( - (STATES(:,18)+10.6600)./11.1000))) }, 3.56000.*exp( 0.0790000.*STATES(:,18))+ 310000..*exp( 0.350000.*STATES(:,18)));
RATES(:,16) = 1000.00.*( ALGEBRAIC(:,4).*(1.00000 - STATES(:,16)) - ALGEBRAIC(:,14).*STATES(:,16));
ALGEBRAIC(:,5) = piecewise({STATES(:,18)>= - 40.0000, 0.00000 }, ( ( - 127140..*exp( 0.244400.*STATES(:,18)) - 3.47400e-05.*exp( - 0.0439100.*STATES(:,18))).*(STATES(:,18)+37.7800))./(1.00000+exp( 0.311000.*(STATES(:,18)+79.2300))));
ALGEBRAIC(:,15) = piecewise({STATES(:,18)>= - 40.0000, ( 0.300000.*exp( - 2.57500e-07.*STATES(:,18)))./(1.00000+exp( - 0.100000.*(STATES(:,18)+32.0000))) }, ( 0.121200.*exp( - 0.0105200.*STATES(:,18)))./(1.00000+exp( - 0.137800.*(STATES(:,18)+40.1400))));
RATES(:,17) = 1000.00.*( ALGEBRAIC(:,5).*(1.00000 - STATES(:,17)) - ALGEBRAIC(:,15).*STATES(:,17));
ALGEBRAIC(:,6) = 400.000.*exp((STATES(:,18)+2.00000)./10.0000);
ALGEBRAIC(:,16) = 50.0000.*exp(( - 1.00000.*(STATES(:,18)+2.00000))./13.0000);
RATES(:,19) = ALGEBRAIC(:,6).*(1.00000 - STATES(:,19)) - ALGEBRAIC(:,16).*STATES(:,19);
RATES(:,20) = 2.00000.*ALGEBRAIC(:,6).*(1.00000 - STATES(:,20)) - (ALGEBRAIC(:,16)./2.00000).*STATES(:,20);
ALGEBRAIC(:,2) = STATES(:,10)./CONSTANTS(:,51);
ALGEBRAIC(:,18) = CONSTANTS(:,73).*(( 0.375000.*ALGEBRAIC(:,2))./CONSTANTS(:,122)+0.625000);
RATES(:,21) = ALGEBRAIC(:,18).*(1.00000 - STATES(:,21)) - CONSTANTS(:,74).*STATES(:,21);
ALGEBRAIC(:,7) = 1.00000./(1.00000+exp((STATES(:,18)+55.0000)./7.50000))+0.100000./(1.00000+exp(( - STATES(:,18)+21.0000)./6.00000));
ALGEBRAIC(:,17) = 0.0200000+0.300000./(1.00000+exp((STATES(:,18)+30.0000)./9.50000));
RATES(:,22) = (ALGEBRAIC(:,7) - STATES(:,22))./ALGEBRAIC(:,17);
ALGEBRAIC(:,8) = 1.00000./(1.00000+exp((STATES(:,18)+10.6000)./ - 11.4200));
ALGEBRAIC(:,19) = 1.00000./( 45.1600.*exp( 0.0357700.*(STATES(:,18)+50.0000))+ 98.9000.*exp( - 0.100000.*(STATES(:,18)+38.0000)));
RATES(:,24) = (ALGEBRAIC(:,8) - STATES(:,24))./ALGEBRAIC(:,19);
ALGEBRAIC(:,9) = 1.00000./(1.00000+exp((STATES(:,18)+43.5000)./6.88410));
ALGEBRAIC(:,20) = 0.350000.*exp( - 1.00000.*power((STATES(:,18)+70.0000)./15.0000, 2.00000))+0.0350000;
RATES(:,25) = (ALGEBRAIC(:,9) - STATES(:,25))./ALGEBRAIC(:,20);
ALGEBRAIC(:,21) = 3.70000.*exp( - 1.00000.*power((STATES(:,18)+70.0000)./30.0000, 2.00000))+0.0350000;
RATES(:,26) = (ALGEBRAIC(:,9) - STATES(:,26))./ALGEBRAIC(:,21);
ALGEBRAIC(:,10) = 1.00000./(1.00000+exp((STATES(:,18)+11.5000)./ - 11.8200));
ALGEBRAIC(:,22) = 10.0000./( 45.1600.*exp( 0.0357700.*(STATES(:,18)+50.0000))+ 98.9000.*exp( - 0.100000.*(STATES(:,18)+38.0000)));
RATES(:,27) = (ALGEBRAIC(:,10) - STATES(:,27))./ALGEBRAIC(:,22);
ALGEBRAIC(:,25) = CONSTANTS(:,13).*STATES(:,4);
ALGEBRAIC(:,27) = CONSTANTS(:,14).*STATES(:,5).*STATES(:,6);
RATES(:,5) = ALGEBRAIC(:,25) - ALGEBRAIC(:,27);
ALGEBRAIC(:,24) = CONSTANTS(:,75).*STATES(:,14);
ALGEBRAIC(:,26) = ALGEBRAIC(:,24).*(power(1.00000 - STATES(:,19), 4.00000)+ 2.00000.*STATES(:,19).*power(1.00000 - STATES(:,19), 3.00000)+ 4.00000.*power(STATES(:,19), 2.00000).*power(1.00000 - STATES(:,19), 2.00000)+ 8.00000.*power(STATES(:,19), 3.00000).*(1.00000 - STATES(:,19))+ 16.0000.*power(STATES(:,19), 4.00000).*(1.00000 - ALGEBRAIC(:,18)./CONSTANTS(:,74)));
ALGEBRAIC(:,28) = CONSTANTS(:,76).*(power(1.00000 - STATES(:,20), 4.00000)+ 0.500000.*STATES(:,20).*power(1.00000 - STATES(:,20), 3.00000)+ 0.250000.*power(STATES(:,20), 2.00000).*power(1.00000 - STATES(:,20), 2.00000)+ 0.125000.*power(STATES(:,20), 3.00000).*(1.00000 - STATES(:,20))+ 0.0625000.*power(STATES(:,20), 4.00000));
RATES(:,23) = ALGEBRAIC(:,28).*(1.00000 - STATES(:,23)) - ALGEBRAIC(:,26).*STATES(:,23);
[CONSTANTS, STATES, ALGEBRAIC] = rootfind_1(VOI, CONSTANTS, STATES, ALGEBRAIC);
ALGEBRAIC(:,80) = CONSTANTS(:,8).*(ALGEBRAIC(:,62)+ALGEBRAIC(:,63));
ALGEBRAIC(:,12) = CONSTANTS(:,9).*STATES(:,2);
RATES(:,2) = ALGEBRAIC(:,80) - ALGEBRAIC(:,12);
ALGEBRAIC(:,81) = CONSTANTS(:,12).*(ALGEBRAIC(:,64)+ALGEBRAIC(:,63));
RATES(:,4) = ALGEBRAIC(:,81) - ALGEBRAIC(:,25);
RATES(:,6) = ALGEBRAIC(:,81) - ALGEBRAIC(:,27);
[CONSTANTS, STATES, ALGEBRAIC] = rootfind_2(VOI, CONSTANTS, STATES, ALGEBRAIC);
ALGEBRAIC(:,48) = ( CONSTANTS(:,21).*ALGEBRAIC(:,44).*CONSTANTS(:,16))./(CONSTANTS(:,25)+CONSTANTS(:,16));
ALGEBRAIC(:,46) = ( CONSTANTS(:,20).*ALGEBRAIC(:,43).*CONSTANTS(:,16))./(CONSTANTS(:,24)+CONSTANTS(:,16));
[CONSTANTS, STATES, ALGEBRAIC] = rootfind_3(VOI, CONSTANTS, STATES, ALGEBRAIC);
ALGEBRAIC(:,53) = ( CONSTANTS(:,22).*ALGEBRAIC(:,51).*CONSTANTS(:,16))./(CONSTANTS(:,26)+CONSTANTS(:,16));
ALGEBRAIC(:,82) = ( CONSTANTS(:,23).*CONSTANTS(:,128).*ALGEBRAIC(:,69))./(CONSTANTS(:,27)+ALGEBRAIC(:,69));
RATES(:,7) = (ALGEBRAIC(:,46)+ALGEBRAIC(:,48)+ALGEBRAIC(:,53)) - ALGEBRAIC(:,82);
ALGEBRAIC(:,29) = CONSTANTS(:,38) - STATES(:,8);
ALGEBRAIC(:,83) = ( CONSTANTS(:,41).*ALGEBRAIC(:,68).*ALGEBRAIC(:,29))./(CONSTANTS(:,42)+ALGEBRAIC(:,29));
[CONSTANTS, STATES, ALGEBRAIC] = rootfind_4(VOI, CONSTANTS, STATES, ALGEBRAIC);
ALGEBRAIC(:,60) = ( CONSTANTS(:,43).*ALGEBRAIC(:,57).*STATES(:,8))./(CONSTANTS(:,44)+STATES(:,8));
RATES(:,8) = ALGEBRAIC(:,83) - ALGEBRAIC(:,60);
ALGEBRAIC(:,30) = CONSTANTS(:,40) - STATES(:,9);
ALGEBRAIC(:,84) = ( CONSTANTS(:,45).*ALGEBRAIC(:,68).*ALGEBRAIC(:,30))./(CONSTANTS(:,46)+ALGEBRAIC(:,30));
ALGEBRAIC(:,32) = ( CONSTANTS(:,47).*STATES(:,9))./(CONSTANTS(:,48)+STATES(:,9));
RATES(:,9) = ALGEBRAIC(:,84) - ALGEBRAIC(:,32);
ALGEBRAIC(:,34) = CONSTANTS(:,51) - STATES(:,10);
ALGEBRAIC(:,85) = ( CONSTANTS(:,50).*CONSTANTS(:,54).*ALGEBRAIC(:,78).*ALGEBRAIC(:,34))./(CONSTANTS(:,55)+ CONSTANTS(:,50).*ALGEBRAIC(:,34));
ALGEBRAIC(:,36) = ( CONSTANTS(:,50).*CONSTANTS(:,58).*CONSTANTS(:,53).*STATES(:,10))./(CONSTANTS(:,59)+ CONSTANTS(:,50).*STATES(:,10));
RATES(:,10) = ALGEBRAIC(:,85) - ALGEBRAIC(:,36);
ALGEBRAIC(:,38) = CONSTANTS(:,51) - STATES(:,11);
ALGEBRAIC(:,86) = ( CONSTANTS(:,50).*CONSTANTS(:,54).*ALGEBRAIC(:,78).*ALGEBRAIC(:,38))./(CONSTANTS(:,55)+ CONSTANTS(:,50).*ALGEBRAIC(:,38));
ALGEBRAIC(:,40) = ( CONSTANTS(:,50).*CONSTANTS(:,56).*CONSTANTS(:,52).*STATES(:,11))./(CONSTANTS(:,57)+ CONSTANTS(:,50).*STATES(:,11));
RATES(:,11) = ALGEBRAIC(:,86) - ALGEBRAIC(:,40);
ALGEBRAIC(:,88) = CONSTANTS(:,10).*ALGEBRAIC(:,68).*STATES(:,1);
ALGEBRAIC(:,23) = CONSTANTS(:,11).*STATES(:,3);
RATES(:,1) = (ALGEBRAIC(:,12) - ALGEBRAIC(:,80))+(ALGEBRAIC(:,23) - ALGEBRAIC(:,88));
RATES(:,3) = ALGEBRAIC(:,88) - ALGEBRAIC(:,23);
ALGEBRAIC(:,35) = (1.00000./( CONSTANTS(:,125).*CONSTANTS(:,115))).*log(CONSTANTS(:,65)./STATES(:,12));
ALGEBRAIC(:,41) = CONSTANTS(:,68).*power(STATES(:,15), 3.00000).*STATES(:,16).*STATES(:,17).*(STATES(:,18) - ALGEBRAIC(:,35));
ALGEBRAIC(:,93) = exp( CONSTANTS(:,85).*STATES(:,18).*CONSTANTS(:,125)).*power(STATES(:,12), 3.00000).*CONSTANTS(:,67);
ALGEBRAIC(:,94) = exp( (CONSTANTS(:,85) - 1.00000).*STATES(:,18).*CONSTANTS(:,125)).*power(CONSTANTS(:,65), 3.00000).*STATES(:,14);
ALGEBRAIC(:,95) = (CONSTANTS(:,81)./( (power(CONSTANTS(:,82), 3.00000)+power(CONSTANTS(:,65), 3.00000)).*(CONSTANTS(:,83)+CONSTANTS(:,67)).*(1.00000+ CONSTANTS(:,84).*exp( (CONSTANTS(:,85) - 1.00000).*STATES(:,18).*CONSTANTS(:,125))))).*(ALGEBRAIC(:,93) - ALGEBRAIC(:,94));
ALGEBRAIC(:,96) = 1.00000./(1.00000+ 0.124500.*exp( - 0.100000.*STATES(:,18).*CONSTANTS(:,125))+ 0.0365000.*CONSTANTS(:,127).*exp( - STATES(:,18).*CONSTANTS(:,125)));
ALGEBRAIC(:,97) = ( (( CONSTANTS(:,86).*ALGEBRAIC(:,96))./(1.00000+power(CONSTANTS(:,87)./STATES(:,12), 1.50000))).*CONSTANTS(:,66))./(CONSTANTS(:,66)+CONSTANTS(:,88));
ALGEBRAIC(:,100) = CONSTANTS(:,92).*(STATES(:,18) - ALGEBRAIC(:,35));
ALGEBRAIC(:,101) = ALGEBRAIC(:,41)+ALGEBRAIC(:,100)+ 3.00000.*ALGEBRAIC(:,95)+ 3.00000.*ALGEBRAIC(:,97);
RATES(:,12) = ( - 1000.00.*ALGEBRAIC(:,101).*CONSTANTS(:,63))./( CONSTANTS(:,60).*CONSTANTS(:,115).*CONSTANTS(:,114));
ALGEBRAIC(:,47) = ( CONSTANTS(:,78).*STATES(:,18).*CONSTANTS(:,114).*CONSTANTS(:,125).*( STATES(:,13).*exp( STATES(:,18).*CONSTANTS(:,125)) - CONSTANTS(:,66)))./(exp( STATES(:,18).*CONSTANTS(:,125)) - 1.00000);
ALGEBRAIC(:,33) = STATES(:,11)./CONSTANTS(:,51);
ALGEBRAIC(:,49) = 0.500000.*(( 0.400000.*ALGEBRAIC(:,33))./CONSTANTS(:,123)+0.600000);
ALGEBRAIC(:,45) = ( CONSTANTS(:,77).*4.00000.*STATES(:,18).*CONSTANTS(:,114).*CONSTANTS(:,125).*( 0.00100000.*exp( 2.00000.*STATES(:,18).*CONSTANTS(:,125)) - 0.341000.*CONSTANTS(:,67)))./(exp( 2.00000.*STATES(:,18).*CONSTANTS(:,125)) - 1.00000);
ALGEBRAIC(:,52) = ALGEBRAIC(:,45).*CONSTANTS(:,79).*ALGEBRAIC(:,49).*power(STATES(:,19), 4.00000).*STATES(:,21).*STATES(:,22).*STATES(:,23);
ALGEBRAIC(:,54) = (ALGEBRAIC(:,47)./(1.00000+ALGEBRAIC(:,52)./CONSTANTS(:,80))).*CONSTANTS(:,79).*ALGEBRAIC(:,49).*power(STATES(:,19), 4.00000).*STATES(:,21).*STATES(:,22).*STATES(:,23);
ALGEBRAIC(:,37) = (1.00000./( CONSTANTS(:,125).*CONSTANTS(:,116))).*log(CONSTANTS(:,66)./STATES(:,13));
ALGEBRAIC(:,59) = CONSTANTS(:,69).*STATES(:,24).*( 0.886000.*STATES(:,25)+ 0.114000.*STATES(:,26)).*(STATES(:,18) - ALGEBRAIC(:,37));
ALGEBRAIC(:,61) = CONSTANTS(:,70).*STATES(:,27).*STATES(:,28).*(STATES(:,18) - ALGEBRAIC(:,37));
ALGEBRAIC(:,79) = 1.02000./(1.00000+exp( 0.238500.*((STATES(:,18) - ALGEBRAIC(:,37)) - 59.2150)));
ALGEBRAIC(:,87) = ( 0.491240.*exp( 0.0803200.*((STATES(:,18)+5.47600) - ALGEBRAIC(:,37)))+exp( 0.0617500.*((STATES(:,18) - ALGEBRAIC(:,37)) - 594.310)))./(1.00000+exp( - 0.514300.*((STATES(:,18) - ALGEBRAIC(:,37))+4.75300)));
ALGEBRAIC(:,89) = ALGEBRAIC(:,79)./(ALGEBRAIC(:,79)+ALGEBRAIC(:,87));
ALGEBRAIC(:,90) = CONSTANTS(:,71).*power((CONSTANTS(:,66)./5.40000), 1.0 ./ 2).*ALGEBRAIC(:,89).*(STATES(:,18) - ALGEBRAIC(:,37));
ALGEBRAIC(:,91) = 1.00000./(1.00000+exp((7.48800 - STATES(:,18))./5.98000));
ALGEBRAIC(:,92) = CONSTANTS(:,72).*ALGEBRAIC(:,91).*(STATES(:,18) - ALGEBRAIC(:,37));
ALGEBRAIC(:,102) = ((ALGEBRAIC(:,59)+ALGEBRAIC(:,61)+ALGEBRAIC(:,90)+ALGEBRAIC(:,92)) - 2.00000.*ALGEBRAIC(:,97))+ALGEBRAIC(:,54);
RATES(:,13) = ( - 1000.00.*ALGEBRAIC(:,102).*CONSTANTS(:,63))./( CONSTANTS(:,60).*CONSTANTS(:,116).*CONSTANTS(:,114));
ALGEBRAIC(:,98) = ( CONSTANTS(:,89).*STATES(:,14))./(CONSTANTS(:,90)+STATES(:,14));
ALGEBRAIC(:,39) = (1.00000./( CONSTANTS(:,125).*CONSTANTS(:,117))).*log(CONSTANTS(:,67)./STATES(:,14));
ALGEBRAIC(:,99) = CONSTANTS(:,91).*(STATES(:,18) - ALGEBRAIC(:,39));
ALGEBRAIC(:,103) = (ALGEBRAIC(:,52)+ALGEBRAIC(:,99)+ALGEBRAIC(:,98)) - 2.00000.*ALGEBRAIC(:,95);
ALGEBRAIC(:,105) = piecewise({VOI>59.1000&VOI<59.5000, CONSTANTS(:,119) }, CONSTANTS(:,118));
ALGEBRAIC(:,107) = piecewise({ rem(VOI+0.900000, 1.00000)<0.00500000, 10.0000 }, 0.00000);
ALGEBRAIC(:,109) = piecewise({CONSTANTS(:,1)==0.00000, 0.00000 , CONSTANTS(:,1)==1.00000, ALGEBRAIC(:,107) }, (ALGEBRAIC(:,105) - STATES(:,18))./0.0200000);
RATES(:,18) = - 1000.00.*((ALGEBRAIC(:,103)+ALGEBRAIC(:,102)+ALGEBRAIC(:,101)) - ALGEBRAIC(:,109));
RATES(:,30) = piecewise({ - 1000.00.*((ALGEBRAIC(:,103)+ALGEBRAIC(:,102)+ALGEBRAIC(:,101)) - ALGEBRAIC(:,109))>30000.0, 1.00000 - 10000.0.*STATES(:,30) }, 1.00000);
ALGEBRAIC(:,31) = ALGEBRAIC(:,29)./CONSTANTS(:,38);
ALGEBRAIC(:,112) = ( CONSTANTS(:,96).*(1.00000+ 2.00000.*ALGEBRAIC(:,31)))./(1.00000+ 2.00000.*CONSTANTS(:,121));
ALGEBRAIC(:,113) = ( CONSTANTS(:,95).*power(STATES(:,14), 2.00000))./(power(ALGEBRAIC(:,112), 2.00000)+power(STATES(:,14), 2.00000));
ALGEBRAIC(:,114) = ( CONSTANTS(:,95).*STATES(:,31))./CONSTANTS(:,97);
ALGEBRAIC(:,115) = (STATES(:,31) - STATES(:,29))./CONSTANTS(:,106);
RATES(:,31) = (ALGEBRAIC(:,113) - ALGEBRAIC(:,114)) - ( ALGEBRAIC(:,115).*CONSTANTS(:,62))./CONSTANTS(:,61);
ALGEBRAIC(:,104) = STATES(:,30)+0.00200000;
ALGEBRAIC(:,106) = 1.00000 - exp( - ALGEBRAIC(:,104)./CONSTANTS(:,98));
ALGEBRAIC(:,108) = exp( - ALGEBRAIC(:,104)./CONSTANTS(:,99));
ALGEBRAIC(:,110) = CONSTANTS(:,100)./(1.00000+exp((ALGEBRAIC(:,103)+5.00000)./0.900000));
ALGEBRAIC(:,111) = ALGEBRAIC(:,110).*ALGEBRAIC(:,106).*ALGEBRAIC(:,108).*(STATES(:,29) - STATES(:,14));
ALGEBRAIC(:,117) = 1.00000./(1.00000+( CONSTANTS(:,104).*CONSTANTS(:,105))./power(CONSTANTS(:,105)+STATES(:,29), 2.00000));
RATES(:,29) = ALGEBRAIC(:,117).*(ALGEBRAIC(:,115) - ALGEBRAIC(:,111));
ALGEBRAIC(:,116) = ( CONSTANTS(:,107).*CONSTANTS(:,110))./power(CONSTANTS(:,110)+STATES(:,14), 2.00000);
ALGEBRAIC(:,118) = ( CONSTANTS(:,108).*CONSTANTS(:,111))./power(CONSTANTS(:,111)+STATES(:,14), 2.00000);
ALGEBRAIC(:,120) = ( CONSTANTS(:,109).*CONSTANTS(:,112))./power(CONSTANTS(:,112)+STATES(:,14), 2.00000);
ALGEBRAIC(:,121) = 1.00000./(1.00000+ALGEBRAIC(:,118)+ALGEBRAIC(:,116)+ALGEBRAIC(:,116)+ALGEBRAIC(:,120));
RATES(:,14) = - ALGEBRAIC(:,121).*((( 1000.00.*ALGEBRAIC(:,103).*CONSTANTS(:,63))./( CONSTANTS(:,60).*CONSTANTS(:,117).*CONSTANTS(:,114))+( (ALGEBRAIC(:,113) - ALGEBRAIC(:,114)).*CONSTANTS(:,61))./CONSTANTS(:,60)) - ( ALGEBRAIC(:,111).*CONSTANTS(:,62))./CONSTANTS(:,60));
RATES = RATES';
end
% Calculate algebraic variables
function ALGEBRAIC = computeAlgebraic(ALGEBRAIC, CONSTANTS, STATES, VOI)
statesSize = size(STATES);
statesColumnCount = statesSize(2);
if ( statesColumnCount == 1)
STATES = STATES';
utilOnes = 1;
else
statesRowCount = statesSize(1);
utilOnes = ones(statesRowCount, 1);
end
ALGEBRAIC(:,11) = 1.00000./(1.00000+exp((STATES(:,18)+87.5000)./10.3000));
ALGEBRAIC(:,3) = ( 0.320000.*(STATES(:,18)+47.1300))./(1.00000 - exp( - 0.100000.*(STATES(:,18)+47.1300)));
ALGEBRAIC(:,13) = 0.0800000.*exp( - STATES(:,18)./11.0000);
ALGEBRAIC(:,4) = piecewise({STATES(:,18)>= - 40.0000, 0.00000 }, 0.135000.*exp((80.0000+STATES(:,18))./ - 6.80000));
ALGEBRAIC(:,14) = piecewise({STATES(:,18)>= - 40.0000, 1.00000./( 0.130000.*(1.00000+exp( - (STATES(:,18)+10.6600)./11.1000))) }, 3.56000.*exp( 0.0790000.*STATES(:,18))+ 310000..*exp( 0.350000.*STATES(:,18)));
ALGEBRAIC(:,5) = piecewise({STATES(:,18)>= - 40.0000, 0.00000 }, ( ( - 127140..*exp( 0.244400.*STATES(:,18)) - 3.47400e-05.*exp( - 0.0439100.*STATES(:,18))).*(STATES(:,18)+37.7800))./(1.00000+exp( 0.311000.*(STATES(:,18)+79.2300))));
ALGEBRAIC(:,15) = piecewise({STATES(:,18)>= - 40.0000, ( 0.300000.*exp( - 2.57500e-07.*STATES(:,18)))./(1.00000+exp( - 0.100000.*(STATES(:,18)+32.0000))) }, ( 0.121200.*exp( - 0.0105200.*STATES(:,18)))./(1.00000+exp( - 0.137800.*(STATES(:,18)+40.1400))));
ALGEBRAIC(:,6) = 400.000.*exp((STATES(:,18)+2.00000)./10.0000);
ALGEBRAIC(:,16) = 50.0000.*exp(( - 1.00000.*(STATES(:,18)+2.00000))./13.0000);
ALGEBRAIC(:,2) = STATES(:,10)./CONSTANTS(:,51);
ALGEBRAIC(:,18) = CONSTANTS(:,73).*(( 0.375000.*ALGEBRAIC(:,2))./CONSTANTS(:,122)+0.625000);
ALGEBRAIC(:,7) = 1.00000./(1.00000+exp((STATES(:,18)+55.0000)./7.50000))+0.100000./(1.00000+exp(( - STATES(:,18)+21.0000)./6.00000));
ALGEBRAIC(:,17) = 0.0200000+0.300000./(1.00000+exp((STATES(:,18)+30.0000)./9.50000));
ALGEBRAIC(:,8) = 1.00000./(1.00000+exp((STATES(:,18)+10.6000)./ - 11.4200));
ALGEBRAIC(:,19) = 1.00000./( 45.1600.*exp( 0.0357700.*(STATES(:,18)+50.0000))+ 98.9000.*exp( - 0.100000.*(STATES(:,18)+38.0000)));
ALGEBRAIC(:,9) = 1.00000./(1.00000+exp((STATES(:,18)+43.5000)./6.88410));
ALGEBRAIC(:,20) = 0.350000.*exp( - 1.00000.*power((STATES(:,18)+70.0000)./15.0000, 2.00000))+0.0350000;
ALGEBRAIC(:,21) = 3.70000.*exp( - 1.00000.*power((STATES(:,18)+70.0000)./30.0000, 2.00000))+0.0350000;
ALGEBRAIC(:,10) = 1.00000./(1.00000+exp((STATES(:,18)+11.5000)./ - 11.8200));
ALGEBRAIC(:,22) = 10.0000./( 45.1600.*exp( 0.0357700.*(STATES(:,18)+50.0000))+ 98.9000.*exp( - 0.100000.*(STATES(:,18)+38.0000)));
ALGEBRAIC(:,25) = CONSTANTS(:,13).*STATES(:,4);
ALGEBRAIC(:,27) = CONSTANTS(:,14).*STATES(:,5).*STATES(:,6);
ALGEBRAIC(:,24) = CONSTANTS(:,75).*STATES(:,14);
ALGEBRAIC(:,26) = ALGEBRAIC(:,24).*(power(1.00000 - STATES(:,19), 4.00000)+ 2.00000.*STATES(:,19).*power(1.00000 - STATES(:,19), 3.00000)+ 4.00000.*power(STATES(:,19), 2.00000).*power(1.00000 - STATES(:,19), 2.00000)+ 8.00000.*power(STATES(:,19), 3.00000).*(1.00000 - STATES(:,19))+ 16.0000.*power(STATES(:,19), 4.00000).*(1.00000 - ALGEBRAIC(:,18)./CONSTANTS(:,74)));
ALGEBRAIC(:,28) = CONSTANTS(:,76).*(power(1.00000 - STATES(:,20), 4.00000)+ 0.500000.*STATES(:,20).*power(1.00000 - STATES(:,20), 3.00000)+ 0.250000.*power(STATES(:,20), 2.00000).*power(1.00000 - STATES(:,20), 2.00000)+ 0.125000.*power(STATES(:,20), 3.00000).*(1.00000 - STATES(:,20))+ 0.0625000.*power(STATES(:,20), 4.00000));
ALGEBRAIC(:,80) = CONSTANTS(:,8).*(ALGEBRAIC(:,62)+ALGEBRAIC(:,63));
ALGEBRAIC(:,12) = CONSTANTS(:,9).*STATES(:,2);
ALGEBRAIC(:,81) = CONSTANTS(:,12).*(ALGEBRAIC(:,64)+ALGEBRAIC(:,63));
ALGEBRAIC(:,48) = ( CONSTANTS(:,21).*ALGEBRAIC(:,44).*CONSTANTS(:,16))./(CONSTANTS(:,25)+CONSTANTS(:,16));
ALGEBRAIC(:,46) = ( CONSTANTS(:,20).*ALGEBRAIC(:,43).*CONSTANTS(:,16))./(CONSTANTS(:,24)+CONSTANTS(:,16));
ALGEBRAIC(:,53) = ( CONSTANTS(:,22).*ALGEBRAIC(:,51).*CONSTANTS(:,16))./(CONSTANTS(:,26)+CONSTANTS(:,16));
ALGEBRAIC(:,82) = ( CONSTANTS(:,23).*CONSTANTS(:,128).*ALGEBRAIC(:,69))./(CONSTANTS(:,27)+ALGEBRAIC(:,69));
ALGEBRAIC(:,29) = CONSTANTS(:,38) - STATES(:,8);
ALGEBRAIC(:,83) = ( CONSTANTS(:,41).*ALGEBRAIC(:,68).*ALGEBRAIC(:,29))./(CONSTANTS(:,42)+ALGEBRAIC(:,29));
ALGEBRAIC(:,60) = ( CONSTANTS(:,43).*ALGEBRAIC(:,57).*STATES(:,8))./(CONSTANTS(:,44)+STATES(:,8));
ALGEBRAIC(:,30) = CONSTANTS(:,40) - STATES(:,9);
ALGEBRAIC(:,84) = ( CONSTANTS(:,45).*ALGEBRAIC(:,68).*ALGEBRAIC(:,30))./(CONSTANTS(:,46)+ALGEBRAIC(:,30));
ALGEBRAIC(:,32) = ( CONSTANTS(:,47).*STATES(:,9))./(CONSTANTS(:,48)+STATES(:,9));
ALGEBRAIC(:,34) = CONSTANTS(:,51) - STATES(:,10);
ALGEBRAIC(:,85) = ( CONSTANTS(:,50).*CONSTANTS(:,54).*ALGEBRAIC(:,78).*ALGEBRAIC(:,34))./(CONSTANTS(:,55)+ CONSTANTS(:,50).*ALGEBRAIC(:,34));
ALGEBRAIC(:,36) = ( CONSTANTS(:,50).*CONSTANTS(:,58).*CONSTANTS(:,53).*STATES(:,10))./(CONSTANTS(:,59)+ CONSTANTS(:,50).*STATES(:,10));
ALGEBRAIC(:,38) = CONSTANTS(:,51) - STATES(:,11);
ALGEBRAIC(:,86) = ( CONSTANTS(:,50).*CONSTANTS(:,54).*ALGEBRAIC(:,78).*ALGEBRAIC(:,38))./(CONSTANTS(:,55)+ CONSTANTS(:,50).*ALGEBRAIC(:,38));
ALGEBRAIC(:,40) = ( CONSTANTS(:,50).*CONSTANTS(:,56).*CONSTANTS(:,52).*STATES(:,11))./(CONSTANTS(:,57)+ CONSTANTS(:,50).*STATES(:,11));
ALGEBRAIC(:,88) = CONSTANTS(:,10).*ALGEBRAIC(:,68).*STATES(:,1);
ALGEBRAIC(:,23) = CONSTANTS(:,11).*STATES(:,3);
ALGEBRAIC(:,35) = (1.00000./( CONSTANTS(:,125).*CONSTANTS(:,115))).*log(CONSTANTS(:,65)./STATES(:,12));
ALGEBRAIC(:,41) = CONSTANTS(:,68).*power(STATES(:,15), 3.00000).*STATES(:,16).*STATES(:,17).*(STATES(:,18) - ALGEBRAIC(:,35));
ALGEBRAIC(:,93) = exp( CONSTANTS(:,85).*STATES(:,18).*CONSTANTS(:,125)).*power(STATES(:,12), 3.00000).*CONSTANTS(:,67);
ALGEBRAIC(:,94) = exp( (CONSTANTS(:,85) - 1.00000).*STATES(:,18).*CONSTANTS(:,125)).*power(CONSTANTS(:,65), 3.00000).*STATES(:,14);
ALGEBRAIC(:,95) = (CONSTANTS(:,81)./( (power(CONSTANTS(:,82), 3.00000)+power(CONSTANTS(:,65), 3.00000)).*(CONSTANTS(:,83)+CONSTANTS(:,67)).*(1.00000+ CONSTANTS(:,84).*exp( (CONSTANTS(:,85) - 1.00000).*STATES(:,18).*CONSTANTS(:,125))))).*(ALGEBRAIC(:,93) - ALGEBRAIC(:,94));
ALGEBRAIC(:,96) = 1.00000./(1.00000+ 0.124500.*exp( - 0.100000.*STATES(:,18).*CONSTANTS(:,125))+ 0.0365000.*CONSTANTS(:,127).*exp( - STATES(:,18).*CONSTANTS(:,125)));
ALGEBRAIC(:,97) = ( (( CONSTANTS(:,86).*ALGEBRAIC(:,96))./(1.00000+power(CONSTANTS(:,87)./STATES(:,12), 1.50000))).*CONSTANTS(:,66))./(CONSTANTS(:,66)+CONSTANTS(:,88));
ALGEBRAIC(:,100) = CONSTANTS(:,92).*(STATES(:,18) - ALGEBRAIC(:,35));
ALGEBRAIC(:,101) = ALGEBRAIC(:,41)+ALGEBRAIC(:,100)+ 3.00000.*ALGEBRAIC(:,95)+ 3.00000.*ALGEBRAIC(:,97);
ALGEBRAIC(:,47) = ( CONSTANTS(:,78).*STATES(:,18).*CONSTANTS(:,114).*CONSTANTS(:,125).*( STATES(:,13).*exp( STATES(:,18).*CONSTANTS(:,125)) - CONSTANTS(:,66)))./(exp( STATES(:,18).*CONSTANTS(:,125)) - 1.00000);
ALGEBRAIC(:,33) = STATES(:,11)./CONSTANTS(:,51);
ALGEBRAIC(:,49) = 0.500000.*(( 0.400000.*ALGEBRAIC(:,33))./CONSTANTS(:,123)+0.600000);
ALGEBRAIC(:,45) = ( CONSTANTS(:,77).*4.00000.*STATES(:,18).*CONSTANTS(:,114).*CONSTANTS(:,125).*( 0.00100000.*exp( 2.00000.*STATES(:,18).*CONSTANTS(:,125)) - 0.341000.*CONSTANTS(:,67)))./(exp( 2.00000.*STATES(:,18).*CONSTANTS(:,125)) - 1.00000);
ALGEBRAIC(:,52) = ALGEBRAIC(:,45).*CONSTANTS(:,79).*ALGEBRAIC(:,49).*power(STATES(:,19), 4.00000).*STATES(:,21).*STATES(:,22).*STATES(:,23);
ALGEBRAIC(:,54) = (ALGEBRAIC(:,47)./(1.00000+ALGEBRAIC(:,52)./CONSTANTS(:,80))).*CONSTANTS(:,79).*ALGEBRAIC(:,49).*power(STATES(:,19), 4.00000).*STATES(:,21).*STATES(:,22).*STATES(:,23);
ALGEBRAIC(:,37) = (1.00000./( CONSTANTS(:,125).*CONSTANTS(:,116))).*log(CONSTANTS(:,66)./STATES(:,13));
ALGEBRAIC(:,59) = CONSTANTS(:,69).*STATES(:,24).*( 0.886000.*STATES(:,25)+ 0.114000.*STATES(:,26)).*(STATES(:,18) - ALGEBRAIC(:,37));
ALGEBRAIC(:,61) = CONSTANTS(:,70).*STATES(:,27).*STATES(:,28).*(STATES(:,18) - ALGEBRAIC(:,37));
ALGEBRAIC(:,79) = 1.02000./(1.00000+exp( 0.238500.*((STATES(:,18) - ALGEBRAIC(:,37)) - 59.2150)));
ALGEBRAIC(:,87) = ( 0.491240.*exp( 0.0803200.*((STATES(:,18)+5.47600) - ALGEBRAIC(:,37)))+exp( 0.0617500.*((STATES(:,18) - ALGEBRAIC(:,37)) - 594.310)))./(1.00000+exp( - 0.514300.*((STATES(:,18) - ALGEBRAIC(:,37))+4.75300)));
ALGEBRAIC(:,89) = ALGEBRAIC(:,79)./(ALGEBRAIC(:,79)+ALGEBRAIC(:,87));
ALGEBRAIC(:,90) = CONSTANTS(:,71).*power((CONSTANTS(:,66)./5.40000), 1.0 ./ 2).*ALGEBRAIC(:,89).*(STATES(:,18) - ALGEBRAIC(:,37));
ALGEBRAIC(:,91) = 1.00000./(1.00000+exp((7.48800 - STATES(:,18))./5.98000));
ALGEBRAIC(:,92) = CONSTANTS(:,72).*ALGEBRAIC(:,91).*(STATES(:,18) - ALGEBRAIC(:,37));
ALGEBRAIC(:,102) = ((ALGEBRAIC(:,59)+ALGEBRAIC(:,61)+ALGEBRAIC(:,90)+ALGEBRAIC(:,92)) - 2.00000.*ALGEBRAIC(:,97))+ALGEBRAIC(:,54);
ALGEBRAIC(:,98) = ( CONSTANTS(:,89).*STATES(:,14))./(CONSTANTS(:,90)+STATES(:,14));
ALGEBRAIC(:,39) = (1.00000./( CONSTANTS(:,125).*CONSTANTS(:,117))).*log(CONSTANTS(:,67)./STATES(:,14));
ALGEBRAIC(:,99) = CONSTANTS(:,91).*(STATES(:,18) - ALGEBRAIC(:,39));
ALGEBRAIC(:,103) = (ALGEBRAIC(:,52)+ALGEBRAIC(:,99)+ALGEBRAIC(:,98)) - 2.00000.*ALGEBRAIC(:,95);
ALGEBRAIC(:,105) = piecewise({VOI>59.1000&VOI<59.5000, CONSTANTS(:,119) }, CONSTANTS(:,118));
ALGEBRAIC(:,107) = piecewise({ rem(VOI+0.900000, 1.00000)<0.00500000, 10.0000 }, 0.00000);
ALGEBRAIC(:,109) = piecewise({CONSTANTS(:,1)==0.00000, 0.00000 , CONSTANTS(:,1)==1.00000, ALGEBRAIC(:,107) }, (ALGEBRAIC(:,105) - STATES(:,18))./0.0200000);
ALGEBRAIC(:,31) = ALGEBRAIC(:,29)./CONSTANTS(:,38);
ALGEBRAIC(:,112) = ( CONSTANTS(:,96).*(1.00000+ 2.00000.*ALGEBRAIC(:,31)))./(1.00000+ 2.00000.*CONSTANTS(:,121));
ALGEBRAIC(:,113) = ( CONSTANTS(:,95).*power(STATES(:,14), 2.00000))./(power(ALGEBRAIC(:,112), 2.00000)+power(STATES(:,14), 2.00000));
ALGEBRAIC(:,114) = ( CONSTANTS(:,95).*STATES(:,31))./CONSTANTS(:,97);
ALGEBRAIC(:,115) = (STATES(:,31) - STATES(:,29))./CONSTANTS(:,106);
ALGEBRAIC(:,104) = STATES(:,30)+0.00200000;
ALGEBRAIC(:,106) = 1.00000 - exp( - ALGEBRAIC(:,104)./CONSTANTS(:,98));
ALGEBRAIC(:,108) = exp( - ALGEBRAIC(:,104)./CONSTANTS(:,99));
ALGEBRAIC(:,110) = CONSTANTS(:,100)./(1.00000+exp((ALGEBRAIC(:,103)+5.00000)./0.900000));
ALGEBRAIC(:,111) = ALGEBRAIC(:,110).*ALGEBRAIC(:,106).*ALGEBRAIC(:,108).*(STATES(:,29) - STATES(:,14));
ALGEBRAIC(:,117) = 1.00000./(1.00000+( CONSTANTS(:,104).*CONSTANTS(:,105))./power(CONSTANTS(:,105)+STATES(:,29), 2.00000));
ALGEBRAIC(:,116) = ( CONSTANTS(:,107).*CONSTANTS(:,110))./power(CONSTANTS(:,110)+STATES(:,14), 2.00000);
ALGEBRAIC(:,118) = ( CONSTANTS(:,108).*CONSTANTS(:,111))./power(CONSTANTS(:,111)+STATES(:,14), 2.00000);
ALGEBRAIC(:,120) = ( CONSTANTS(:,109).*CONSTANTS(:,112))./power(CONSTANTS(:,112)+STATES(:,14), 2.00000);
ALGEBRAIC(:,121) = 1.00000./(1.00000+ALGEBRAIC(:,118)+ALGEBRAIC(:,116)+ALGEBRAIC(:,116)+ALGEBRAIC(:,120));
ALGEBRAIC(:,1) = STATES(:,8)./CONSTANTS(:,38);
ALGEBRAIC(:,58) = ALGEBRAIC(:,52)+ALGEBRAIC(:,54);
ALGEBRAIC(:,119) = 1000.00.*(( (STATES(:,29)+STATES(:,29)./ALGEBRAIC(:,117)).*CONSTANTS(:,62))./CONSTANTS(:,60)+( STATES(:,31).*CONSTANTS(:,61))./CONSTANTS(:,60));
end
% Functions required for solving differential algebraic equation
function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_0(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN)
ALGEBRAIC = ALGEBRAIC_IN;
CONSTANTS = CONSTANTS_IN;
STATES = STATES_IN;
global initialGuess_0;
if (length(initialGuess_0) ~= 3), initialGuess_0 = [0.0389,0,0.1];, end
options = optimset('Display', 'off', 'TolX', 1E-6);
if length(VOI) == 1
residualfn = @(algebraicCandidate)residualSN_0(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES);
soln = fsolve(residualfn, initialGuess_0, options);
initialGuess_0 = soln;
CONSTANTS(:,128) = soln(1);
CONSTANTS(:,129) = soln(2);
CONSTANTS(:,130) = soln(3);
else
SET_CONSTANTS(:,128) = logical(1);
SET_CONSTANTS(:,129) = logical(1);
SET_CONSTANTS(:,130) = logical(1);
for i=1:length(VOI)
residualfn = @(algebraicCandidate)residualSN_0(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:));
soln = fsolve(residualfn, initialGuess_0, options);
initialGuess_0 = soln;
TEMP_CONSTANTS(:,128) = soln(1);
TEMP_CONSTANTS(:,129) = soln(2);
TEMP_CONSTANTS(:,130) = soln(3);
ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC);
end
end
end
function resid = residualSN_0(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES)
CONSTANTS(:,128) = algebraicCandidate(1);
CONSTANTS(:,129) = algebraicCandidate(2);
CONSTANTS(:,130) = algebraicCandidate(3);
resid(1) = CONSTANTS(:,130) - ( CONSTANTS(:,128).*CONSTANTS(:,129))./CONSTANTS(:,30);
resid(2) = CONSTANTS(:,128) - (CONSTANTS(:,17) - CONSTANTS(:,130));
resid(3) = CONSTANTS(:,129) - (CONSTANTS(:,18) - CONSTANTS(:,130));
end
% Functions required for solving differential algebraic equation
function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_1(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN)
ALGEBRAIC = ALGEBRAIC_IN;
CONSTANTS = CONSTANTS_IN;
STATES = STATES_IN;
global initialGuess_1;
if (length(initialGuess_1) ~= 17), initialGuess_1 = [0.000191562,0.0117971,6.39347e-06,0.988,5.5258e-05,3.8182,0.0588,0.227,0.000534718,0.00290213,0.21597,5.83396e-05,0.0306208,0.116856,0.00234908,3.91219,0.0083];, end
options = optimset('Display', 'off', 'TolX', 1E-6);
if length(VOI) == 1
residualfn = @(algebraicCandidate)residualSN_1(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES);
soln = fsolve(residualfn, initialGuess_1, options);
initialGuess_1 = soln;
ALGEBRAIC(:,62) = soln(1);
ALGEBRAIC(:,63) = soln(2);
ALGEBRAIC(:,64) = soln(3);
ALGEBRAIC(:,65) = soln(4);
ALGEBRAIC(:,66) = soln(5);
ALGEBRAIC(:,67) = soln(6);
ALGEBRAIC(:,68) = soln(7);
ALGEBRAIC(:,69) = soln(8);
ALGEBRAIC(:,70) = soln(9);
ALGEBRAIC(:,71) = soln(10);
ALGEBRAIC(:,72) = soln(11);
ALGEBRAIC(:,73) = soln(12);
ALGEBRAIC(:,74) = soln(13);
ALGEBRAIC(:,75) = soln(14);
ALGEBRAIC(:,76) = soln(15);
ALGEBRAIC(:,77) = soln(16);
ALGEBRAIC(:,78) = soln(17);
else
SET_ALGEBRAIC(:,62) = logical(1);
SET_ALGEBRAIC(:,63) = logical(1);
SET_ALGEBRAIC(:,64) = logical(1);
SET_ALGEBRAIC(:,65) = logical(1);
SET_ALGEBRAIC(:,66) = logical(1);
SET_ALGEBRAIC(:,67) = logical(1);
SET_ALGEBRAIC(:,68) = logical(1);
SET_ALGEBRAIC(:,69) = logical(1);
SET_ALGEBRAIC(:,70) = logical(1);
SET_ALGEBRAIC(:,71) = logical(1);
SET_ALGEBRAIC(:,72) = logical(1);
SET_ALGEBRAIC(:,73) = logical(1);
SET_ALGEBRAIC(:,74) = logical(1);
SET_ALGEBRAIC(:,75) = logical(1);
SET_ALGEBRAIC(:,76) = logical(1);
SET_ALGEBRAIC(:,77) = logical(1);
SET_ALGEBRAIC(:,78) = logical(1);
for i=1:length(VOI)
residualfn = @(algebraicCandidate)residualSN_1(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:));
soln = fsolve(residualfn, initialGuess_1, options);
initialGuess_1 = soln;
TEMP_ALGEBRAIC(:,62) = soln(1);
TEMP_ALGEBRAIC(:,63) = soln(2);
TEMP_ALGEBRAIC(:,64) = soln(3);
TEMP_ALGEBRAIC(:,65) = soln(4);
TEMP_ALGEBRAIC(:,66) = soln(5);
TEMP_ALGEBRAIC(:,67) = soln(6);
TEMP_ALGEBRAIC(:,68) = soln(7);
TEMP_ALGEBRAIC(:,69) = soln(8);
TEMP_ALGEBRAIC(:,70) = soln(9);
TEMP_ALGEBRAIC(:,71) = soln(10);
TEMP_ALGEBRAIC(:,72) = soln(11);
TEMP_ALGEBRAIC(:,73) = soln(12);
TEMP_ALGEBRAIC(:,74) = soln(13);
TEMP_ALGEBRAIC(:,75) = soln(14);
TEMP_ALGEBRAIC(:,76) = soln(15);
TEMP_ALGEBRAIC(:,77) = soln(16);
TEMP_ALGEBRAIC(:,78) = soln(17);
ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC);
end
end
end
function resid = residualSN_1(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES)
ALGEBRAIC(:,62) = algebraicCandidate(1);
ALGEBRAIC(:,63) = algebraicCandidate(2);
ALGEBRAIC(:,64) = algebraicCandidate(3);
ALGEBRAIC(:,65) = algebraicCandidate(4);
ALGEBRAIC(:,66) = algebraicCandidate(5);
ALGEBRAIC(:,67) = algebraicCandidate(6);
ALGEBRAIC(:,68) = algebraicCandidate(7);
ALGEBRAIC(:,69) = algebraicCandidate(8);
ALGEBRAIC(:,70) = algebraicCandidate(9);
ALGEBRAIC(:,71) = algebraicCandidate(10);
ALGEBRAIC(:,72) = algebraicCandidate(11);
ALGEBRAIC(:,73) = algebraicCandidate(12);
ALGEBRAIC(:,74) = algebraicCandidate(13);
ALGEBRAIC(:,75) = algebraicCandidate(14);
ALGEBRAIC(:,76) = algebraicCandidate(15);
ALGEBRAIC(:,77) = algebraicCandidate(16);
ALGEBRAIC(:,78) = algebraicCandidate(17);
resid(1) = ALGEBRAIC(:,62) - ( ALGEBRAIC(:,65).*ALGEBRAIC(:,66))./CONSTANTS(:,5);
resid(2) = ALGEBRAIC(:,63) - ( ALGEBRAIC(:,62).*ALGEBRAIC(:,67))./CONSTANTS(:,6);
resid(3) = ALGEBRAIC(:,64) - ( ALGEBRAIC(:,66).*ALGEBRAIC(:,67))./CONSTANTS(:,7);
resid(4) = ALGEBRAIC(:,65) - ((CONSTANTS(:,2) - ALGEBRAIC(:,62)) - ALGEBRAIC(:,63));
resid(5) = ALGEBRAIC(:,66) - (((STATES(:,1) - ALGEBRAIC(:,62)) - ALGEBRAIC(:,63)) - ALGEBRAIC(:,64));
resid(6) = ALGEBRAIC(:,67) - ((CONSTANTS(:,4) - ALGEBRAIC(:,63)) - ALGEBRAIC(:,64));
resid(7) = ALGEBRAIC(:,70) - ( CONSTANTS(:,33).*CONSTANTS(:,37))./(CONSTANTS(:,37)+ALGEBRAIC(:,68)+ALGEBRAIC(:,78));
resid(8) = ALGEBRAIC(:,71) - (ALGEBRAIC(:,68)./CONSTANTS(:,36)).*ALGEBRAIC(:,68).*(1.00000+ALGEBRAIC(:,70)./CONSTANTS(:,37));
resid(9) = ALGEBRAIC(:,72) - ALGEBRAIC(:,68).*(1.00000+ALGEBRAIC(:,70)./CONSTANTS(:,37));
resid(10) = ALGEBRAIC(:,73) - (ALGEBRAIC(:,78)./CONSTANTS(:,36)).*ALGEBRAIC(:,78).*(1.00000+ALGEBRAIC(:,70)./CONSTANTS(:,37));
resid(11) = ALGEBRAIC(:,74) - ALGEBRAIC(:,78).*(1.00000+ALGEBRAIC(:,70)./CONSTANTS(:,37));
resid(12) = ALGEBRAIC(:,75) - (CONSTANTS(:,34)./ALGEBRAIC(:,69)).*ALGEBRAIC(:,71);
resid(13) = ALGEBRAIC(:,76) - (CONSTANTS(:,34)./ALGEBRAIC(:,69)).*ALGEBRAIC(:,73);
resid(14) = ALGEBRAIC(:,77) - (( CONSTANTS(:,34).*CONSTANTS(:,35))./CONSTANTS(:,36)+( CONSTANTS(:,34).*ALGEBRAIC(:,69))./CONSTANTS(:,36)+power(ALGEBRAIC(:,69), 2.00000)./CONSTANTS(:,36));
resid(15) = ALGEBRAIC(:,69) - ((STATES(:,7) - (ALGEBRAIC(:,75)+ 2.00000.*ALGEBRAIC(:,71)+ 2.00000.*ALGEBRAIC(:,72))) - (ALGEBRAIC(:,76)+ 2.00000.*ALGEBRAIC(:,73)+ 2.00000.*ALGEBRAIC(:,74)));
resid(16) = 0.00000 - ( 2.00000.*CONSTANTS(:,31).*power(ALGEBRAIC(:,69), 2.00000) - ALGEBRAIC(:,68).*(1.00000+ALGEBRAIC(:,70)./CONSTANTS(:,37)).*( ALGEBRAIC(:,77).*ALGEBRAIC(:,68)+power(ALGEBRAIC(:,69), 2.00000)));
resid(17) = 0.00000 - ( 2.00000.*CONSTANTS(:,32).*power(ALGEBRAIC(:,69), 2.00000) - ALGEBRAIC(:,78).*(1.00000+ALGEBRAIC(:,70)./CONSTANTS(:,37)).*( ALGEBRAIC(:,77).*ALGEBRAIC(:,78)+power(ALGEBRAIC(:,69), 2.00000)));
end
% Functions required for solving differential algebraic equation
function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_2(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN)
ALGEBRAIC = ALGEBRAIC_IN;
CONSTANTS = CONSTANTS_IN;
STATES = STATES_IN;
global initialGuess_2;
if (length(initialGuess_2) ~= 3), initialGuess_2 = [0.0224,0.0471,0.00263705];, end
options = optimset('Display', 'off', 'TolX', 1E-6);
if length(VOI) == 1
residualfn = @(algebraicCandidate)residualSN_2(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES);
soln = fsolve(residualfn, initialGuess_2, options);
initialGuess_2 = soln;
ALGEBRAIC(:,42) = soln(1);
ALGEBRAIC(:,43) = soln(2);
ALGEBRAIC(:,44) = soln(3);
else
SET_ALGEBRAIC(:,42) = logical(1);
SET_ALGEBRAIC(:,43) = logical(1);
SET_ALGEBRAIC(:,44) = logical(1);
for i=1:length(VOI)
residualfn = @(algebraicCandidate)residualSN_2(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:));
soln = fsolve(residualfn, initialGuess_2, options);
initialGuess_2 = soln;
TEMP_ALGEBRAIC(:,42) = soln(1);
TEMP_ALGEBRAIC(:,43) = soln(2);
TEMP_ALGEBRAIC(:,44) = soln(3);
ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC);
end
end
end
function resid = residualSN_2(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES)
ALGEBRAIC(:,42) = algebraicCandidate(1);
ALGEBRAIC(:,43) = algebraicCandidate(2);
ALGEBRAIC(:,44) = algebraicCandidate(3);
resid(1) = ALGEBRAIC(:,44) - ( ALGEBRAIC(:,42).*ALGEBRAIC(:,43))./CONSTANTS(:,28);
resid(2) = ALGEBRAIC(:,42) - (STATES(:,4) - ALGEBRAIC(:,44));
resid(3) = ALGEBRAIC(:,43) - (CONSTANTS(:,15) - ALGEBRAIC(:,44));
end
% Functions required for solving differential algebraic equation
function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_3(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN)
ALGEBRAIC = ALGEBRAIC_IN;
CONSTANTS = CONSTANTS_IN;
STATES = STATES_IN;
global initialGuess_3;
if (length(initialGuess_3) ~= 2), initialGuess_3 = [0,0];, end
options = optimset('Display', 'off', 'TolX', 1E-6);
if length(VOI) == 1
residualfn = @(algebraicCandidate)residualSN_3(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES);
soln = fsolve(residualfn, initialGuess_3, options);
initialGuess_3 = soln;
ALGEBRAIC(:,50) = soln(1);
ALGEBRAIC(:,51) = soln(2);
else
SET_ALGEBRAIC(:,50) = logical(1);
SET_ALGEBRAIC(:,51) = logical(1);
for i=1:length(VOI)
residualfn = @(algebraicCandidate)residualSN_3(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:));
soln = fsolve(residualfn, initialGuess_3, options);
initialGuess_3 = soln;
TEMP_ALGEBRAIC(:,50) = soln(1);
TEMP_ALGEBRAIC(:,51) = soln(2);
ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC);
end
end
end
function resid = residualSN_3(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES)
ALGEBRAIC(:,50) = algebraicCandidate(1);
ALGEBRAIC(:,51) = algebraicCandidate(2);
resid(1) = ALGEBRAIC(:,51) - ( ALGEBRAIC(:,50).*ALGEBRAIC(:,43))./CONSTANTS(:,29);
resid(2) = ALGEBRAIC(:,50) - (CONSTANTS(:,19) - ALGEBRAIC(:,51));
end
% Functions required for solving differential algebraic equation
function [CONSTANTS, STATES, ALGEBRAIC] = rootfind_4(VOI, CONSTANTS_IN, STATES_IN, ALGEBRAIC_IN)
ALGEBRAIC = ALGEBRAIC_IN;
CONSTANTS = CONSTANTS_IN;
STATES = STATES_IN;
global initialGuess_4;
if (length(initialGuess_4) ~= 3), initialGuess_4 = [0.0525373,6.2734e-05,0.8375];, end
options = optimset('Display', 'off', 'TolX', 1E-6);
if length(VOI) == 1
residualfn = @(algebraicCandidate)residualSN_4(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES);
soln = fsolve(residualfn, initialGuess_4, options);
initialGuess_4 = soln;
ALGEBRAIC(:,55) = soln(1);
ALGEBRAIC(:,56) = soln(2);
ALGEBRAIC(:,57) = soln(3);
else
SET_ALGEBRAIC(:,55) = logical(1);
SET_ALGEBRAIC(:,56) = logical(1);
SET_ALGEBRAIC(:,57) = logical(1);
for i=1:length(VOI)
residualfn = @(algebraicCandidate)residualSN_4(algebraicCandidate, ALGEBRAIC(i,:), VOI(i), CONSTANTS, STATES(i,:));
soln = fsolve(residualfn, initialGuess_4, options);
initialGuess_4 = soln;
TEMP_ALGEBRAIC(:,55) = soln(1);
TEMP_ALGEBRAIC(:,56) = soln(2);
TEMP_ALGEBRAIC(:,57) = soln(3);
ALGEBRAIC(i,SET_ALGEBRAIC) = TEMP_ALGEBRAIC(SET_ALGEBRAIC);
end
end
end
function resid = residualSN_4(algebraicCandidate, ALGEBRAIC, VOI, CONSTANTS, STATES)
ALGEBRAIC(:,55) = algebraicCandidate(1);
ALGEBRAIC(:,56) = algebraicCandidate(2);
ALGEBRAIC(:,57) = algebraicCandidate(3);
resid(1) = ALGEBRAIC(:,55) - ( ALGEBRAIC(:,56).*ALGEBRAIC(:,57))./CONSTANTS(:,49);
resid(2) = ALGEBRAIC(:,56) - (STATES(:,9) - ALGEBRAIC(:,55));
resid(3) = ALGEBRAIC(:,57) - (CONSTANTS(:,39) - ALGEBRAIC(:,55));
end
% Compute result of a piecewise function
function x = piecewise(cases, default)
set = [0];
for i = 1:2:length(cases)
if (length(cases{i+1}) == 1)
x(cases{i} & ~set,:) = cases{i+1};
else
x(cases{i} & ~set,:) = cases{i+1}(cases{i} & ~set);
end
set = set | cases{i};
if(set), break, end
end
if (length(default) == 1)
x(~set,:) = default;
else
x(~set,:) = default(~set);
end
end
% Pad out or shorten strings to a set length
function strout = strpad(strin)
req_length = 160;
insize = size(strin,2);
if insize > req_length
strout = strin(1:req_length);
else
strout = [strin, blanks(req_length - insize)];
end
end