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 =37;
end
% There are a total of 20 entries in each of the rate and state variable arrays.
% There are a total of 144 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 (hour)');
    LEGEND_ALGEBRAIC(:,5) = strpad('V_K in component membrane (nanoM)');
    LEGEND_ALGEBRAIC(:,37) = strpad('fr in component membrane (hertz)');
    LEGEND_ALGEBRAIC(:,1) = strpad('v_sPC in component membrane (flux)');
    LEGEND_ALGEBRAIC(:,25) = strpad('Res in component membrane (gigaohm)');
    LEGEND_ALGEBRAIC(:,19) = strpad('V_rest in component membrane (milliV)');
    LEGEND_ALGEBRAIC(:,24) = strpad('V_reset in component membrane (milliV)');
    LEGEND_CONSTANTS(:,1) = strpad('V_theta in component membrane (milliV)');
    LEGEND_CONSTANTS(:,2) = strpad('Cm in component membrane (nanoF)');
    LEGEND_CONSTANTS(:,3) = strpad('R in component membrane (joule_per_mole_kelvin)');
    LEGEND_CONSTANTS(:,4) = strpad('T in component membrane (kelvin)');
    LEGEND_CONSTANTS(:,5) = strpad('T_abs in component membrane (kelvin)');
    LEGEND_CONSTANTS(:,6) = strpad('T_room in component membrane (kelvin)');
    LEGEND_CONSTANTS(:,7) = strpad('F in component membrane (coulomb_per_mole)');
    LEGEND_ALGEBRAIC(:,6) = strpad('alpha in component membrane (dimensionless)');
    LEGEND_ALGEBRAIC(:,2) = strpad('beta in component membrane (dimensionless)');
    LEGEND_ALGEBRAIC(:,16) = strpad('beta_a in component membrane (dimensionless)');
    LEGEND_ALGEBRAIC(:,17) = strpad('c in component membrane (dimensionless)');
    LEGEND_CONSTANTS(:,8) = strpad('P_Ca in component membrane (per_milliM)');
    LEGEND_ALGEBRAIC(:,3) = strpad('P_K in component membrane (per_milliM)');
    LEGEND_CONSTANTS(:,9) = strpad('P_Na in component membrane (per_milliM)');
    LEGEND_CONSTANTS(:,10) = strpad('P_Cl in component membrane (per_milliM)');
    LEGEND_CONSTANTS(:,11) = strpad('v_PK in component membrane (per_milliM)');
    LEGEND_CONSTANTS(:,12) = strpad('npk in component membrane (dimensionless)');
    LEGEND_CONSTANTS(:,13) = strpad('K_PK in component membrane (per_nanoM2)');
    LEGEND_CONSTANTS(:,14) = strpad('V_R in component membrane (gigaohm)');
    LEGEND_CONSTANTS(:,15) = strpad('K_R in component membrane (milliV)');
    LEGEND_CONSTANTS(:,16) = strpad('Ca_ex in component membrane (microM)');
    LEGEND_CONSTANTS(:,17) = strpad('Na_ex in component membrane (milliM)');
    LEGEND_CONSTANTS(:,18) = strpad('K_ex in component membrane (milliM)');
    LEGEND_CONSTANTS(:,19) = strpad('Cl_ex in component membrane (milliM)');
    LEGEND_CONSTANTS(:,143) = strpad('Na_in in component membrane (milliM)');
    LEGEND_CONSTANTS(:,140) = strpad('K_in in component membrane (milliM)');
    LEGEND_CONSTANTS(:,142) = strpad('theta_Na in component membrane (dimensionless)');
    LEGEND_CONSTANTS(:,139) = strpad('theta_K in component membrane (dimensionless)');
    LEGEND_CONSTANTS(:,20) = strpad('E_Na_0 in component membrane (milliV)');
    LEGEND_CONSTANTS(:,141) = strpad('E_Na in component membrane (milliV)');
    LEGEND_CONSTANTS(:,21) = strpad('g_Na in component membrane (nanoS)');
    LEGEND_ALGEBRAIC(:,8) = strpad('E_Ca in component membrane (milliV)');
    LEGEND_ALGEBRAIC(:,12) = strpad('g_Ca in component membrane (nanoS)');
    LEGEND_ALGEBRAIC(:,10) = strpad('g_K in component membrane (nanoS)');
    LEGEND_CONSTANTS(:,22) = strpad('g_K0 in component membrane (nanoS)');
    LEGEND_CONSTANTS(:,23) = strpad('v_gk in component membrane (nanoS)');
    LEGEND_CONSTANTS(:,24) = strpad('K_gk in component membrane (nanoM)');
    LEGEND_ALGEBRAIC(:,13) = strpad('g_KCa in component membrane (nanoS)');
    LEGEND_CONSTANTS(:,25) = strpad('v_KCa in component membrane (nanoS)');
    LEGEND_CONSTANTS(:,26) = strpad('K_KCa in component membrane (per_nanoM)');
    LEGEND_CONSTANTS(:,27) = strpad('nkca in component membrane (dimensionless)');
    LEGEND_ALGEBRAIC(:,30) = strpad('g_L in component membrane (nanoS)');
    LEGEND_CONSTANTS(:,28) = strpad('E_L_0 in component membrane (milliV)');
    LEGEND_CONSTANTS(:,144) = strpad('E_L in component membrane (milliV)');
    LEGEND_CONSTANTS(:,29) = strpad('E_K_0 in component membrane (milliV)');
    LEGEND_CONSTANTS(:,137) = strpad('E_K in component membrane (milliV)');
    LEGEND_ALGEBRAIC(:,31) = strpad('g_ex in component membrane (nanoS)');
    LEGEND_CONSTANTS(:,30) = strpad('E_ex in component membrane (milliV)');
    LEGEND_CONSTANTS(:,31) = strpad('v_ex1 in component membrane (nanoS)');
    LEGEND_CONSTANTS(:,32) = strpad('v_ex2 in component membrane (nanoS)');
    LEGEND_CONSTANTS(:,33) = strpad('K_ex1 in component membrane (picoA_2_5)');
    LEGEND_CONSTANTS(:,34) = strpad('K_ex2 in component membrane (per_microM)');
    LEGEND_CONSTANTS(:,35) = strpad('nex1 in component membrane (dimensionless)');
    LEGEND_CONSTANTS(:,36) = strpad('nex2 in component membrane (dimensionless)');
    LEGEND_CONSTANTS(:,37) = strpad('g_inhib in component membrane (nanoS)');
    LEGEND_ALGEBRAIC(:,26) = strpad('E_inhib in component membrane (milliV)');
    LEGEND_ALGEBRAIC(:,14) = strpad('GABA in component membrane (nanoM)');
    LEGEND_CONSTANTS(:,38) = strpad('GABA0 in component membrane (nanoM)');
    LEGEND_CONSTANTS(:,39) = strpad('v_GABA in component membrane (nanoM)');
    LEGEND_CONSTANTS(:,40) = strpad('K_GABA in component membrane (nanoM)');
    LEGEND_ALGEBRAIC(:,15) = strpad('Cl_in in component membrane (milliM)');
    LEGEND_CONSTANTS(:,41) = strpad('Cl0 in component membrane (milliM)');
    LEGEND_CONSTANTS(:,42) = strpad('v_Cl1 in component membrane (milliM)');
    LEGEND_CONSTANTS(:,43) = strpad('v_Cl2 in component membrane (milliM)');
    LEGEND_CONSTANTS(:,44) = strpad('K_Cl1 in component membrane (nanoM)');
    LEGEND_CONSTANTS(:,45) = strpad('K_Cl2 in component membrane (per_nanoM_0_2)');
    LEGEND_CONSTANTS(:,46) = strpad('nCl in component membrane (dimensionless)');
    LEGEND_ALGEBRAIC(:,27) = strpad('i_Na in component membrane (picoA)');
    LEGEND_ALGEBRAIC(:,29) = strpad('i_Na_abs in component membrane (picoA)');
    LEGEND_ALGEBRAIC(:,21) = strpad('i_Ca in component membrane (picoA)');
    LEGEND_ALGEBRAIC(:,22) = strpad('i_K in component membrane (picoA)');
    LEGEND_ALGEBRAIC(:,23) = strpad('i_KCa in component membrane (picoA)');
    LEGEND_ALGEBRAIC(:,32) = strpad('i_L in component membrane (picoA)');
    LEGEND_ALGEBRAIC(:,33) = strpad('i_ex in component membrane (picoA)');
    LEGEND_ALGEBRAIC(:,28) = strpad('i_inhib in component membrane (picoA)');
    LEGEND_ALGEBRAIC(:,34) = strpad('I_star in component membrane (picoA)');
    LEGEND_ALGEBRAIC(:,35) = strpad('R_star in component membrane (milliohm)');
    LEGEND_ALGEBRAIC(:,36) = strpad('tau_m in component membrane (second)');
    LEGEND_ALGEBRAIC(:,20) = strpad('theta in component membrane (milliV)');
    LEGEND_ALGEBRAIC(:,18) = strpad('psi in component membrane (dimensionless)');
    LEGEND_CONSTANTS(:,47) = strpad('k_q in component membrane (per_kelvin)');
    LEGEND_CONSTANTS(:,48) = strpad('k_MK in component membrane (nanoM)');
    LEGEND_CONSTANTS(:,49) = strpad('V_MK in component membrane (nanoM)');
    LEGEND_CONSTANTS(:,50) = strpad('V_b in component membrane (nanoM)');
    LEGEND_CONSTANTS(:,51) = strpad('k_b in component membrane (dimensionless)');
    LEGEND_CONSTANTS(:,52) = strpad('K_D in component membrane (nanoM)');
    LEGEND_CONSTANTS(:,53) = strpad('v_sP0 in component membrane (flux)');
    LEGEND_CONSTANTS(:,54) = strpad('C_T in component membrane (flux)');
    LEGEND_CONSTANTS(:,55) = strpad('nca in component membrane (dimensionless)');
    LEGEND_CONSTANTS(:,56) = strpad('v_Ca in component membrane (nanoS)');
    LEGEND_CONSTANTS(:,57) = strpad('K_Ca in component membrane (nanoM_2_2)');
    LEGEND_CONSTANTS(:,58) = strpad('K_C in component membrane (nanoM)');
    LEGEND_STATES(:,1) = strpad('Ca_in in component Ca_in (microM)');
    LEGEND_STATES(:,2) = strpad('BC in component BC (nanoM)');
    LEGEND_STATES(:,3) = strpad('MP in component MP (nanoM)');
    LEGEND_STATES(:,4) = strpad('VIP in component VIP (nanoM)');
    LEGEND_STATES(:,5) = strpad('CC in component CC (nanoM)');
    LEGEND_STATES(:,6) = strpad('CB in component CB (nanoM)');
    LEGEND_CONSTANTS(:,59) = strpad('beta_IP3 in component Ca_in (dimensionless)');
    LEGEND_STATES(:,7) = strpad('Ca_store in component Ca_in (nanoM)');
    LEGEND_ALGEBRAIC(:,7) = strpad('v0 in component Ca_in (microM_per_hour)');
    LEGEND_CONSTANTS(:,138) = strpad('v1 in component Ca_in (microM_per_hour)');
    LEGEND_ALGEBRAIC(:,9) = strpad('v2 in component Ca_in (microM_per_hour)');
    LEGEND_ALGEBRAIC(:,11) = strpad('v3 in component Ca_in (microM_per_hour)');
    LEGEND_CONSTANTS(:,60) = strpad('v_kk in component Ca_in (per_microM_per_hour)');
    LEGEND_CONSTANTS(:,61) = strpad('v_v0 in component Ca_in (microM_per_hour)');
    LEGEND_CONSTANTS(:,62) = strpad('kf in component Ca_in (per_hour)');
    LEGEND_ALGEBRAIC(:,4) = strpad('k in component Ca_in (per_microM_per_hour)');
    LEGEND_CONSTANTS(:,63) = strpad('nkk in component Ca_in (dimensionless)');
    LEGEND_CONSTANTS(:,64) = strpad('nv0 in component Ca_in (dimensionless)');
    LEGEND_CONSTANTS(:,65) = strpad('m in component Ca_in (dimensionless)');
    LEGEND_CONSTANTS(:,66) = strpad('n in component Ca_in (dimensionless)');
    LEGEND_CONSTANTS(:,67) = strpad('p in component Ca_in (dimensionless)');
    LEGEND_CONSTANTS(:,68) = strpad('v in component Ca_in (dimensionless)');
    LEGEND_CONSTANTS(:,69) = strpad('K_kk in component Ca_in (nanoM_0_1)');
    LEGEND_CONSTANTS(:,70) = strpad('K_v0 in component Ca_in (nanoM_4_5)');
    LEGEND_CONSTANTS(:,71) = strpad('K_2 in component Ca_in (microM)');
    LEGEND_CONSTANTS(:,72) = strpad('K_R in component Ca_in (microM)');
    LEGEND_CONSTANTS(:,73) = strpad('K_A in component Ca_in (microM)');
    LEGEND_CONSTANTS(:,74) = strpad('v_M1 in component Ca_in (microM_per_hour)');
    LEGEND_CONSTANTS(:,75) = strpad('v_M2 in component Ca_in (microM_per_hour)');
    LEGEND_CONSTANTS(:,76) = strpad('v_M3 in component Ca_in (microM_per_hour)');
    LEGEND_CONSTANTS(:,77) = strpad('v_VIP in component VIP (flux)');
    LEGEND_CONSTANTS(:,78) = strpad('K_VIP in component VIP (hertz_1_9)');
    LEGEND_CONSTANTS(:,79) = strpad('kd_VIP in component VIP (nanoM_0_8_per_hour)');
    LEGEND_CONSTANTS(:,80) = strpad('nd_VIP in component VIP (dimensionless)');
    LEGEND_CONSTANTS(:,81) = strpad('n_VIP in component VIP (dimensionless)');
    LEGEND_CONSTANTS(:,82) = strpad('vmP in component MP (flux)');
    LEGEND_CONSTANTS(:,83) = strpad('kdmp in component MP (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,84) = strpad('KAP in component MP (nanoM)');
    LEGEND_CONSTANTS(:,85) = strpad('KmP in component MP (nanoM)');
    LEGEND_CONSTANTS(:,86) = strpad('n in component model_parameters (dimensionless)');
    LEGEND_STATES(:,8) = strpad('BN in component BN (nanoM)');
    LEGEND_STATES(:,9) = strpad('MC in component MC (nanoM)');
    LEGEND_CONSTANTS(:,87) = strpad('vsC in component MC (flux)');
    LEGEND_CONSTANTS(:,88) = strpad('vmC in component MC (flux)');
    LEGEND_CONSTANTS(:,89) = strpad('kdmc in component MC (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,90) = strpad('KAC in component MC (nanoM)');
    LEGEND_CONSTANTS(:,91) = strpad('KmC in component MC (nanoM)');
    LEGEND_STATES(:,10) = strpad('MB in component MB (nanoM)');
    LEGEND_CONSTANTS(:,92) = strpad('vsB in component MB (flux)');
    LEGEND_CONSTANTS(:,93) = strpad('vmB in component MB (flux)');
    LEGEND_CONSTANTS(:,94) = strpad('kdmb in component MB (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,95) = strpad('KIB in component MB (nanoM)');
    LEGEND_CONSTANTS(:,96) = strpad('KmB in component MB (nanoM)');
    LEGEND_CONSTANTS(:,97) = strpad('m in component model_parameters (dimensionless)');
    LEGEND_STATES(:,11) = strpad('PC in component PC (nanoM)');
    LEGEND_CONSTANTS(:,98) = strpad('ksP in component model_parameters (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,99) = strpad('Kp in component model_parameters (nanoM)');
    LEGEND_CONSTANTS(:,100) = strpad('Kdp in component model_parameters (nanoM)');
    LEGEND_CONSTANTS(:,101) = strpad('k3 in component model_parameters (second_order_rate_constant)');
    LEGEND_CONSTANTS(:,102) = strpad('k4 in component model_parameters (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,103) = strpad('kdn in component model_parameters (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,104) = strpad('V1P in component model_parameters (flux)');
    LEGEND_CONSTANTS(:,105) = strpad('V2P in component model_parameters (flux)');
    LEGEND_STATES(:,12) = strpad('PCP in component PCP (nanoM)');
    LEGEND_STATES(:,13) = strpad('PCC in component PCC (nanoM)');
    LEGEND_CONSTANTS(:,106) = strpad('ksC in component model_parameters (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,107) = strpad('kdnc in component model_parameters (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,108) = strpad('V1C in component model_parameters (flux)');
    LEGEND_CONSTANTS(:,109) = strpad('V2C in component model_parameters (flux)');
    LEGEND_STATES(:,14) = strpad('CCP in component CCP (nanoM)');
    LEGEND_CONSTANTS(:,110) = strpad('vdPC in component model_parameters (flux)');
    LEGEND_CONSTANTS(:,111) = strpad('Kd in component model_parameters (nanoM)');
    LEGEND_CONSTANTS(:,112) = strpad('vdCC in component model_parameters (flux)');
    LEGEND_CONSTANTS(:,113) = strpad('k1 in component model_parameters (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,114) = strpad('k2 in component model_parameters (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,115) = strpad('V1PC in component model_parameters (flux)');
    LEGEND_CONSTANTS(:,116) = strpad('V2PC in component model_parameters (flux)');
    LEGEND_STATES(:,15) = strpad('PCCP in component PCCP (nanoM)');
    LEGEND_STATES(:,16) = strpad('PCN in component PCN (nanoM)');
    LEGEND_CONSTANTS(:,117) = strpad('k7 in component model_parameters (second_order_rate_constant)');
    LEGEND_CONSTANTS(:,118) = strpad('k8 in component model_parameters (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,119) = strpad('V3PC in component model_parameters (flux)');
    LEGEND_CONSTANTS(:,120) = strpad('V4PC in component model_parameters (flux)');
    LEGEND_STATES(:,17) = strpad('PCNP in component PCNP (nanoM)');
    LEGEND_STATES(:,18) = strpad('IN in component IN (nanoM)');
    LEGEND_CONSTANTS(:,121) = strpad('vdPCC in component model_parameters (flux)');
    LEGEND_CONSTANTS(:,122) = strpad('vdPCN in component model_parameters (flux)');
    LEGEND_CONSTANTS(:,123) = strpad('ksB in component model_parameters (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,124) = strpad('k5 in component model_parameters (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,125) = strpad('k6 in component model_parameters (first_order_rate_constant)');
    LEGEND_CONSTANTS(:,126) = strpad('V1B in component model_parameters (flux)');
    LEGEND_CONSTANTS(:,127) = strpad('V2B in component model_parameters (flux)');
    LEGEND_STATES(:,19) = strpad('BCP in component BCP (nanoM)');
    LEGEND_CONSTANTS(:,128) = strpad('vdBC in component model_parameters (flux)');
    LEGEND_CONSTANTS(:,129) = strpad('V3B in component model_parameters (flux)');
    LEGEND_CONSTANTS(:,130) = strpad('V4B in component model_parameters (flux)');
    LEGEND_STATES(:,20) = strpad('BNP in component BNP (nanoM)');
    LEGEND_CONSTANTS(:,131) = strpad('vdBN in component model_parameters (flux)');
    LEGEND_CONSTANTS(:,132) = strpad('vdIN in component model_parameters (flux)');
    LEGEND_CONSTANTS(:,133) = strpad('K1CB in component CB (nanoM)');
    LEGEND_CONSTANTS(:,134) = strpad('K2CB in component CB (nanoM)');
    LEGEND_CONSTANTS(:,135) = strpad('WT in component CB (hour)');
    LEGEND_CONSTANTS(:,136) = strpad('vP in component CB (nanoM)');
    LEGEND_RATES(:,1) = strpad('d/dt Ca_in in component Ca_in (microM)');
    LEGEND_RATES(:,7) = strpad('d/dt Ca_store in component Ca_in (nanoM)');
    LEGEND_RATES(:,4) = strpad('d/dt VIP in component VIP (nanoM)');
    LEGEND_RATES(:,3) = strpad('d/dt MP in component MP (nanoM)');
    LEGEND_RATES(:,9) = strpad('d/dt MC in component MC (nanoM)');
    LEGEND_RATES(:,10) = strpad('d/dt MB in component MB (nanoM)');
    LEGEND_RATES(:,11) = strpad('d/dt PC in component PC (nanoM)');
    LEGEND_RATES(:,5) = strpad('d/dt CC in component CC (nanoM)');
    LEGEND_RATES(:,12) = strpad('d/dt PCP in component PCP (nanoM)');
    LEGEND_RATES(:,14) = strpad('d/dt CCP in component CCP (nanoM)');
    LEGEND_RATES(:,13) = strpad('d/dt PCC in component PCC (nanoM)');
    LEGEND_RATES(:,16) = strpad('d/dt PCN in component PCN (nanoM)');
    LEGEND_RATES(:,15) = strpad('d/dt PCCP in component PCCP (nanoM)');
    LEGEND_RATES(:,17) = strpad('d/dt PCNP in component PCNP (nanoM)');
    LEGEND_RATES(:,2) = strpad('d/dt BC in component BC (nanoM)');
    LEGEND_RATES(:,19) = strpad('d/dt BCP in component BCP (nanoM)');
    LEGEND_RATES(:,8) = strpad('d/dt BN in component BN (nanoM)');
    LEGEND_RATES(:,20) = strpad('d/dt BNP in component BNP (nanoM)');
    LEGEND_RATES(:,18) = strpad('d/dt IN in component IN (nanoM)');
    LEGEND_RATES(:,6) = strpad('d/dt CB in component CB (nanoM)');
    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) = 20.0;
    CONSTANTS(:,2) = 5.0;
    CONSTANTS(:,3) = 8.314;
    CONSTANTS(:,4) = 37;
    CONSTANTS(:,5) = 273.15;
    CONSTANTS(:,6) = 22;
    CONSTANTS(:,7) = 96485;
    CONSTANTS(:,8) = 0.05;
    CONSTANTS(:,9) = 0.036;
    CONSTANTS(:,10) = 0.3;
    CONSTANTS(:,11) = 1.9;
    CONSTANTS(:,12) = -2.0;
    CONSTANTS(:,13) = 1.0;
    CONSTANTS(:,14) = 0.41;
    CONSTANTS(:,15) = 34.0;
    CONSTANTS(:,16) = 5.0;
    CONSTANTS(:,17) = 145.0;
    CONSTANTS(:,18) = 1.0;
    CONSTANTS(:,19) = 114.5;
    CONSTANTS(:,20) = 45.0;
    CONSTANTS(:,21) = 36.0;
    CONSTANTS(:,22) = 9.7;
    CONSTANTS(:,23) = 10.0;
    CONSTANTS(:,24) = 10.0;
    CONSTANTS(:,25) = 3.0;
    CONSTANTS(:,26) = 0.16;
    CONSTANTS(:,27) = -1.0;
    CONSTANTS(:,28) = -29.0;
    CONSTANTS(:,29) = -97.0;
    CONSTANTS(:,30) = 0.0;
    CONSTANTS(:,31) = 105.0;
    CONSTANTS(:,32) = 4.4;
    CONSTANTS(:,33) = 574050000.0;
    CONSTANTS(:,34) = 1.0;
    CONSTANTS(:,35) = 2.5;
    CONSTANTS(:,36) = -1.0;
    CONSTANTS(:,37) = 12.3;
    CONSTANTS(:,38) = 0.2;
    CONSTANTS(:,39) = 19.0;
    CONSTANTS(:,40) = 3.0;
    CONSTANTS(:,41) = 1.0;
    CONSTANTS(:,42) = 15.5;
    CONSTANTS(:,43) = 19.0;
    CONSTANTS(:,44) = 4.0;
    CONSTANTS(:,45) = 1.0;
    CONSTANTS(:,46) = -0.2;
    CONSTANTS(:,47) = 8.75e-5;
    CONSTANTS(:,48) = 2.9;
    CONSTANTS(:,49) = 5.0;
    CONSTANTS(:,50) = 2.0;
    CONSTANTS(:,51) = 2.0;
    CONSTANTS(:,52) = 0.08;
    CONSTANTS(:,53) = 1.0;
    CONSTANTS(:,54) = 1.6;
    CONSTANTS(:,55) = 2.2;
    CONSTANTS(:,56) = 12.3;
    CONSTANTS(:,57) = 22.0;
    CONSTANTS(:,58) = 0.15;
    STATES(:,1) = 0.1;
    STATES(:,2) = 2.41;
    STATES(:,3) = 2.8;
    STATES(:,4) = 0.0;
    STATES(:,5) = 12.0;
    STATES(:,6) = 0.12;
    CONSTANTS(:,59) = 0.5;
    STATES(:,7) = 0.1;
    CONSTANTS(:,60) = 3.3;
    CONSTANTS(:,61) = 0.09;
    CONSTANTS(:,62) = 0.001;
    CONSTANTS(:,63) = 0.1;
    CONSTANTS(:,64) = 4.5;
    CONSTANTS(:,65) = 6.0;
    CONSTANTS(:,66) = 2.2;
    CONSTANTS(:,67) = 4.2;
    CONSTANTS(:,68) = 2.0;
    CONSTANTS(:,69) = 0.02;
    CONSTANTS(:,70) = 4.5;
    CONSTANTS(:,71) = 5.0;
    CONSTANTS(:,72) = 3.0;
    CONSTANTS(:,73) = 0.67;
    CONSTANTS(:,74) = 0.0003;
    CONSTANTS(:,75) = 149.5;
    CONSTANTS(:,76) = 400.0;
    CONSTANTS(:,77) = 0.5;
    CONSTANTS(:,78) = 15.0;
    CONSTANTS(:,79) = 0.5;
    CONSTANTS(:,80) = 0.2;
    CONSTANTS(:,81) = 1.9;
    CONSTANTS(:,82) = 1.1;
    CONSTANTS(:,83) = 0.01;
    CONSTANTS(:,84) = 0.6;
    CONSTANTS(:,85) = 0.31;
    CONSTANTS(:,86) = 4.0;
    STATES(:,8) = 1.94;
    STATES(:,9) = 2.0;
    CONSTANTS(:,87) = 1.1;
    CONSTANTS(:,88) = 1.0;
    CONSTANTS(:,89) = 0.01;
    CONSTANTS(:,90) = 0.6;
    CONSTANTS(:,91) = 0.4;
    STATES(:,10) = 7.94;
    CONSTANTS(:,92) = 1.0;
    CONSTANTS(:,93) = 0.8;
    CONSTANTS(:,94) = 0.01;
    CONSTANTS(:,95) = 2.2;
    CONSTANTS(:,96) = 0.4;
    CONSTANTS(:,97) = 2.0;
    STATES(:,11) = 0.4;
    CONSTANTS(:,98) = 0.6;
    CONSTANTS(:,99) = 0.1;
    CONSTANTS(:,100) = 0.1;
    CONSTANTS(:,101) = 0.4;
    CONSTANTS(:,102) = 0.2;
    CONSTANTS(:,103) = 0.01;
    CONSTANTS(:,104) = 0.4;
    CONSTANTS(:,105) = 0.3;
    STATES(:,12) = 0.13;
    STATES(:,13) = 1.26;
    CONSTANTS(:,106) = 1.6;
    CONSTANTS(:,107) = 0.12;
    CONSTANTS(:,108) = 0.6;
    CONSTANTS(:,109) = 0.1;
    STATES(:,14) = 9.0;
    CONSTANTS(:,110) = 0.7;
    CONSTANTS(:,111) = 0.3;
    CONSTANTS(:,112) = 0.7;
    CONSTANTS(:,113) = 0.45;
    CONSTANTS(:,114) = 0.2;
    CONSTANTS(:,115) = 0.4;
    CONSTANTS(:,116) = 0.1;
    STATES(:,15) = 0.2;
    STATES(:,16) = 0.16;
    CONSTANTS(:,117) = 0.5;
    CONSTANTS(:,118) = 0.1;
    CONSTANTS(:,119) = 0.4;
    CONSTANTS(:,120) = 0.1;
    STATES(:,17) = 0.091;
    STATES(:,18) = 0.05;
    CONSTANTS(:,121) = 0.7;
    CONSTANTS(:,122) = 0.7;
    CONSTANTS(:,123) = 0.12;
    CONSTANTS(:,124) = 0.4;
    CONSTANTS(:,125) = 0.2;
    CONSTANTS(:,126) = 0.5;
    CONSTANTS(:,127) = 0.1;
    STATES(:,19) = 0.48;
    CONSTANTS(:,128) = 0.5;
    CONSTANTS(:,129) = 0.5;
    CONSTANTS(:,130) = 0.2;
    STATES(:,20) = 0.32;
    CONSTANTS(:,131) = 0.6;
    CONSTANTS(:,132) = 0.8;
    CONSTANTS(:,133) = 0.01;
    CONSTANTS(:,134) = 0.01;
    CONSTANTS(:,135) = 1.0;
    CONSTANTS(:,136) = 1.0;
    CONSTANTS(:,137) = ( CONSTANTS(:,29).*(CONSTANTS(:,4)+CONSTANTS(:,5)))./(CONSTANTS(:,6)+CONSTANTS(:,5));
    CONSTANTS(:,138) =  CONSTANTS(:,74).*CONSTANTS(:,59);
    CONSTANTS(:,139) = exp(CONSTANTS(:,137)./( CONSTANTS(:,47).*(CONSTANTS(:,4)+CONSTANTS(:,5)).*1000.00));
    CONSTANTS(:,140) = CONSTANTS(:,18)./CONSTANTS(:,139);
    CONSTANTS(:,141) = ( CONSTANTS(:,20).*(CONSTANTS(:,4)+CONSTANTS(:,5)))./(CONSTANTS(:,6)+CONSTANTS(:,5));
    CONSTANTS(:,142) = exp(CONSTANTS(:,141)./( CONSTANTS(:,47).*(CONSTANTS(:,4)+CONSTANTS(:,5)).*1000.00));
    CONSTANTS(:,143) = CONSTANTS(:,17)./CONSTANTS(:,142);
    CONSTANTS(:,144) = ( CONSTANTS(:,28).*(CONSTANTS(:,4)+CONSTANTS(:,5)))./(CONSTANTS(:,6)+CONSTANTS(:,5));
    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
    RATES(:,9) =  CONSTANTS(:,87).*(power(STATES(:,8), CONSTANTS(:,86))./(power(CONSTANTS(:,90), CONSTANTS(:,86))+power(STATES(:,8), CONSTANTS(:,86)))) - ( CONSTANTS(:,88).*(STATES(:,9)./(CONSTANTS(:,91)+STATES(:,9)))+ CONSTANTS(:,89).*STATES(:,9));
    RATES(:,10) =  CONSTANTS(:,92).*(power(CONSTANTS(:,95), CONSTANTS(:,97))./(power(CONSTANTS(:,95), CONSTANTS(:,97))+power(STATES(:,8), CONSTANTS(:,97)))) - ( CONSTANTS(:,93).*(STATES(:,10)./(CONSTANTS(:,96)+STATES(:,10)))+ CONSTANTS(:,94).*STATES(:,10));
    RATES(:,11) = ( CONSTANTS(:,98).*STATES(:,3)+ CONSTANTS(:,105).*(STATES(:,12)./(CONSTANTS(:,100)+STATES(:,12)))+ CONSTANTS(:,102).*STATES(:,13)) - ( CONSTANTS(:,104).*(STATES(:,11)./(CONSTANTS(:,99)+STATES(:,11)))+ CONSTANTS(:,101).*STATES(:,11).*STATES(:,5)+ CONSTANTS(:,103).*STATES(:,11));
    RATES(:,5) = ( CONSTANTS(:,106).*STATES(:,9)+ CONSTANTS(:,109).*(STATES(:,14)./(CONSTANTS(:,100)+STATES(:,14)))+ CONSTANTS(:,102).*STATES(:,13)) - ( CONSTANTS(:,108).*(STATES(:,5)./(CONSTANTS(:,99)+STATES(:,5)))+ CONSTANTS(:,101).*STATES(:,11).*STATES(:,5)+ CONSTANTS(:,107).*STATES(:,5));
    RATES(:,12) =  CONSTANTS(:,104).*(STATES(:,11)./(CONSTANTS(:,99)+STATES(:,11))) - ( CONSTANTS(:,105).*(STATES(:,12)./(CONSTANTS(:,100)+STATES(:,12)))+ CONSTANTS(:,110).*(STATES(:,12)./(CONSTANTS(:,111)+STATES(:,12)))+ CONSTANTS(:,103).*STATES(:,12));
    RATES(:,14) =  CONSTANTS(:,108).*(STATES(:,5)./(CONSTANTS(:,99)+STATES(:,5))) - ( CONSTANTS(:,109).*(STATES(:,14)./(CONSTANTS(:,100)+STATES(:,14)))+ CONSTANTS(:,112).*(STATES(:,14)./(CONSTANTS(:,111)+STATES(:,14)))+ CONSTANTS(:,103).*STATES(:,14));
    RATES(:,13) = ( CONSTANTS(:,116).*(STATES(:,15)./(CONSTANTS(:,100)+STATES(:,15)))+ CONSTANTS(:,101).*STATES(:,11).*STATES(:,5)+ CONSTANTS(:,114).*STATES(:,16)) - ( CONSTANTS(:,115).*(STATES(:,13)./(CONSTANTS(:,99)+STATES(:,13)))+ CONSTANTS(:,102).*STATES(:,13)+ CONSTANTS(:,113).*STATES(:,13)+ CONSTANTS(:,103).*STATES(:,13));
    RATES(:,16) = ( CONSTANTS(:,120).*(STATES(:,17)./(CONSTANTS(:,100)+STATES(:,17)))+ CONSTANTS(:,113).*STATES(:,13)+ CONSTANTS(:,118).*STATES(:,18)) - ( CONSTANTS(:,119).*(STATES(:,16)./(CONSTANTS(:,99)+STATES(:,16)))+ CONSTANTS(:,114).*STATES(:,16)+ CONSTANTS(:,117).*STATES(:,8).*STATES(:,16)+ CONSTANTS(:,103).*STATES(:,16));
    RATES(:,15) =  CONSTANTS(:,115).*(STATES(:,13)./(CONSTANTS(:,99)+STATES(:,13))) - ( CONSTANTS(:,116).*(STATES(:,15)./(CONSTANTS(:,100)+STATES(:,15)))+ CONSTANTS(:,121).*(STATES(:,15)./(CONSTANTS(:,111)+STATES(:,15)))+ CONSTANTS(:,103).*STATES(:,15));
    RATES(:,17) =  CONSTANTS(:,119).*(STATES(:,16)./(CONSTANTS(:,99)+STATES(:,16))) - ( CONSTANTS(:,120).*(STATES(:,17)./(CONSTANTS(:,100)+STATES(:,17)))+ CONSTANTS(:,122).*(STATES(:,17)./(CONSTANTS(:,111)+STATES(:,17)))+ CONSTANTS(:,103).*STATES(:,17));
    RATES(:,2) = ( CONSTANTS(:,127).*(STATES(:,19)./(CONSTANTS(:,100)+STATES(:,19)))+ CONSTANTS(:,125).*STATES(:,8)+ CONSTANTS(:,123).*STATES(:,10)) - ( CONSTANTS(:,126).*(STATES(:,2)./(CONSTANTS(:,99)+STATES(:,2)))+ CONSTANTS(:,124).*STATES(:,2)+ CONSTANTS(:,103).*STATES(:,2));
    RATES(:,19) =  CONSTANTS(:,126).*(STATES(:,2)./(CONSTANTS(:,99)+STATES(:,2))) - ( CONSTANTS(:,127).*(STATES(:,19)./(CONSTANTS(:,100)+STATES(:,19)))+ CONSTANTS(:,128).*(STATES(:,19)./(CONSTANTS(:,111)+STATES(:,19)))+ CONSTANTS(:,103).*STATES(:,19));
    RATES(:,8) = ( CONSTANTS(:,130).*(STATES(:,20)./(CONSTANTS(:,100)+STATES(:,20)))+ CONSTANTS(:,124).*STATES(:,2)+ CONSTANTS(:,118).*STATES(:,18)) - ( CONSTANTS(:,129).*(STATES(:,8)./(CONSTANTS(:,99)+STATES(:,8)))+ CONSTANTS(:,125).*STATES(:,8)+ CONSTANTS(:,117).*STATES(:,8).*STATES(:,16)+ CONSTANTS(:,103).*STATES(:,8));
    RATES(:,20) =  CONSTANTS(:,129).*(STATES(:,8)./(CONSTANTS(:,99)+STATES(:,8))) - ( CONSTANTS(:,130).*(STATES(:,20)./(CONSTANTS(:,100)+STATES(:,20)))+ CONSTANTS(:,131).*(STATES(:,20)./(CONSTANTS(:,111)+STATES(:,20)))+ CONSTANTS(:,103).*STATES(:,20));
    RATES(:,18) =  CONSTANTS(:,117).*STATES(:,8).*STATES(:,16) - ( CONSTANTS(:,118).*STATES(:,18)+ CONSTANTS(:,132).*(STATES(:,18)./(CONSTANTS(:,111)+STATES(:,18)))+ CONSTANTS(:,103).*STATES(:,18));
    ALGEBRAIC(:,1) = CONSTANTS(:,53)+( CONSTANTS(:,54).*STATES(:,6))./(CONSTANTS(:,58)+STATES(:,6));
    RATES(:,3) =  ALGEBRAIC(:,1).*(power(STATES(:,8), CONSTANTS(:,86))./(power(CONSTANTS(:,84), CONSTANTS(:,86))+power(STATES(:,8), CONSTANTS(:,86)))) - ( CONSTANTS(:,82).*(STATES(:,3)./(CONSTANTS(:,85)+STATES(:,3)))+ CONSTANTS(:,83).*STATES(:,3));
    ALGEBRAIC(:,2) = STATES(:,4)./(STATES(:,4)+CONSTANTS(:,52));
    ALGEBRAIC(:,5) = ( CONSTANTS(:,49).*STATES(:,1))./(CONSTANTS(:,48)+STATES(:,1))+( CONSTANTS(:,50).*ALGEBRAIC(:,2))./(CONSTANTS(:,51)+ALGEBRAIC(:,2));
    RATES(:,6) = (( ALGEBRAIC(:,5).*(1.00000 - STATES(:,6)))./((CONSTANTS(:,133)+1.00000) - STATES(:,6)) - ( CONSTANTS(:,136).*STATES(:,6))./(CONSTANTS(:,134)+STATES(:,6)))./CONSTANTS(:,135);
    ALGEBRAIC(:,7) =  CONSTANTS(:,61).*(power(STATES(:,2), CONSTANTS(:,64))./(CONSTANTS(:,70)+power(STATES(:,2), CONSTANTS(:,64))));
    ALGEBRAIC(:,9) =  CONSTANTS(:,75).*(power(STATES(:,1), CONSTANTS(:,66))./(power(CONSTANTS(:,71), CONSTANTS(:,66))+power(STATES(:,1), CONSTANTS(:,66))));
    ALGEBRAIC(:,11) =  CONSTANTS(:,76).*(power(STATES(:,7), CONSTANTS(:,65))./(power(CONSTANTS(:,72), CONSTANTS(:,65))+power(STATES(:,7), CONSTANTS(:,65)))).*(power(STATES(:,1), CONSTANTS(:,67))./(power(CONSTANTS(:,73), CONSTANTS(:,67))+power(STATES(:,1), CONSTANTS(:,67))));
    ALGEBRAIC(:,4) =  CONSTANTS(:,60).*(power(STATES(:,5), CONSTANTS(:,63))./(CONSTANTS(:,69)+power(STATES(:,5), CONSTANTS(:,63))));
    RATES(:,1) = (ALGEBRAIC(:,7)+CONSTANTS(:,138)+ALGEBRAIC(:,11)+ CONSTANTS(:,62).*STATES(:,7)) - (ALGEBRAIC(:,9)+ ALGEBRAIC(:,4).*power(STATES(:,1), CONSTANTS(:,68)));
    RATES(:,7) = ALGEBRAIC(:,9) - (ALGEBRAIC(:,11)+ CONSTANTS(:,62).*STATES(:,7));
    ALGEBRAIC(:,3) =  CONSTANTS(:,11).*(power(STATES(:,2), CONSTANTS(:,12))./(CONSTANTS(:,13)+power(STATES(:,2), CONSTANTS(:,12))));
    ALGEBRAIC(:,6) =  4.00000.*CONSTANTS(:,8).*STATES(:,1).*0.00100000+ ALGEBRAIC(:,3).*CONSTANTS(:,140)+ CONSTANTS(:,9).*CONSTANTS(:,143)+ CONSTANTS(:,10).*CONSTANTS(:,19);
    ALGEBRAIC(:,14) = CONSTANTS(:,38)+( CONSTANTS(:,39).*STATES(:,4))./(CONSTANTS(:,40)+STATES(:,4));
    ALGEBRAIC(:,15) = CONSTANTS(:,41)+ CONSTANTS(:,42).*(STATES(:,3)./(CONSTANTS(:,44)+STATES(:,3)))+ CONSTANTS(:,43).*(power(ALGEBRAIC(:,14), CONSTANTS(:,46))./(CONSTANTS(:,45)+power(ALGEBRAIC(:,14), CONSTANTS(:,46))));
    ALGEBRAIC(:,16) = ( ALGEBRAIC(:,3).*CONSTANTS(:,140)+ CONSTANTS(:,9).*CONSTANTS(:,143)+ CONSTANTS(:,10).*CONSTANTS(:,19)) - ( ALGEBRAIC(:,3).*CONSTANTS(:,18)+ CONSTANTS(:,9).*CONSTANTS(:,17)+ CONSTANTS(:,10).*ALGEBRAIC(:,15));
    ALGEBRAIC(:,17) =  - ( ALGEBRAIC(:,3).*CONSTANTS(:,18)+ 4.00000.*CONSTANTS(:,8).*CONSTANTS(:,16).*0.00100000+ CONSTANTS(:,9).*CONSTANTS(:,17)+ CONSTANTS(:,10).*ALGEBRAIC(:,15));
    ALGEBRAIC(:,18) = (power((power(ALGEBRAIC(:,16), 2.00000) -  4.00000.*ALGEBRAIC(:,6).*ALGEBRAIC(:,17)), 1.0 ./ 2) - ALGEBRAIC(:,16))./( 2.00000.*ALGEBRAIC(:,6));
    ALGEBRAIC(:,19) =  (( CONSTANTS(:,3).*(CONSTANTS(:,4)+CONSTANTS(:,5)))./CONSTANTS(:,7)).*log(ALGEBRAIC(:,18)).*1000.00;
    ALGEBRAIC(:,24) = ALGEBRAIC(:,19)+4.00000;
    ALGEBRAIC(:,8) =  (( CONSTANTS(:,47).*(CONSTANTS(:,4)+CONSTANTS(:,5)))./2.00000).*log(CONSTANTS(:,16)./STATES(:,1)).*1000.00;
    ALGEBRAIC(:,12) =  CONSTANTS(:,56).*(power(STATES(:,3), CONSTANTS(:,55))./(CONSTANTS(:,57)+power(STATES(:,3), CONSTANTS(:,55))));
    ALGEBRAIC(:,10) = CONSTANTS(:,22)+ CONSTANTS(:,23).*(STATES(:,3)./(CONSTANTS(:,24)+STATES(:,3)));
    ALGEBRAIC(:,13) =  CONSTANTS(:,25).*(power(STATES(:,5), CONSTANTS(:,27))./(CONSTANTS(:,26)+power(STATES(:,5), CONSTANTS(:,27))));
    ALGEBRAIC(:,25) =  CONSTANTS(:,14).*(ALGEBRAIC(:,19)./(CONSTANTS(:,15)+ALGEBRAIC(:,19)));
    ALGEBRAIC(:,30) = 1.00000./ALGEBRAIC(:,25);
    ALGEBRAIC(:,27) =  CONSTANTS(:,21).*(ALGEBRAIC(:,19) - CONSTANTS(:,141));
    ALGEBRAIC(:,29) = power(power(ALGEBRAIC(:,27), 2.00000), 1.0 ./ 2);
    ALGEBRAIC(:,31) =  CONSTANTS(:,31).*(power(ALGEBRAIC(:,29), CONSTANTS(:,35))./(CONSTANTS(:,33)+power(ALGEBRAIC(:,29), CONSTANTS(:,35))))+ CONSTANTS(:,32).*(power(STATES(:,1), CONSTANTS(:,36))./(CONSTANTS(:,34)+power(STATES(:,1), CONSTANTS(:,36))));
    ALGEBRAIC(:,26) =   - CONSTANTS(:,47).*(CONSTANTS(:,4)+CONSTANTS(:,5)).*log(CONSTANTS(:,19)./ALGEBRAIC(:,15)).*1000.00;
    ALGEBRAIC(:,34) = ( CONSTANTS(:,21).*CONSTANTS(:,141)+ ALGEBRAIC(:,12).*ALGEBRAIC(:,8)+ ALGEBRAIC(:,10).*CONSTANTS(:,137)+ ALGEBRAIC(:,30).*CONSTANTS(:,144)+ ALGEBRAIC(:,13).*CONSTANTS(:,137)) - ( CONSTANTS(:,37).*ALGEBRAIC(:,26)+ ALGEBRAIC(:,31).*CONSTANTS(:,30));
    ALGEBRAIC(:,35) = 1.00000./((CONSTANTS(:,21)+ALGEBRAIC(:,10)+ALGEBRAIC(:,30)+ALGEBRAIC(:,12)+ALGEBRAIC(:,13)) - (CONSTANTS(:,37)+ALGEBRAIC(:,31)));
    ALGEBRAIC(:,36) =  CONSTANTS(:,2).*ALGEBRAIC(:,35);
    ALGEBRAIC(:,20) = ALGEBRAIC(:,19)+CONSTANTS(:,1);
    ALGEBRAIC(:,37) = -1.00000./( ALGEBRAIC(:,36).*log((ALGEBRAIC(:,20) -  ALGEBRAIC(:,35).*ALGEBRAIC(:,34))./(ALGEBRAIC(:,24) -  ALGEBRAIC(:,35).*ALGEBRAIC(:,34))));
    RATES(:,4) =  CONSTANTS(:,77).*(power(ALGEBRAIC(:,37), CONSTANTS(:,81))./(CONSTANTS(:,78)+power(ALGEBRAIC(:,37), CONSTANTS(:,81)))) -  CONSTANTS(:,79).*power(STATES(:,4), CONSTANTS(:,80));
   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(:,1) = CONSTANTS(:,53)+( CONSTANTS(:,54).*STATES(:,6))./(CONSTANTS(:,58)+STATES(:,6));
    ALGEBRAIC(:,2) = STATES(:,4)./(STATES(:,4)+CONSTANTS(:,52));
    ALGEBRAIC(:,5) = ( CONSTANTS(:,49).*STATES(:,1))./(CONSTANTS(:,48)+STATES(:,1))+( CONSTANTS(:,50).*ALGEBRAIC(:,2))./(CONSTANTS(:,51)+ALGEBRAIC(:,2));
    ALGEBRAIC(:,7) =  CONSTANTS(:,61).*(power(STATES(:,2), CONSTANTS(:,64))./(CONSTANTS(:,70)+power(STATES(:,2), CONSTANTS(:,64))));
    ALGEBRAIC(:,9) =  CONSTANTS(:,75).*(power(STATES(:,1), CONSTANTS(:,66))./(power(CONSTANTS(:,71), CONSTANTS(:,66))+power(STATES(:,1), CONSTANTS(:,66))));
    ALGEBRAIC(:,11) =  CONSTANTS(:,76).*(power(STATES(:,7), CONSTANTS(:,65))./(power(CONSTANTS(:,72), CONSTANTS(:,65))+power(STATES(:,7), CONSTANTS(:,65)))).*(power(STATES(:,1), CONSTANTS(:,67))./(power(CONSTANTS(:,73), CONSTANTS(:,67))+power(STATES(:,1), CONSTANTS(:,67))));
    ALGEBRAIC(:,4) =  CONSTANTS(:,60).*(power(STATES(:,5), CONSTANTS(:,63))./(CONSTANTS(:,69)+power(STATES(:,5), CONSTANTS(:,63))));
    ALGEBRAIC(:,3) =  CONSTANTS(:,11).*(power(STATES(:,2), CONSTANTS(:,12))./(CONSTANTS(:,13)+power(STATES(:,2), CONSTANTS(:,12))));
    ALGEBRAIC(:,6) =  4.00000.*CONSTANTS(:,8).*STATES(:,1).*0.00100000+ ALGEBRAIC(:,3).*CONSTANTS(:,140)+ CONSTANTS(:,9).*CONSTANTS(:,143)+ CONSTANTS(:,10).*CONSTANTS(:,19);
    ALGEBRAIC(:,14) = CONSTANTS(:,38)+( CONSTANTS(:,39).*STATES(:,4))./(CONSTANTS(:,40)+STATES(:,4));
    ALGEBRAIC(:,15) = CONSTANTS(:,41)+ CONSTANTS(:,42).*(STATES(:,3)./(CONSTANTS(:,44)+STATES(:,3)))+ CONSTANTS(:,43).*(power(ALGEBRAIC(:,14), CONSTANTS(:,46))./(CONSTANTS(:,45)+power(ALGEBRAIC(:,14), CONSTANTS(:,46))));
    ALGEBRAIC(:,16) = ( ALGEBRAIC(:,3).*CONSTANTS(:,140)+ CONSTANTS(:,9).*CONSTANTS(:,143)+ CONSTANTS(:,10).*CONSTANTS(:,19)) - ( ALGEBRAIC(:,3).*CONSTANTS(:,18)+ CONSTANTS(:,9).*CONSTANTS(:,17)+ CONSTANTS(:,10).*ALGEBRAIC(:,15));
    ALGEBRAIC(:,17) =  - ( ALGEBRAIC(:,3).*CONSTANTS(:,18)+ 4.00000.*CONSTANTS(:,8).*CONSTANTS(:,16).*0.00100000+ CONSTANTS(:,9).*CONSTANTS(:,17)+ CONSTANTS(:,10).*ALGEBRAIC(:,15));
    ALGEBRAIC(:,18) = (power((power(ALGEBRAIC(:,16), 2.00000) -  4.00000.*ALGEBRAIC(:,6).*ALGEBRAIC(:,17)), 1.0 ./ 2) - ALGEBRAIC(:,16))./( 2.00000.*ALGEBRAIC(:,6));
    ALGEBRAIC(:,19) =  (( CONSTANTS(:,3).*(CONSTANTS(:,4)+CONSTANTS(:,5)))./CONSTANTS(:,7)).*log(ALGEBRAIC(:,18)).*1000.00;
    ALGEBRAIC(:,24) = ALGEBRAIC(:,19)+4.00000;
    ALGEBRAIC(:,8) =  (( CONSTANTS(:,47).*(CONSTANTS(:,4)+CONSTANTS(:,5)))./2.00000).*log(CONSTANTS(:,16)./STATES(:,1)).*1000.00;
    ALGEBRAIC(:,12) =  CONSTANTS(:,56).*(power(STATES(:,3), CONSTANTS(:,55))./(CONSTANTS(:,57)+power(STATES(:,3), CONSTANTS(:,55))));
    ALGEBRAIC(:,10) = CONSTANTS(:,22)+ CONSTANTS(:,23).*(STATES(:,3)./(CONSTANTS(:,24)+STATES(:,3)));
    ALGEBRAIC(:,13) =  CONSTANTS(:,25).*(power(STATES(:,5), CONSTANTS(:,27))./(CONSTANTS(:,26)+power(STATES(:,5), CONSTANTS(:,27))));
    ALGEBRAIC(:,25) =  CONSTANTS(:,14).*(ALGEBRAIC(:,19)./(CONSTANTS(:,15)+ALGEBRAIC(:,19)));
    ALGEBRAIC(:,30) = 1.00000./ALGEBRAIC(:,25);
    ALGEBRAIC(:,27) =  CONSTANTS(:,21).*(ALGEBRAIC(:,19) - CONSTANTS(:,141));
    ALGEBRAIC(:,29) = power(power(ALGEBRAIC(:,27), 2.00000), 1.0 ./ 2);
    ALGEBRAIC(:,31) =  CONSTANTS(:,31).*(power(ALGEBRAIC(:,29), CONSTANTS(:,35))./(CONSTANTS(:,33)+power(ALGEBRAIC(:,29), CONSTANTS(:,35))))+ CONSTANTS(:,32).*(power(STATES(:,1), CONSTANTS(:,36))./(CONSTANTS(:,34)+power(STATES(:,1), CONSTANTS(:,36))));
    ALGEBRAIC(:,26) =   - CONSTANTS(:,47).*(CONSTANTS(:,4)+CONSTANTS(:,5)).*log(CONSTANTS(:,19)./ALGEBRAIC(:,15)).*1000.00;
    ALGEBRAIC(:,34) = ( CONSTANTS(:,21).*CONSTANTS(:,141)+ ALGEBRAIC(:,12).*ALGEBRAIC(:,8)+ ALGEBRAIC(:,10).*CONSTANTS(:,137)+ ALGEBRAIC(:,30).*CONSTANTS(:,144)+ ALGEBRAIC(:,13).*CONSTANTS(:,137)) - ( CONSTANTS(:,37).*ALGEBRAIC(:,26)+ ALGEBRAIC(:,31).*CONSTANTS(:,30));
    ALGEBRAIC(:,35) = 1.00000./((CONSTANTS(:,21)+ALGEBRAIC(:,10)+ALGEBRAIC(:,30)+ALGEBRAIC(:,12)+ALGEBRAIC(:,13)) - (CONSTANTS(:,37)+ALGEBRAIC(:,31)));
    ALGEBRAIC(:,36) =  CONSTANTS(:,2).*ALGEBRAIC(:,35);
    ALGEBRAIC(:,20) = ALGEBRAIC(:,19)+CONSTANTS(:,1);
    ALGEBRAIC(:,37) = -1.00000./( ALGEBRAIC(:,36).*log((ALGEBRAIC(:,20) -  ALGEBRAIC(:,35).*ALGEBRAIC(:,34))./(ALGEBRAIC(:,24) -  ALGEBRAIC(:,35).*ALGEBRAIC(:,34))));
    ALGEBRAIC(:,21) =  ALGEBRAIC(:,12).*(ALGEBRAIC(:,19) - ALGEBRAIC(:,8));
    ALGEBRAIC(:,22) =  ALGEBRAIC(:,10).*(ALGEBRAIC(:,19) - CONSTANTS(:,137));
    ALGEBRAIC(:,23) =  ALGEBRAIC(:,13).*(ALGEBRAIC(:,19) - CONSTANTS(:,137));
    ALGEBRAIC(:,28) =  CONSTANTS(:,37).*(ALGEBRAIC(:,19) - ALGEBRAIC(:,26));
    ALGEBRAIC(:,32) =  ALGEBRAIC(:,30).*(ALGEBRAIC(:,19) - CONSTANTS(:,144));
    ALGEBRAIC(:,33) =  ALGEBRAIC(:,31).*(ALGEBRAIC(:,19) - CONSTANTS(:,30));
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