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 =0;
end
% There are a total of 0 entries in each of the rate and state variable arrays.
% There are a total of 156 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_CONSTANTS(:,1) = strpad('elementary_charge in component codata_2006_electromagnetic (coulomb)');
    LEGEND_CONSTANTS(:,2) = strpad('magnetic_flux_quantum in component codata_2006_electromagnetic (weber)');
    LEGEND_CONSTANTS(:,3) = strpad('conductance_quantum in component codata_2006_electromagnetic (siemens)');
    LEGEND_CONSTANTS(:,4) = strpad('inverse_of_conductance_quantum in component codata_2006_electromagnetic (ohm)');
    LEGEND_CONSTANTS(:,5) = strpad('Josephson_constant in component codata_2006_electromagnetic (hertz_per_volt)');
    LEGEND_CONSTANTS(:,6) = strpad('von_Klitzing_constant in component codata_2006_electromagnetic (ohm)');
    LEGEND_CONSTANTS(:,7) = strpad('Bohr_magneton in component codata_2006_electromagnetic (joule_per_tesla)');
    LEGEND_CONSTANTS(:,8) = strpad('nuclear_magneton in component codata_2006_electromagnetic (joule_per_tesla)');
    LEGEND_CONSTANTS(:,9) = strpad('speed_of_light_in_vacuum in component codata_2006_universal (metre_per_second)');
    LEGEND_CONSTANTS(:,10) = strpad('magnetic_constant in component codata_2006_universal (newton_per_square_ampere)');
    LEGEND_CONSTANTS(:,11) = strpad('electric_constant in component codata_2006_universal (farad_per_metre)');
    LEGEND_CONSTANTS(:,12) = strpad('characteristic_impedance_of_vacuum in component codata_2006_universal (ohm)');
    LEGEND_CONSTANTS(:,13) = strpad('Newtonian_constant_of_gravitation in component codata_2006_universal (cubic_metre_per_kilogram_per_square_second)');
    LEGEND_CONSTANTS(:,14) = strpad('Planck_constant in component codata_2006_universal (joule_second)');
    LEGEND_CONSTANTS(:,15) = strpad('Planck_constant_over_2_pi in component codata_2006_universal (joule_second)');
    LEGEND_CONSTANTS(:,16) = strpad('Planck_mass in component codata_2006_universal (kilogram)');
    LEGEND_CONSTANTS(:,17) = strpad('Planck_temperature in component codata_2006_universal (kelvin)');
    LEGEND_CONSTANTS(:,18) = strpad('Planck_length in component codata_2006_universal (metre)');
    LEGEND_CONSTANTS(:,19) = strpad('Planck_time in component codata_2006_universal (second)');
    LEGEND_CONSTANTS(:,20) = strpad('fine_structure_constant in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,21) = strpad('inverse_fine_structure_constant in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,22) = strpad('Rydberg_constant in component codata_2006_atomic_and_nuclear (per_metre)');
    LEGEND_CONSTANTS(:,23) = strpad('Bohr_radius in component codata_2006_atomic_and_nuclear (metre)');
    LEGEND_CONSTANTS(:,24) = strpad('Hartree_energy in component codata_2006_atomic_and_nuclear (joule)');
    LEGEND_CONSTANTS(:,25) = strpad('quantum_of_circulation in component codata_2006_atomic_and_nuclear (square_metre_per_second)');
    LEGEND_CONSTANTS(:,26) = strpad('Fermi_coupling_constant in component codata_2006_atomic_and_nuclear (per_square_gigaelectronvolt)');
    LEGEND_CONSTANTS(:,27) = strpad('weak_mixing_angle in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,28) = strpad('electron_mass in component codata_2006_atomic_and_nuclear (kilogram)');
    LEGEND_CONSTANTS(:,29) = strpad('electron_muon_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,30) = strpad('electron_tau_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,31) = strpad('electron_proton_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,32) = strpad('electron_neutron_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,33) = strpad('electron_deuteron_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,34) = strpad('electron_alpha_particle_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,35) = strpad('electron_charge_to_mass_quotient in component codata_2006_atomic_and_nuclear (coulomb_per_kilogram)');
    LEGEND_CONSTANTS(:,36) = strpad('electron_molar_mass in component codata_2006_atomic_and_nuclear (kilogram_per_mole)');
    LEGEND_CONSTANTS(:,37) = strpad('Compton_wavelength in component codata_2006_atomic_and_nuclear (metre)');
    LEGEND_CONSTANTS(:,38) = strpad('classic_electron_radius in component codata_2006_atomic_and_nuclear (metre)');
    LEGEND_CONSTANTS(:,39) = strpad('Thomson_cross_section in component codata_2006_atomic_and_nuclear (square_metre)');
    LEGEND_CONSTANTS(:,40) = strpad('electron_magnetic_moment in component codata_2006_atomic_and_nuclear (joule_per_tesla)');
    LEGEND_CONSTANTS(:,41) = strpad('electron_magnetic_moment_to_Bohr_magneton_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,42) = strpad('electron_magnetic_moment_to_nuclear_magneton_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,43) = strpad('electron_magnetic_moment_anomaly in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,44) = strpad('electron_g_factor in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,45) = strpad('electron_muon_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,46) = strpad('electron_proton_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,47) = strpad('electron_to_shielded_proton_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,48) = strpad('electron_neutron_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,49) = strpad('electron_deuteron_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,50) = strpad('electron_to_shielded_helion_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,51) = strpad('electron_gyromagnetic_ratio in component codata_2006_atomic_and_nuclear (per_second_per_tesla)');
    LEGEND_CONSTANTS(:,52) = strpad('muon_mass in component codata_2006_atomic_and_nuclear (kilogram)');
    LEGEND_CONSTANTS(:,53) = strpad('muon_electron_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,54) = strpad('muon_tau_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,55) = strpad('muon_proton_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,56) = strpad('muon_neutron_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,57) = strpad('muon_molar_mass in component codata_2006_atomic_and_nuclear (kilogram_per_mole)');
    LEGEND_CONSTANTS(:,58) = strpad('muon_Compton_wavelength in component codata_2006_atomic_and_nuclear (metre)');
    LEGEND_CONSTANTS(:,59) = strpad('muon_magnetic_moment in component codata_2006_atomic_and_nuclear (joule_per_tesla)');
    LEGEND_CONSTANTS(:,60) = strpad('muon_magnetic_moment_to_Bohr_magneton_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,61) = strpad('muon_magnetic_moment_to_nuclear_magneton_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,62) = strpad('muon_magnetic_moment_anomaly in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,63) = strpad('muon_g_factor in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,64) = strpad('muon_proton_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,65) = strpad('tau_mass in component codata_2006_atomic_and_nuclear (kilogram)');
    LEGEND_CONSTANTS(:,66) = strpad('tau_electron_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,67) = strpad('tau_muon_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,68) = strpad('tau_proton_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,69) = strpad('tau_neutron_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,70) = strpad('tau_molar_mass in component codata_2006_atomic_and_nuclear (kilogram_per_mole)');
    LEGEND_CONSTANTS(:,71) = strpad('tau_Compton_wavelength in component codata_2006_atomic_and_nuclear (metre)');
    LEGEND_CONSTANTS(:,72) = strpad('proton_mass in component codata_2006_atomic_and_nuclear (kilogram)');
    LEGEND_CONSTANTS(:,73) = strpad('proton_electron_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,74) = strpad('proton_muon_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,75) = strpad('proton_tau_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,76) = strpad('proton_neutron_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,77) = strpad('proton_charge_to_mass_quotient in component codata_2006_atomic_and_nuclear (coulomb_per_kilogram)');
    LEGEND_CONSTANTS(:,78) = strpad('proton_molar_mass in component codata_2006_atomic_and_nuclear (kilogram_per_mole)');
    LEGEND_CONSTANTS(:,79) = strpad('proton_Compton_wavelength in component codata_2006_atomic_and_nuclear (metre)');
    LEGEND_CONSTANTS(:,80) = strpad('proton_rms_electron_radius in component codata_2006_atomic_and_nuclear (metre)');
    LEGEND_CONSTANTS(:,81) = strpad('proton_magnetic_moment in component codata_2006_atomic_and_nuclear (joule_per_tesla)');
    LEGEND_CONSTANTS(:,82) = strpad('proton_magnetic_moment_to_Bohr_magneton_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,83) = strpad('proton_magnetic_moment_to_nuclear_magneton_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,84) = strpad('proton_g_factor in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,85) = strpad('proton_neutron_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,86) = strpad('shielded_proton_magnetic_moment in component codata_2006_atomic_and_nuclear (joule_per_tesla)');
    LEGEND_CONSTANTS(:,87) = strpad('shielded_proton_magnetic_moment_to_Bohr_magneton_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,88) = strpad('shielded_proton_magnetic_moment_to_nuclear_magneton_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,89) = strpad('proton_magnetic_shielding_correction in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,90) = strpad('proton_gyromagnetic_ratio in component codata_2006_atomic_and_nuclear (per_second_per_tesla)');
    LEGEND_CONSTANTS(:,91) = strpad('shielded_proton_gyromagnetic_ratio in component codata_2006_atomic_and_nuclear (per_second_per_tesla)');
    LEGEND_CONSTANTS(:,92) = strpad('neutron_mass in component codata_2006_atomic_and_nuclear (kilogram)');
    LEGEND_CONSTANTS(:,93) = strpad('neutron_electron_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,94) = strpad('neutron_muon_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,95) = strpad('neutron_tau_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,96) = strpad('neutron_proton_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,97) = strpad('neutron_molar_mass in component codata_2006_atomic_and_nuclear (kilogram_per_mole)');
    LEGEND_CONSTANTS(:,98) = strpad('neutron_Compton_wavelength in component codata_2006_atomic_and_nuclear (metre)');
    LEGEND_CONSTANTS(:,99) = strpad('neutron_magnetic_moment in component codata_2006_atomic_and_nuclear (joule_per_tesla)');
    LEGEND_CONSTANTS(:,100) = strpad('neutron_magnetic_moment_to_Bohr_magneton_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,101) = strpad('neutron_magnetic_moment_to_nuclear_magneton_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,102) = strpad('neutron_g_factor in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,103) = strpad('neutron_electron_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,104) = strpad('neutron_proton_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,105) = strpad('neutron_to_shielded_proton_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,106) = strpad('neutron_gyromagnetic_ratio in component codata_2006_atomic_and_nuclear (per_second_per_tesla)');
    LEGEND_CONSTANTS(:,107) = strpad('deuteron_mass in component codata_2006_atomic_and_nuclear (kilogram)');
    LEGEND_CONSTANTS(:,108) = strpad('deuteron_electron_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,109) = strpad('deuteron_proton_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,110) = strpad('deuteron_molar_mass in component codata_2006_atomic_and_nuclear (kilogram_per_mole)');
    LEGEND_CONSTANTS(:,111) = strpad('deuteron_rms_charge_radius in component codata_2006_atomic_and_nuclear (metre)');
    LEGEND_CONSTANTS(:,112) = strpad('deuteron_magnetic_moment in component codata_2006_atomic_and_nuclear (joule_per_tesla)');
    LEGEND_CONSTANTS(:,113) = strpad('deuteron_magnetic_moment_to_Bohr_magneton_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,114) = strpad('deuteron_magnetic_moment_to_nuclear_magneton_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,115) = strpad('deuteron_g_factor in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,116) = strpad('deuteron_electron_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,117) = strpad('deuteron_proton_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,118) = strpad('deuteron_neutron_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,119) = strpad('triton_mass in component codata_2006_atomic_and_nuclear (kilogram)');
    LEGEND_CONSTANTS(:,120) = strpad('triton_electron_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,121) = strpad('triton_proton_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,122) = strpad('triton_molar_mass in component codata_2006_atomic_and_nuclear (kilogram_per_mole)');
    LEGEND_CONSTANTS(:,123) = strpad('triton_magnetic_moment in component codata_2006_atomic_and_nuclear (joule_per_tesla)');
    LEGEND_CONSTANTS(:,124) = strpad('triton_magnetic_moment_to_Bohr_magneton_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,125) = strpad('triton_magnetic_moment_to_nuclear_magneton_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,126) = strpad('triton_g_factor in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,127) = strpad('triton_electron_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,128) = strpad('triton_proton_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,129) = strpad('triton_neutron_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,130) = strpad('helion_mass in component codata_2006_atomic_and_nuclear (kilogram)');
    LEGEND_CONSTANTS(:,131) = strpad('helion_electron_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,132) = strpad('helion_proton_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,133) = strpad('helion_molar_mass in component codata_2006_atomic_and_nuclear (kilogram_per_mole)');
    LEGEND_CONSTANTS(:,134) = strpad('shielded_helion_magnetic_moment in component codata_2006_atomic_and_nuclear (joule_per_tesla)');
    LEGEND_CONSTANTS(:,135) = strpad('shielded_helion_magnetic_moment_to_Bohr_magneton_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,136) = strpad('shielded_helion_magnetic_moment_to_nuclear_magneton_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,137) = strpad('shielded_helion_to_proton_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,138) = strpad('shielded_helion_to_shielded_proton_magnetic_moment_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,139) = strpad('shielded_helion_gyromagnetic_ratio in component codata_2006_atomic_and_nuclear (per_second_per_tesla)');
    LEGEND_CONSTANTS(:,140) = strpad('alpha_particle_mass in component codata_2006_atomic_and_nuclear (kilogram)');
    LEGEND_CONSTANTS(:,141) = strpad('alpha_particle_to_electron_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,142) = strpad('alpha_particle_to_proton_mass_ratio in component codata_2006_atomic_and_nuclear (dimensionless)');
    LEGEND_CONSTANTS(:,143) = strpad('alpha_particle_molar_mass in component codata_2006_atomic_and_nuclear (kilogram_per_mole)');
    LEGEND_CONSTANTS(:,144) = strpad('Avogadro_constant in component codata_2006_physicochemical (per_mole)');
    LEGEND_CONSTANTS(:,145) = strpad('atomic_mass_constant in component codata_2006_physicochemical (kilogram)');
    LEGEND_CONSTANTS(:,146) = strpad('Faraday_constant in component codata_2006_physicochemical (coulomb_per_mole)');
    LEGEND_CONSTANTS(:,147) = strpad('molar_Planck_constant in component codata_2006_physicochemical (joule_second_per_mole)');
    LEGEND_CONSTANTS(:,148) = strpad('molar_gas_constant in component codata_2006_physicochemical (joule_per_mole_per_kelvin)');
    LEGEND_CONSTANTS(:,149) = strpad('Boltzmann_constant in component codata_2006_physicochemical (joule_per_kelvin)');
    LEGEND_CONSTANTS(:,150) = strpad('molar_volume_of_ideal_gas in component codata_2006_physicochemical (cubic_metre_per_mole)');
    LEGEND_CONSTANTS(:,151) = strpad('Sackur_Tetrode_constant in component codata_2006_physicochemical (dimensionless)');
    LEGEND_CONSTANTS(:,152) = strpad('Stefan_Boltzmann_constant in component codata_2006_physicochemical (watt_per_square_meter_per_quartic_kelvin)');
    LEGEND_CONSTANTS(:,153) = strpad('first_radiation_constant in component codata_2006_physicochemical (watt_per_square_metre)');
    LEGEND_CONSTANTS(:,154) = strpad('first_radiation_constant_for_spectral_radiance in component codata_2006_physicochemical (watt_per_square_metre_per_steradian)');
    LEGEND_CONSTANTS(:,155) = strpad('second_radiation_constant in component codata_2006_physicochemical (metre_kelvin)');
    LEGEND_CONSTANTS(:,156) = strpad('Wien_displacement_law_constant in component codata_2006_physicochemical (metre_kelvin)');
    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) = 1.602176487e-19;
    CONSTANTS(:,2) = 2.067833667e-15;
    CONSTANTS(:,3) = 7.7480917004e-5;
    CONSTANTS(:,4) = 12906.4037787;
    CONSTANTS(:,5) = 483597.891e9;
    CONSTANTS(:,6) = 25812.807557;
    CONSTANTS(:,7) = 927.400915e-26;
    CONSTANTS(:,8) = 5.05078324e-27;
    CONSTANTS(:,9) = 299792458;
    CONSTANTS(:,10) = 1.2566370614359172953850573533118011536788677597500423283899778369231265625145e-6;
    CONSTANTS(:,11) = 8.8541878176203898505365630317107502606083701665994498081024171524053950954600e-12;
    CONSTANTS(:,12) = 376.73031346177065546819840042031930826862350835241865523207463829670726922131;
    CONSTANTS(:,13) = 6.67428e-11;
    CONSTANTS(:,14) = 6.62606896e-34;
    CONSTANTS(:,15) = 1.054571628e-34;
    CONSTANTS(:,16) = 6.62606896e-34;
    CONSTANTS(:,17) = 1.416785e32;
    CONSTANTS(:,18) = 1.616252e-35;
    CONSTANTS(:,19) = 5.39124e-44;
    CONSTANTS(:,20) = 7.2973525376e-3;
    CONSTANTS(:,21) = 137.035999679;
    CONSTANTS(:,22) = 10973731.568527;
    CONSTANTS(:,23) = 0.52917720859e-10;
    CONSTANTS(:,24) = 4.35974394e-18;
    CONSTANTS(:,25) = 3.6369475199e-4;
    CONSTANTS(:,26) = 1.6637e-5;
    CONSTANTS(:,27) = 0.22255;
    CONSTANTS(:,28) = 9.10938215e-31;
    CONSTANTS(:,29) = 4.83633171e-3;
    CONSTANTS(:,30) = 2.87564e-4;
    CONSTANTS(:,31) = 5.4461702177e-4;
    CONSTANTS(:,32) = 5.4386734459e-4;
    CONSTANTS(:,33) = 2.7244371093e-4;
    CONSTANTS(:,34) = 1.37093355570e-4;
    CONSTANTS(:,35) = -1.758820150e11;
    CONSTANTS(:,36) = 5.4857990943e-7;
    CONSTANTS(:,37) = 2.4263102175e-12;
    CONSTANTS(:,38) = 2.4263102175e-12;
    CONSTANTS(:,39) = 0.6652458558e-28;
    CONSTANTS(:,40) = -928.476377e-26;
    CONSTANTS(:,41) = -1.00115965218111;
    CONSTANTS(:,42) = -1838.28197092;
    CONSTANTS(:,43) = 1.15965218111e-3;
    CONSTANTS(:,44) = -2.0023193043622;
    CONSTANTS(:,45) = 206.7669877;
    CONSTANTS(:,46) = -658.2106848;
    CONSTANTS(:,47) = -658.2275971;
    CONSTANTS(:,48) = 960.92050;
    CONSTANTS(:,49) = -2143.923498;
    CONSTANTS(:,50) = 864.058257;
    CONSTANTS(:,51) = 1.760859770e11;
    CONSTANTS(:,52) = 1.88353130e-28;
    CONSTANTS(:,53) = 206.7682823;
    CONSTANTS(:,54) = 5.94592e-2;
    CONSTANTS(:,55) = 0.1126095261;
    CONSTANTS(:,56) = 0.1124545167;
    CONSTANTS(:,57) = 0.1134289256e-3;
    CONSTANTS(:,58) = 11.73444104e-15;
    CONSTANTS(:,59) = -4.49044786e-26;
    CONSTANTS(:,60) = -4.84197049e-3;
    CONSTANTS(:,61) = -8.89059705;
    CONSTANTS(:,62) = 1.16592069e-3;
    CONSTANTS(:,63) = -2.0023318414;
    CONSTANTS(:,64) = -3.183345137;
    CONSTANTS(:,65) = 3.16777e-27;
    CONSTANTS(:,66) = 3477.48;
    CONSTANTS(:,67) = 16.8183;
    CONSTANTS(:,68) = 1.89390;
    CONSTANTS(:,69) = 1.89129;
    CONSTANTS(:,70) = 1.90768e-3;
    CONSTANTS(:,71) = 0.69772e-15;
    CONSTANTS(:,72) = 1.672621637e-27;
    CONSTANTS(:,73) = 1836.15267247;
    CONSTANTS(:,74) = 8.88024339;
    CONSTANTS(:,75) = 0.528012;
    CONSTANTS(:,76) = 0.99862347824;
    CONSTANTS(:,77) = 9.57883392e7;
    CONSTANTS(:,78) = 1.00727646677e-3;
    CONSTANTS(:,79) = 1.3214098446e-15;
    CONSTANTS(:,80) = 0.8768e-15;
    CONSTANTS(:,81) = 1.410606662e-26;
    CONSTANTS(:,82) = 1.521032209e-3;
    CONSTANTS(:,83) = 2.792847356;
    CONSTANTS(:,84) = 5.585694713;
    CONSTANTS(:,85) = -1.45989806;
    CONSTANTS(:,86) = 1.410570419e-26;
    CONSTANTS(:,87) = 1.520993128e-3;
    CONSTANTS(:,88) = 2.792775598;
    CONSTANTS(:,89) = 25.694e-6;
    CONSTANTS(:,90) = 2.675222099e8;
    CONSTANTS(:,91) = 2.675153362e8;
    CONSTANTS(:,92) = 1.674927211e-27;
    CONSTANTS(:,93) = 1838.6836605;
    CONSTANTS(:,94) = 8.89248409;
    CONSTANTS(:,95) = 0.528740;
    CONSTANTS(:,96) = 1.00137841918;
    CONSTANTS(:,97) = 1.00866491597e-3;
    CONSTANTS(:,98) = 1.3195908951e-15;
    CONSTANTS(:,99) = -0.96623641e-26;
    CONSTANTS(:,100) = -1.04187563e-3;
    CONSTANTS(:,101) = -1.91304273;
    CONSTANTS(:,102) = -3.82608545;
    CONSTANTS(:,103) = 1.04066882e-3;
    CONSTANTS(:,104) = -0.68497934;
    CONSTANTS(:,105) = -0.68499694;
    CONSTANTS(:,106) = 1.83247185e8;
    CONSTANTS(:,107) = 3.34358320e-27;
    CONSTANTS(:,108) = 3670.4829654;
    CONSTANTS(:,109) = 1.99900750108;
    CONSTANTS(:,110) = 2.013553212724e-3;
    CONSTANTS(:,111) = 2.1402e-15;
    CONSTANTS(:,112) = 0.433073465e-26;
    CONSTANTS(:,113) = -0.4669754556e-3;
    CONSTANTS(:,114) = 0.8574382308;
    CONSTANTS(:,115) = 0.8574382308;
    CONSTANTS(:,116) = -4.664345537e-4;
    CONSTANTS(:,117) = 0.3070122070;
    CONSTANTS(:,118) = -0.44820652;
    CONSTANTS(:,119) = 5.00735588e-27;
    CONSTANTS(:,120) = 5496.9215269;
    CONSTANTS(:,121) = 2.9937170309;
    CONSTANTS(:,122) = 3.0155007134e-3;
    CONSTANTS(:,123) = 1.504609361e-26;
    CONSTANTS(:,124) = 1.622393657e-3;
    CONSTANTS(:,125) = 2.978962448;
    CONSTANTS(:,126) = 5.957924896;
    CONSTANTS(:,127) = -1.620514423e-3;
    CONSTANTS(:,128) = 1.066639908;
    CONSTANTS(:,129) = -1.55718553;
    CONSTANTS(:,130) = 5.00641192e-27;
    CONSTANTS(:,131) = 5495.8852765;
    CONSTANTS(:,132) = 2.9931526713;
    CONSTANTS(:,133) = 3.0149322473e-3;
    CONSTANTS(:,134) = -1.074552982e-26;
    CONSTANTS(:,135) = -1.158671471e-3;
    CONSTANTS(:,136) = -2.127497718;
    CONSTANTS(:,137) = -0.761766558;
    CONSTANTS(:,138) = -0.7617861313;
    CONSTANTS(:,139) = 2.037894730e8;
    CONSTANTS(:,140) = 6.64465620e-27;
    CONSTANTS(:,141) = 7294.2995365;
    CONSTANTS(:,142) = 3.97259968951;
    CONSTANTS(:,143) = 4.001506179127e-3;
    CONSTANTS(:,144) = 6.02214179e23;
    CONSTANTS(:,145) = 1.660538782e-27;
    CONSTANTS(:,146) = 96485.3399;
    CONSTANTS(:,147) = 3.9903126821e-10;
    CONSTANTS(:,148) = 8.314472;
    CONSTANTS(:,149) = 1.3806504e-23;
    CONSTANTS(:,150) = 22.413996e-3;
    CONSTANTS(:,151) = -1.1648677;
    CONSTANTS(:,152) = 5.670400e-8;
    CONSTANTS(:,153) = 3.74177118e-16;
    CONSTANTS(:,154) = 1.191042759e-16;
    CONSTANTS(:,155) = 1.4387752e-2;
    CONSTANTS(:,156) = 2.8977685e-3;
    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 = 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
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