tic;
[voi, states, algebraic, constants, rates] = solveModel();
toc;

function [algebraicVariableCount] = getalgebraicVariableCount()
    % Used later when setting a global variable with the number of algebraic variables.
    % Note: This is not the "main method".
    algebraicVariableCount =110;
end
% There are a total of 46 enp.stimDuries in each of the rate and state variable arrays.
% There are a total of 116 enp.stimDuries in the constant variable array.
%

function [voi, states, algebraic, constants, rates] = solveModel()

    [legend_states, legend_algebraic, legend_voi, legend_constants] = createlegends();

    % Create algebraic of correct size
    global algebraicVariableCount;  algebraicVariableCount = getalgebraicVariableCount();

    % change amoiunt of total givben chemical species available (labels) by a
    % multiplication factor (Gmult)
    labels = {'q_AC_init in comp env', ...
        'q_cAMP_init in comp env', ...
        'q_R_B1_init in comp env', ...
        'q_Gs_init in comp env', ...
        'q_R_M2_init in comp env', ...
        'q_Gi_init in comp env'};
%     labels = {'q_R_B1_init in comp env', ...
%         'q_R_M2_init in comp env'};
    [igs,~,~] = find_indices(labels, cellstr(legend_constants), cellstr(legend_states), cellstr(legend_algebraic));
    Gmult = 1e3; %1e3;
    [INIT_states, constants] = initConsts(Gmult, igs);

    % Set timespan to solve over
    voi = [0,6.5];%linspace(0, 0.5,2000);
%     voi = [0 1e-3];

    % Set numerical accuracy options for ODE solver
    options = odeset('RelTol', 1e-07, 'AbsTol', 1e-07, 'MaxStep', 1e-4); % step: 1e-4

    % extra options for stimuli [ L_B1, L_M2 ] settings 
    p = struct;
%     p.stimTimes = [[1,20];[3,40]];
	p.stimThreshold = [-0.95, -0.95]; %[1e-3,1e-3];
    p.stimHolding = [1e-9,1e-9];
    p.stimMag = [1.888e1,1.888e1]; %[1e6,1e6];
    p.stimStart = [1,2];
%     p.stimEnd = [3,4]; % not used
%     p.stimOffset = 0.0; % seconds % not used
    p.tRamp = [1,1];
    p.stimDur = [1,1];
    if false
        p.stimMag = p.stimHolding(1);
    end
    p.freq = [1,1]; %[0.25,0.25]; % per second
    % Solve model with ODE solver
    [voi, states] = ode15s(@(voi, states)computerates(voi, states, constants,p), voi, INIT_states, options);

    % Compute algebraic variables
    [rates, algebraic] = computerates(voi, states, constants,p);
    algebraic = computealgebraic(algebraic, constants, states, voi,p);

    labels = {'q_L_B1_init in comp env', ... 
                'q_L_M2_init in comp env', ...
                'q_LR_B1Gs in comp env', ...
                'q_LR_M2Gi in comp env', ...
                'q_aGs_GTP in comp env', ...
                'q_aGi_GTP in comp env', ...
                'q_aGs_GTP_AC in comp cAMP', ...
                'q_AC in comp env', ...
                'q_cAMP in comp env', ...
                'q_ATP in comp env'};

    [i_con, i_st, i_alg] = find_indices(labels, cellstr(legend_constants), cellstr(legend_states), cellstr(legend_algebraic));
    
    if true
        plot_selected(i_st,voi,states,legend_voi,legend_states,sprintf('states\nGmult=%g', Gmult),1)
        plot_selected(i_alg,voi,algebraic,legend_voi,legend_algebraic,sprintf('algebraic\nGmult=%g', Gmult),ceil(sqrt(length(i_alg))))
    end
    
    % debugging for q_LR_B1
    labels = {'q_LR_B1 in comp env',...
        'q_LR_B1 in comp LRGbinding_B1AR',...
        'q_LR_B1 in comp GsProtein'};

    if false
        [~, i_st, i_alg] = find_indices(labels, cellstr(legend_constants), cellstr(legend_states), cellstr(legend_algebraic));
        plot_selected(i_st,voi,states,legend_voi,legend_states,sprintf('LR_B1\nGmult=%g', Gmult),ceil(sqrt(length(i_st))))
        plot_selected(i_alg,voi,algebraic,legend_voi,legend_algebraic,sprintf('LR_B1\nGmult=%g', Gmult),ceil(sqrt(length(i_st))))
    end
    
    ids = [80,85,44,39,55];
    figure,
    for i = 1:2
        subplot(2,2,i)
        plot(voi(30:end), algebraic(30:end,ids(i)));
        xlabel('time (s)');
        legend(legend_algebraic(ids(i),:));
    end
    subplot(2,2,3)
    plot(voi(30:end), -algebraic(30:end,ids(1))+algebraic(30:end,ids(2)));
    xlabel('time (s)');
    title(legend_states(22,:));
    subplot(2,2,4)
    for j=1:3
        plot(voi(30:end), algebraic(30:end,ids(j+2))); hold on;
    end
    legend(legend_algebraic(ids(3),:),legend_algebraic(ids(4),:),legend_algebraic(ids(5),:));
    xlabel('time (s)');
    
    % debugging aGsGTP spike
    if true
        istart = 1; %55250;
        iend = length(voi); %57850;
        figure,
        subplot(2,2,1)
        plot(voi(istart:iend), states(istart:iend,14),'.'); % hold on;
        title('aGsGTP cAMP');
        subplot(2,2,2)
        istart2 = 1;
        iend2 = length(voi);
        plot(voi(istart2:iend2), states(istart2:iend2,36),'.'); 
        % legend('comp cAMP','comp GsProtein');
        title('aGsGTP GsProtein');
        subplot(2,2,3)
        plot(voi(istart:iend), rates(14,istart:iend),'.'); % hold on;
        title('d/dtaGsGTP cAMP');
        subplot(2,2,4)
        plot(voi(istart:iend), rates(36,istart:iend),'.'); 
        % legend('comp cAMP','comp GsProtein');
        title('d/dt aGsGTP GsProtein');
    end
end


function [] = plot_selected(ids,x,y,legend_x,legend_y,titlestr,ns)
    figure();
%     plot stimuli
    for i = 1:length(ids)
        subplot(ns,ns,i)
        plot(x(30:end), y(30:end,ids(i)));
        xlabel('time (s)');
        l = legend(legend_y(ids(i),:));
        set(l,'Interpreter','none');
    end
    suptitle(titlestr)    
end

function [i_con, i_st, i_alg] = find_indices(labels, legend_constants, legend_states, legend_algebraic)
% return the indices for the selected labels
    all_legends = [legend_constants; legend_states; legend_algebraic];
    
    i_con = [];
    for i = 1:length(labels)
        i_con = [i_con; find(strcmp(labels{i},legend_constants))];
    end
    i_st = [];
    for i = 1:length(labels)
        i_st = [i_st; find(strcmp(labels{i},legend_states))];
    end
    i_alg = [];
    for i = 1:length(labels)
        i_alg = [i_alg; find(strcmp(labels{i},legend_algebraic))];
    end
end

function [legend_states, legend_algebraic, legend_voi, legend_constants] = createlegends()
    legend_states = ''; legend_algebraic = ''; legend_voi = ''; legend_constants = '';
    legend_constants(:,1) = strpad('kappa_1a in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,2) = strpad('kappa_1b in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,3) = strpad('kappa_2a in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,4) = strpad('kappa_2b in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,5) = strpad('kappa_3a in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,6) = strpad('kappa_3b in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,7) = strpad('kappa_4a in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,8) = strpad('kappa_4b in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,9) = strpad('kappa_5 in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,10) = strpad('kappa_6 in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,11) = strpad('kappa_7 in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,12) = strpad('kappa_GiAC in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,13) = strpad('kappa_sig1_B1 in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,14) = strpad('kappa_sig2_B1 in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,15) = strpad('kappa_sig3_B1 in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,16) = strpad('kappa_sig4_B1 in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,17) = strpad('kappa_sig1_M2 in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,18) = strpad('kappa_sig2_M2 in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,19) = strpad('kappa_sig3_M2 in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,20) = strpad('kappa_sig4_M2 in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,21) = strpad('kappa_act1_Gs in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,22) = strpad('kappa_act2_Gs in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,23) = strpad('kappa_hyd_Gs in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,24) = strpad('kappa_reassoc_Gs in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,25) = strpad('kappa_act1_Gi in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,26) = strpad('kappa_act2_Gi in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,27) = strpad('kappa_hyd_Gi in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,28) = strpad('kappa_reassoc_Gi in comp BG_parameters (fmol_per_sec)');
    legend_constants(:,29) = strpad('K_ATP in comp BG_parameters (per_fmol)');
    legend_constants(:,30) = strpad('K_cAMP in comp BG_parameters (per_fmol)');
    legend_constants(:,31) = strpad('K_AC in comp BG_parameters (per_fmol)');
    legend_constants(:,32) = strpad('K_AC_ATP in comp BG_parameters (per_fmol)');
    legend_constants(:,33) = strpad('K_aGs_GTP_AC in comp BG_parameters (per_fmol)');
    legend_constants(:,34) = strpad('K_aGs_GTP_AC_ATP in comp BG_parameters (per_fmol)');
    legend_constants(:,35) = strpad('K_FSK_AC in comp BG_parameters (per_fmol)');
    legend_constants(:,36) = strpad('K_FSK_AC_ATP in comp BG_parameters (per_fmol)');
    legend_constants(:,37) = strpad('K_PDE in comp BG_parameters (per_fmol)');
    legend_constants(:,38) = strpad('K_PDE_cAMP in comp BG_parameters (per_fmol)');
    legend_constants(:,39) = strpad('K_five_AMP in comp BG_parameters (per_fmol)');
    legend_constants(:,40) = strpad('K_IBMX in comp BG_parameters (per_fmol)');
    legend_constants(:,41) = strpad('K_PDEinh in comp BG_parameters (per_fmol)');
    legend_constants(:,42) = strpad('K_aGs_GTP in comp BG_parameters (per_fmol)');
    legend_constants(:,43) = strpad('K_FSK in comp BG_parameters (per_fmol)');
    legend_constants(:,44) = strpad('K_aGi_GTP in comp BG_parameters (per_fmol)');
    legend_constants(:,45) = strpad('K_ACinh in comp BG_parameters (per_fmol)');
    legend_constants(:,46) = strpad('K_PPi in comp BG_parameters (per_fmol)');
    legend_constants(:,47) = strpad('K_L_B1 in comp BG_parameters (per_fmol)');
    legend_constants(:,48) = strpad('K_R_B1 in comp BG_parameters (per_fmol)');
    legend_constants(:,49) = strpad('K_Gs in comp BG_parameters (per_fmol)');
    legend_constants(:,50) = strpad('K_LR_B1 in comp BG_parameters (per_fmol)');
    legend_constants(:,51) = strpad('K_R_B1Gs in comp BG_parameters (per_fmol)');
    legend_constants(:,52) = strpad('K_LR_B1Gs in comp BG_parameters (per_fmol)');
    legend_constants(:,53) = strpad('K_L_M2 in comp BG_parameters (per_fmol)');
    legend_constants(:,54) = strpad('K_R_M2 in comp BG_parameters (per_fmol)');
    legend_constants(:,55) = strpad('K_Gi in comp BG_parameters (per_fmol)');
    legend_constants(:,56) = strpad('K_LR_M2 in comp BG_parameters (per_fmol)');
    legend_constants(:,57) = strpad('K_R_M2Gi in comp BG_parameters (per_fmol)');
    legend_constants(:,58) = strpad('K_LR_M2Gi in comp BG_parameters (per_fmol)');
    legend_constants(:,59) = strpad('K_beta_gamma_Gs in comp BG_parameters (per_fmol)');
    legend_constants(:,60) = strpad('K_aGs_GDP in comp BG_parameters (per_fmol)');
    legend_constants(:,61) = strpad('K_Pi in comp BG_parameters (per_fmol)');
    legend_constants(:,62) = strpad('K_beta_gamma_Gi in comp BG_parameters (per_fmol)');
    legend_constants(:,63) = strpad('K_aGi_GDP in comp BG_parameters (per_fmol)');
    legend_voi = strpad('time in comp env (second)');
    legend_constants(:,64) = strpad('vol_myo in comp env (pL)');
    legend_constants(:,65) = strpad('q_ATP_init in comp env');
    legend_constants(:,66) = strpad('q_AC_init in comp env');
    legend_constants(:,67) = strpad('q_cAMP_init in comp env');
    legend_constants(:,68) = strpad('q_AC_ATP_init in comp env');
    legend_constants(:,69) = strpad('q_FSK_init in comp env');
    legend_constants(:,70) = strpad('q_FSK_AC_init in comp env');
    legend_constants(:,71) = strpad('q_FSK_AC_ATP_init in comp env');
    legend_constants(:,72) = strpad('q_aGs_GTP_init in comp env');
    legend_constants(:,73) = strpad('q_aGs_GTP_AC_init in comp env');
    legend_constants(:,74) = strpad('q_aGs_GTP_AC_ATP_init in comp env');
    legend_constants(:,75) = strpad('q_PDE_init in comp env');
    legend_constants(:,76) = strpad('q_PDEinh_init in comp env');
    legend_constants(:,77) = strpad('q_PDE_cAMP_init in comp env');
    legend_constants(:,78) = strpad('q_IBMX_init in comp env');
    legend_constants(:,79) = strpad('q_five_AMP_init in comp env');
    legend_constants(:,80) = strpad('q_aGi_GTP_init in comp env');
    legend_constants(:,81) = strpad('q_ACinh_init in comp env');
    legend_constants(:,82) = strpad('q_PPi_init in comp env');
    legend_constants(:,83) = strpad('q_R_B1_init in comp env');
    legend_constants(:,84) = strpad('q_Gs_init in comp env');
    legend_algebraic(:,1) = strpad('q_L_B1_init in comp env');
    legend_constants(:,85) = strpad('q_LR_B1_init in comp env');
    legend_constants(:,86) = strpad('q_R_B1Gs_init in comp env');
    legend_constants(:,87) = strpad('q_LR_B1Gs_init in comp env');
    legend_algebraic(:,2) = strpad('q_L_M2_init in comp env');
    legend_constants(:,88) = strpad('q_R_M2_init in comp env');
    legend_constants(:,89) = strpad('q_Gi_init in comp env');
    legend_constants(:,90) = strpad('q_LR_M2_init in comp env');
    legend_constants(:,91) = strpad('q_R_M2Gi_init in comp env');
    legend_constants(:,92) = strpad('q_LR_M2Gi_init in comp env');
    legend_constants(:,93) = strpad('q_beta_gamma_Gs_init in comp env');
    legend_constants(:,94) = strpad('q_aGs_GDP_init in comp env');
    legend_constants(:,95) = strpad('q_Pi_init in comp env');
    legend_constants(:,96) = strpad('q_beta_gamma_Gi_init in comp env');
    legend_constants(:,97) = strpad('q_aGi_GDP_init in comp env');
    legend_algebraic(:,3) = strpad('q_ATP in comp env');
    legend_algebraic(:,4) = strpad('q_cAMP in comp env');
    legend_algebraic(:,5) = strpad('q_AC in comp env');
    legend_algebraic(:,6) = strpad('q_AC_ATP in comp env');
    legend_algebraic(:,7) = strpad('q_aGs_GTP_AC in comp env');
    legend_algebraic(:,8) = strpad('q_aGs_GTP_AC_ATP in comp env');
    legend_algebraic(:,9) = strpad('q_FSK_AC in comp env');
    legend_algebraic(:,10) = strpad('q_FSK_AC_ATP in comp env');
    legend_algebraic(:,11) = strpad('q_PDE in comp env');
    legend_algebraic(:,12) = strpad('q_PDE_cAMP in comp env');
    legend_algebraic(:,13) = strpad('q_five_AMP in comp env');
    legend_algebraic(:,14) = strpad('q_IBMX in comp env');
    legend_algebraic(:,15) = strpad('q_PDEinh in comp env');
    legend_algebraic(:,16) = strpad('q_aGs_GTP in comp env');
    legend_algebraic(:,17) = strpad('q_FSK in comp env');
    legend_algebraic(:,19) = strpad('q_aGi_GTP in comp env');
    legend_algebraic(:,21) = strpad('q_ACinh in comp env');
    legend_algebraic(:,24) = strpad('q_PPi in comp env');
    legend_algebraic(:,18) = strpad('q_L_B1 in comp env');
    legend_algebraic(:,20) = strpad('q_R_B1 in comp env');
    legend_algebraic(:,22) = strpad('q_Gs in comp env');
    legend_algebraic(:,26) = strpad('q_LR_B1 in comp env');
    legend_algebraic(:,27) = strpad('q_R_B1Gs in comp env');
    legend_algebraic(:,31) = strpad('q_LR_B1Gs in comp env');
    legend_algebraic(:,23) = strpad('q_L_M2 in comp env');
    legend_algebraic(:,25) = strpad('q_R_M2 in comp env');
    legend_algebraic(:,28) = strpad('q_Gi in comp env');
    legend_algebraic(:,32) = strpad('q_LR_M2 in comp env');
    legend_algebraic(:,36) = strpad('q_R_M2Gi in comp env');
    legend_algebraic(:,42) = strpad('q_LR_M2Gi in comp env');
    legend_algebraic(:,29) = strpad('q_beta_gamma_Gs in comp env');
    legend_algebraic(:,34) = strpad('q_aGs_GDP in comp env');
    legend_constants(:,106) = strpad('q_Pi in comp env');
    legend_algebraic(:,37) = strpad('q_beta_gamma_Gi in comp env');
    legend_algebraic(:,46) = strpad('q_aGi_GDP in comp env');
    legend_states(:,1) = strpad('q_ATP in comp cAMP');
    legend_states(:,2) = strpad('q_cAMP in comp cAMP');
    legend_states(:,3) = strpad('q_AC in comp cAMP');
    legend_states(:,4) = strpad('q_AC_ATP in comp cAMP');
    legend_states(:,5) = strpad('q_aGs_GTP_AC in comp cAMP');
    legend_states(:,6) = strpad('q_aGs_GTP_AC_ATP in comp cAMP');
    legend_states(:,7) = strpad('q_FSK_AC in comp cAMP');
    legend_states(:,8) = strpad('q_FSK_AC_ATP in comp cAMP');
    legend_states(:,9) = strpad('q_PDE in comp cAMP');
    legend_states(:,10) = strpad('q_PDE_cAMP in comp cAMP');
    legend_states(:,11) = strpad('q_five_AMP in comp cAMP');
    legend_states(:,12) = strpad('q_IBMX in comp cAMP');
    legend_states(:,13) = strpad('q_PDEinh in comp cAMP');
    legend_states(:,14) = strpad('q_aGs_GTP in comp cAMP');
    legend_states(:,15) = strpad('q_FSK in comp cAMP');
    legend_states(:,16) = strpad('q_aGi_GTP in comp cAMP');
    legend_states(:,17) = strpad('q_ACinh in comp cAMP');
    legend_states(:,18) = strpad('q_PPi in comp cAMP');
    legend_states(:,19) = strpad('q_L_B1 in comp LRGbinding_B1AR');
    legend_states(:,20) = strpad('q_R_B1 in comp LRGbinding_B1AR');
    legend_states(:,21) = strpad('q_Gs in comp LRGbinding_B1AR');
    legend_states(:,22) = strpad('q_LR_B1 in comp LRGbinding_B1AR');
    legend_states(:,23) = strpad('q_R_B1Gs in comp LRGbinding_B1AR');
    legend_states(:,24) = strpad('q_LR_B1Gs in comp LRGbinding_B1AR');
    legend_states(:,25) = strpad('q_L_M2 in comp LRGbinding_M2');
    legend_states(:,26) = strpad('q_R_M2 in comp LRGbinding_M2');
    legend_states(:,27) = strpad('q_Gi in comp LRGbinding_M2');
    legend_states(:,28) = strpad('q_LR_M2 in comp LRGbinding_M2');
    legend_states(:,29) = strpad('q_R_M2Gi in comp LRGbinding_M2');
    legend_states(:,30) = strpad('q_LR_M2Gi in comp LRGbinding_M2');
    legend_states(:,31) = strpad('q_R_B1 in comp GsProtein');
    legend_states(:,32) = strpad('q_Gs in comp GsProtein');
    legend_states(:,33) = strpad('q_R_B1Gs in comp GsProtein');
    legend_states(:,34) = strpad('q_LR_B1 in comp GsProtein');
    legend_states(:,35) = strpad('q_LR_B1Gs in comp GsProtein');
    legend_states(:,36) = strpad('q_aGs_GTP in comp GsProtein');
    legend_states(:,37) = strpad('q_beta_gamma_Gs in comp GsProtein');
    legend_states(:,38) = strpad('q_aGs_GDP in comp GsProtein');
    legend_constants(:,98) = strpad('q_Pi in comp GsProtein');
    legend_states(:,39) = strpad('q_R_M2 in comp GiProtein');
    legend_states(:,40) = strpad('q_Gi in comp GiProtein');
    legend_states(:,41) = strpad('q_R_M2Gi in comp GiProtein');
    legend_states(:,42) = strpad('q_LR_M2 in comp GiProtein');
    legend_states(:,43) = strpad('q_LR_M2Gi in comp GiProtein');
    legend_states(:,44) = strpad('q_aGi_GTP in comp GiProtein');
    legend_states(:,45) = strpad('q_beta_gamma_Gi in comp GiProtein');
    legend_states(:,46) = strpad('q_aGi_GDP in comp GiProtein');
    legend_constants(:,99) = strpad('q_Pi in comp GiProtein');
    legend_constants(:,100) = strpad('R in comp constants (J_per_K_per_mol)');
    legend_constants(:,101) = strpad('T in comp constants (kelvin)');
    legend_constants(:,102) = strpad('F in comp constants (C_per_mol)');
    legend_algebraic(:,99) = strpad('v1a in comp cAMP (fmol_per_sec)');
    legend_algebraic(:,100) = strpad('v1b in comp cAMP (fmol_per_sec)');
    legend_algebraic(:,101) = strpad('v2a in comp cAMP (fmol_per_sec)');
    legend_algebraic(:,102) = strpad('v2b in comp cAMP (fmol_per_sec)');
    legend_algebraic(:,103) = strpad('v3a in comp cAMP (fmol_per_sec)');
    legend_algebraic(:,104) = strpad('v3b in comp cAMP (fmol_per_sec)');
    legend_algebraic(:,105) = strpad('v4a in comp cAMP (fmol_per_sec)');
    legend_algebraic(:,107) = strpad('v4b in comp cAMP (fmol_per_sec)');
    legend_algebraic(:,109) = strpad('v5 in comp cAMP (fmol_per_sec)');
    legend_algebraic(:,106) = strpad('v6 in comp cAMP (fmol_per_sec)');
    legend_algebraic(:,108) = strpad('v7 in comp cAMP (fmol_per_sec)');
    legend_algebraic(:,110) = strpad('vGiAC in comp cAMP (fmol_per_sec)');
    legend_algebraic(:,30) = strpad('mu_ATP in comp cAMP (J_per_mol)');
    legend_algebraic(:,38) = strpad('mu_AC in comp cAMP (J_per_mol)');
    legend_algebraic(:,33) = strpad('mu_cAMP in comp cAMP (J_per_mol)');
    legend_algebraic(:,43) = strpad('mu_AC_ATP in comp cAMP (J_per_mol)');
    legend_algebraic(:,95) = strpad('mu_FSK in comp cAMP (J_per_mol)');
    legend_algebraic(:,59) = strpad('mu_FSK_AC in comp cAMP (J_per_mol)');
    legend_algebraic(:,64) = strpad('mu_FSK_AC_ATP in comp cAMP (J_per_mol)');
    legend_algebraic(:,92) = strpad('mu_aGs_GTP in comp cAMP (J_per_mol)');
    legend_algebraic(:,48) = strpad('mu_aGs_GTP_AC in comp cAMP (J_per_mol)');
    legend_algebraic(:,54) = strpad('mu_aGs_GTP_AC_ATP in comp cAMP (J_per_mol)');
    legend_algebraic(:,69) = strpad('mu_PDE in comp cAMP (J_per_mol)');
    legend_algebraic(:,89) = strpad('mu_PDEinh in comp cAMP (J_per_mol)');
    legend_algebraic(:,74) = strpad('mu_PDE_cAMP in comp cAMP (J_per_mol)');
    legend_algebraic(:,84) = strpad('mu_IBMX in comp cAMP (J_per_mol)');
    legend_algebraic(:,79) = strpad('mu_five_AMP in comp cAMP (J_per_mol)');
    legend_algebraic(:,96) = strpad('mu_aGi_GTP in comp cAMP (J_per_mol)');
    legend_algebraic(:,97) = strpad('mu_ACinh in comp cAMP (J_per_mol)');
    legend_algebraic(:,98) = strpad('mu_PPi in comp cAMP (J_per_mol)');
    legend_constants(:,103) = strpad('vol in comp cAMP (pL)');
    legend_algebraic(:,39) = strpad('mu_L_B1 in comp LRGbinding_B1AR (J_per_mol)');
    legend_algebraic(:,44) = strpad('mu_R_B1 in comp LRGbinding_B1AR (J_per_mol)');
    legend_algebraic(:,49) = strpad('mu_Gs in comp LRGbinding_B1AR (J_per_mol)');
    legend_algebraic(:,55) = strpad('mu_LR_B1 in comp LRGbinding_B1AR (J_per_mol)');
    legend_algebraic(:,60) = strpad('mu_R_B1Gs in comp LRGbinding_B1AR (J_per_mol)');
    legend_algebraic(:,65) = strpad('mu_LR_B1Gs in comp LRGbinding_B1AR (J_per_mol)');
    legend_algebraic(:,70) = strpad('vsig1_B1 in comp LRGbinding_B1AR (fmol_per_sec)');
    legend_algebraic(:,75) = strpad('vsig2_B1 in comp LRGbinding_B1AR (fmol_per_sec)');
    legend_algebraic(:,80) = strpad('vsig3_B1 in comp LRGbinding_B1AR (fmol_per_sec)');
    legend_algebraic(:,85) = strpad('vsig4_B1 in comp LRGbinding_B1AR (fmol_per_sec)');
    legend_constants(:,104) = strpad('vol in comp LRGbinding_B1AR (pL)');
    legend_algebraic(:,50) = strpad('mu_L_M2 in comp LRGbinding_M2 (J_per_mol)');
    legend_algebraic(:,56) = strpad('mu_R_M2 in comp LRGbinding_M2 (J_per_mol)');
    legend_algebraic(:,61) = strpad('mu_Gi in comp LRGbinding_M2 (J_per_mol)');
    legend_algebraic(:,66) = strpad('mu_LR_M2 in comp LRGbinding_M2 (J_per_mol)');
    legend_algebraic(:,71) = strpad('mu_R_M2Gi in comp LRGbinding_M2 (J_per_mol)');
    legend_algebraic(:,76) = strpad('mu_LR_M2Gi in comp LRGbinding_M2 (J_per_mol)');
    legend_algebraic(:,81) = strpad('vsig1_M2 in comp LRGbinding_M2 (fmol_per_sec)');
    legend_algebraic(:,86) = strpad('vsig2_M2 in comp LRGbinding_M2 (fmol_per_sec)');
    legend_algebraic(:,90) = strpad('vsig3_M2 in comp LRGbinding_M2 (fmol_per_sec)');
    legend_algebraic(:,93) = strpad('vsig4_M2 in comp LRGbinding_M2 (fmol_per_sec)');
    legend_constants(:,105) = strpad('vol in comp LRGbinding_M2 (pL)');
    legend_algebraic(:,72) = strpad('vact1_Gs in comp GsProtein (fmol_per_sec)');
    legend_algebraic(:,77) = strpad('vact2_Gs in comp GsProtein (fmol_per_sec)');
    legend_algebraic(:,82) = strpad('vhyd_Gs in comp GsProtein (fmol_per_sec)');
    legend_algebraic(:,87) = strpad('vreassoc_Gs in comp GsProtein (fmol_per_sec)');
    legend_algebraic(:,40) = strpad('mu_R_B1 in comp GsProtein (J_per_mol)');
    legend_algebraic(:,45) = strpad('mu_Gs in comp GsProtein (J_per_mol)');
    legend_algebraic(:,35) = strpad('mu_R_B1Gs in comp GsProtein (J_per_mol)');
    legend_algebraic(:,51) = strpad('mu_LR_B1 in comp GsProtein (J_per_mol)');
    legend_algebraic(:,41) = strpad('mu_LR_B1Gs in comp GsProtein (J_per_mol)');
    legend_algebraic(:,57) = strpad('mu_aGs_GTP in comp GsProtein (J_per_mol)');
    legend_algebraic(:,62) = strpad('mu_beta_gamma_Gs in comp GsProtein (J_per_mol)');
    legend_algebraic(:,67) = strpad('mu_aGs_GDP in comp GsProtein (J_per_mol)');
    legend_constants(:,107) = strpad('mu_Pi in comp GsProtein (J_per_mol)');
    legend_algebraic(:,83) = strpad('vact1_Gi in comp GiProtein (fmol_per_sec)');
    legend_algebraic(:,88) = strpad('vact2_Gi in comp GiProtein (fmol_per_sec)');
    legend_algebraic(:,91) = strpad('vhyd_Gi in comp GiProtein (fmol_per_sec)');
    legend_algebraic(:,94) = strpad('vreassoc_Gi in comp GiProtein (fmol_per_sec)');
    legend_algebraic(:,52) = strpad('mu_R_M2 in comp GiProtein (J_per_mol)');
    legend_algebraic(:,58) = strpad('mu_Gi in comp GiProtein (J_per_mol)');
    legend_algebraic(:,47) = strpad('mu_R_M2Gi in comp GiProtein (J_per_mol)');
    legend_algebraic(:,63) = strpad('mu_LR_M2 in comp GiProtein (J_per_mol)');
    legend_algebraic(:,53) = strpad('mu_LR_M2Gi in comp GiProtein (J_per_mol)');
    legend_algebraic(:,68) = strpad('mu_aGi_GTP in comp GiProtein (J_per_mol)');
    legend_algebraic(:,73) = strpad('mu_beta_gamma_Gi in comp GiProtein (J_per_mol)');
    legend_algebraic(:,78) = strpad('mu_aGi_GDP in comp GiProtein (J_per_mol)');
    legend_constants(:,108) = strpad('mu_Pi in comp GiProtein (J_per_mol)');
    legend_rates(:,1) = strpad('d/dt q_ATP in comp cAMP');
    legend_rates(:,3) = strpad('d/dt q_AC in comp cAMP');
    legend_rates(:,4) = strpad('d/dt q_AC_ATP in comp cAMP');
    legend_rates(:,2) = strpad('d/dt q_cAMP in comp cAMP');
    legend_rates(:,15) = strpad('d/dt q_FSK in comp cAMP');
    legend_rates(:,7) = strpad('d/dt q_FSK_AC in comp cAMP');
    legend_rates(:,8) = strpad('d/dt q_FSK_AC_ATP in comp cAMP');
    legend_rates(:,14) = strpad('d/dt q_aGs_GTP in comp cAMP');
    legend_rates(:,5) = strpad('d/dt q_aGs_GTP_AC in comp cAMP');
    legend_rates(:,6) = strpad('d/dt q_aGs_GTP_AC_ATP in comp cAMP');
    legend_rates(:,10) = strpad('d/dt q_PDE_cAMP in comp cAMP');
    legend_rates(:,9) = strpad('d/dt q_PDE in comp cAMP');
    legend_rates(:,12) = strpad('d/dt q_IBMX in comp cAMP');
    legend_rates(:,13) = strpad('d/dt q_PDEinh in comp cAMP');
    legend_rates(:,11) = strpad('d/dt q_five_AMP in comp cAMP');
    legend_rates(:,16) = strpad('d/dt q_aGi_GTP in comp cAMP');
    legend_rates(:,17) = strpad('d/dt q_ACinh in comp cAMP');
    legend_rates(:,18) = strpad('d/dt q_PPi in comp cAMP');
    legend_rates(:,19) = strpad('d/dt q_L_B1 in comp LRGbinding_B1AR');
    legend_rates(:,20) = strpad('d/dt q_R_B1 in comp LRGbinding_B1AR');
    legend_rates(:,21) = strpad('d/dt q_Gs in comp LRGbinding_B1AR');
    legend_rates(:,22) = strpad('d/dt q_LR_B1 in comp LRGbinding_B1AR');
    legend_rates(:,23) = strpad('d/dt q_R_B1Gs in comp LRGbinding_B1AR');
    legend_rates(:,24) = strpad('d/dt q_LR_B1Gs in comp LRGbinding_B1AR');
    legend_rates(:,25) = strpad('d/dt q_L_M2 in comp LRGbinding_M2');
    legend_rates(:,26) = strpad('d/dt q_R_M2 in comp LRGbinding_M2');
    legend_rates(:,27) = strpad('d/dt q_Gi in comp LRGbinding_M2');
    legend_rates(:,28) = strpad('d/dt q_LR_M2 in comp LRGbinding_M2');
    legend_rates(:,29) = strpad('d/dt q_R_M2Gi in comp LRGbinding_M2');
    legend_rates(:,30) = strpad('d/dt q_LR_M2Gi in comp LRGbinding_M2');
    legend_rates(:,31) = strpad('d/dt q_R_B1 in comp GsProtein');
    legend_rates(:,33) = strpad('d/dt q_R_B1Gs in comp GsProtein');
    legend_rates(:,32) = strpad('d/dt q_Gs in comp GsProtein');
    legend_rates(:,34) = strpad('d/dt q_LR_B1 in comp GsProtein');
    legend_rates(:,35) = strpad('d/dt q_LR_B1Gs in comp GsProtein');
    legend_rates(:,36) = strpad('d/dt q_aGs_GTP in comp GsProtein');
    legend_rates(:,37) = strpad('d/dt q_beta_gamma_Gs in comp GsProtein');
    legend_rates(:,38) = strpad('d/dt q_aGs_GDP in comp GsProtein');
    legend_rates(:,39) = strpad('d/dt q_R_M2 in comp GiProtein');
    legend_rates(:,41) = strpad('d/dt q_R_M2Gi in comp GiProtein');
    legend_rates(:,40) = strpad('d/dt q_Gi in comp GiProtein');
    legend_rates(:,42) = strpad('d/dt q_LR_M2 in comp GiProtein');
    legend_rates(:,43) = strpad('d/dt q_LR_M2Gi in comp GiProtein');
    legend_rates(:,44) = strpad('d/dt q_aGi_GTP in comp GiProtein');
    legend_rates(:,45) = strpad('d/dt q_beta_gamma_Gi in comp GiProtein');
    legend_rates(:,46) = strpad('d/dt q_aGi_GDP in comp GiProtein');
    legend_states  = legend_states';
    legend_algebraic = legend_algebraic';
    legend_rates = legend_rates';
    legend_constants = legend_constants';
end

function [states, constants] = initConsts(Gmult, igs)
    voi = 0; constants = []; states = []; algebraic = [];
    constants(:,1) = 3.92757e+06;
    constants(:,2) = 0.000820038;
    constants(:,3) = 310777;
    constants(:,4) = 0.0838605;
    constants(:,5) = 2.41805e+08;
    constants(:,6) = 2.81169e-17;
    constants(:,7) = 30855;
    constants(:,8) = 0.118673;
    constants(:,9) = 825.588;
    constants(:,10) = 5856.59;
    constants(:,11) = 455682;
    constants(:,12) = 1648.62;
    constants(:,13) = 1083.42;
    constants(:,14) = 16.7671;
    constants(:,15) = 3.64758;
    constants(:,16) = 3.34378e+06;
    constants(:,17) = 2855.54;
    constants(:,18) = 68.8777;
    constants(:,19) = 14.9839;
    constants(:,20) = 2.83927e+06;
    constants(:,21) = 2.19297;
    constants(:,22) = 0.00738313;
    constants(:,23) = 0.235482;
    constants(:,24) = 2.56289e+06;
    constants(:,25) = 0.903114;
    constants(:,26) = 9.47787e-05;
    constants(:,27) = 1.65719;
    constants(:,28) = 1.93732e+06;
    constants(:,29) = 5.4782e-05;
    constants(:,30) = 0.0159048;
    constants(:,31) = 3.65262;
    constants(:,32) = 7.08986;
    constants(:,33) = 4.9636;
    constants(:,34) = 2.94648;
    constants(:,35) = 0.063794;
    constants(:,36) = 0.103389;
    constants(:,37) = 1.72199;
    constants(:,38) = 1.22478;
    constants(:,39) = 0.0159048;
    constants(:,40) = 0.0198138;
    constants(:,41) = 35.211;
    constants(:,42) = 0.0987584;
    constants(:,43) = 1.15389e-05;
    constants(:,44) = 0.0140332;
    constants(:,45) = 17.6328;
    constants(:,46) = 0.000100127;
    constants(:,47) = 0.0569202;
    constants(:,48) = 0.00443997;
    constants(:,49) = 175.673;
    constants(:,50) = 1.31878;
    constants(:,51) = 885.438;
    constants(:,52) = 494.115;
    constants(:,53) = 0.0365205;
    constants(:,54) = 0.00814968;
    constants(:,55) = 36.3124;
    constants(:,56) = 1.55311;
    constants(:,57) = 335.945;
    constants(:,58) = 120.284;
    constants(:,59) = 2.55469e-05;
    constants(:,60) = 0.015617;
    constants(:,61) = 0.000229788;
    constants(:,62) = 0.000237839;
    constants(:,63) = 0.00221913;
    constants(:,64) = 34.4;
    constants(:,65) = 190;
    constants(:,66) = 0.0006143333; %*Gmult;        % q_AC;
    constants(:,67) = 1.889E-03; %*Gmult; % q_camp_init
    constants(:,68) = 1e-15;
    constants(:,69) = 3.800E-05;
    constants(:,70) = 1e-15;
    constants(:,71) = 1e-15;
    constants(:,72) = 0.0009519;
    constants(:,73) = 1e-15;
    constants(:,74) = 1e-15;
    constants(:,75) = 1.434E-03;
    constants(:,76) = 1e-15;
    constants(:,77) = 1e-15;
    constants(:,78) = 3.80E-02;
    constants(:,79) = 1e-15;
    constants(:,80) = 1e-15;
    constants(:,81) = 1e-15;
    constants(:,82) = 0.038;
    constants(:,83) = 0.0004579000; %*Gmult; % q_R_B1_init
    constants(:,84) = 0.1455400000; %*Gmult; % q_Gs_init
    constants(:,85) = 1e-1; % q_LR_B1_init % init, to get a positive mu
    constants(:,86) = 1e-15;
    constants(:,87) = 1e-15;
    constants(:,88) = 0.00072; %*Gmult; % q_R_M2_init
    constants(:,89) = 0.00836; %*Gmult; % q_Gi_init
    constants(:,90) = 1e-15; %q_LR_M2_init
    constants(:,91) = 1e-15;
    constants(:,92) = 1e-15;
    constants(:,93) = 1e-15;
    constants(:,94) = 1e-15;
    constants(:,95) = 570;
    constants(:,96) = 1e-15;
    constants(:,97) = 1e-15;
    states(:,1) = 1e-16;
    states(:,2) = 1e-16;
    states(:,3) = 1e-16;
    states(:,4) = 1e-16;
    states(:,5) = 1e-16;
    states(:,6) = 1e-16;
    states(:,7) = 1e-16;
    states(:,8) = 1e-16;
    states(:,9) = 1e-16;
    states(:,10) = 1e-16;
    states(:,11) = 1e-16;
    states(:,12) = 1e-16;
    states(:,13) = 1e-16;
    states(:,14) = 1e-16;
    states(:,15) = 1e-16;
    states(:,16) = 1e-16;
    states(:,17) = 1e-16;
    states(:,18) = 1e-16;
    states(:,19) = 1e-16;
    states(:,20) = 1e-16;
    states(:,21) = 1e-16;
    states(:,22) = 1e-16;
    states(:,23) = 1e-16;
    states(:,24) = 1e-16;
    states(:,25) = 1e-16;
    states(:,26) = 1e-16;
    states(:,27) = 1e-16;
    states(:,28) = 1e-16;
    states(:,29) = 1e-16;
    states(:,30) = 1e-16;
    states(:,31) = 1e-16;
    states(:,32) = 1e-16;
    states(:,33) = 1e-16;
    states(:,34) = 1e-16;
    states(:,35) = 1e-16;
    states(:,36) = 1e-16;
    states(:,37) = 1e-16;
    states(:,38) = 1e-16;
    constants(:,98) = 1e-16;
    states(:,39) = 1e-16;
    states(:,40) = 1e-16;
    states(:,41) = 1e-16;
    states(:,42) = 1e-16;
    states(:,43) = 1e-16;
    states(:,44) = 1e-16;
    states(:,45) = 1e-16;
    states(:,46) = 1e-16;
    constants(:,99) = 1e-16;
    constants(:,100) = 8.31;
    constants(:,101) = 310;
    constants(:,102) = 96485;
    constants(:,103) = 38.0;
    constants(:,104) = 34.4;
    constants(:,105) = 34.4;
    constants(:,106) = constants(:,95)+constants(:,98)+constants(:,99);
    constants(:,108) = 0.000000;
    constants(:,109) = 0.000000;
    constants(:,110) = 0.000000;
    constants(:,111) = 0.000000;
    constants(:,112) = 0.000000;
    constants(:,113) = 0.000000;
    constants(:,114) = 0.000000;
    constants(:,115) = 0.000000;
    constants(:,107) =  constants(:,100).*constants(:,101).*log( constants(:,61).*constants(:,106));
    constants(:,108) =  constants(:,100).*constants(:,101).*log( constants(:,61).*constants(:,106));
    if (isempty(states)), warning('Initial values for states not set');, end
    
    % multiply constant values by Gmult, with indices given in igs
    for j=1:length(igs)
        constants(:,igs(j)) = constants(:,igs(j))*Gmult;
    end
end

function [rates, algebraic] = computerates(voi, states, constants,p)
%     disp(voi);
    
    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(:,31) = constants(:,108);
    rates(:,33) = constants(:,109);
    rates(:,34) = constants(:,110);
    rates(:,35) = constants(:,111);
    rates(:,39) = constants(:,112);
    rates(:,41) = constants(:,113);
    rates(:,42) = constants(:,114);
    rates(:,43) = constants(:,115);
    algebraic(:,20) = constants(:,83)+states(:,20)+states(:,31);
    algebraic(:,44) =  constants(:,100).*constants(:,101).*log( constants(:,48).*algebraic(:,20));
    algebraic(:,22) = constants(:,84)+states(:,21)+states(:,32);
    algebraic(:,49) =  constants(:,100).*constants(:,101).*log( constants(:,49).*algebraic(:,22));
    algebraic(:,27) = constants(:,86)+states(:,23)+states(:,33);
    algebraic(:,60) =  constants(:,100).*constants(:,101).*log( constants(:,51).*algebraic(:,27));
    algebraic(:,70) =  constants(:,13).*exp((algebraic(:,44)+algebraic(:,49))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,60)./( constants(:,100).*constants(:,101)));
    if voi > 3
        j = 10;
    end

    % oldStim = piecewise({voi>p.stimStart(1)-1 & voi<=p.stimStart(1),  p.stimMag(1).*voi , voi>p.stimStart(1)&voi<p.stimEnd(1), p.stimMag(1)}, p.stimHolding(1)(1));
    m = p.stimMag(1)/p.tRamp(1);
    algebraic(:,1) = piecewise({voi<p.stimStart(1)&voi>p.stimStart(1) - p.tRamp(1), p.stimHolding(1)+ m.*((voi - p.stimStart(1))+p.tRamp(1)) , voi>=p.stimStart(1)&voi<p.stimStart(1)+p.stimDur(1), p.stimMag(1)+p.stimHolding(1) , voi<p.stimStart(1)+p.tRamp(1)+p.stimDur(1)&voi>=p.stimStart(1)+p.stimDur(1), p.stimHolding(1)+  - m.*(((voi - p.stimStart(1)) - p.tRamp(1)) - p.stimDur(1)) }, p.stimHolding(1));

    algebraic(:,18) = algebraic(:,1)+states(:,19);
    algebraic(:,39) =  constants(:,100).*constants(:,101).*log( constants(:,47).*algebraic(:,18));
    algebraic(:,31) = constants(:,87)+states(:,24)+states(:,35);
    algebraic(:,65) =  constants(:,100).*constants(:,101).*log( constants(:,52).*algebraic(:,31));
    algebraic(:,75) =  constants(:,14).*exp((algebraic(:,60)+algebraic(:,39))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,65)./( constants(:,100).*constants(:,101)));
    rates(:,23) = algebraic(:,70) - algebraic(:,75);
    algebraic(:,26) = constants(:,85)+states(:,22)+states(:,34);
    algebraic(:,55) =  constants(:,100).*constants(:,101).*log( constants(:,50).*algebraic(:,26));
    algebraic(:,80) =  constants(:,15).*exp((algebraic(:,55)+algebraic(:,49))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,65)./( constants(:,100).*constants(:,101)));
    rates(:,21) =  - algebraic(:,70) - algebraic(:,80);
    rates(:,24) = algebraic(:,75)+algebraic(:,80);
    algebraic(:,40) =  constants(:,100).*constants(:,101).*log( constants(:,48).*algebraic(:,20));
    algebraic(:,45) =  constants(:,100).*constants(:,101).*log( constants(:,49).*algebraic(:,22));
    algebraic(:,16) = constants(:,72)+states(:,14)+states(:,36);
    algebraic(:,57) =  constants(:,100).*constants(:,101).*log( constants(:,42).*algebraic(:,16));
    algebraic(:,29) = constants(:,93)+states(:,37);
    algebraic(:,62) =  constants(:,100).*constants(:,101).*log( constants(:,59).*algebraic(:,29));
    algebraic(:,72) =  constants(:,21).*(exp((algebraic(:,40)+algebraic(:,45)+constants(:,107))./( constants(:,100).*constants(:,101))) - exp((algebraic(:,57)+algebraic(:,62)+algebraic(:,40))./( constants(:,100).*constants(:,101))));
    algebraic(:,51) =  constants(:,100).*constants(:,101).*log( constants(:,50).*algebraic(:,26));
    algebraic(:,77) =  constants(:,22).*(exp((algebraic(:,51)+algebraic(:,45)+constants(:,107))./( constants(:,100).*constants(:,101))) - exp((algebraic(:,57)+algebraic(:,62)+algebraic(:,51))./( constants(:,100).*constants(:,101))));
    algebraic(:,34) = constants(:,94)+states(:,38);
    algebraic(:,67) =  constants(:,100).*constants(:,101).*log( constants(:,60).*algebraic(:,34));
    algebraic(:,82) =  constants(:,23).*(exp(algebraic(:,57)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,67)+constants(:,107))./( constants(:,100).*constants(:,101))));
    rates(:,36) = (algebraic(:,72)+algebraic(:,77)) - algebraic(:,82);
    algebraic(:,85) =  constants(:,16).*exp((algebraic(:,44)+algebraic(:,39))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,55)./( constants(:,100).*constants(:,101)));
    rates(:,19) =  - algebraic(:,75) - algebraic(:,85);
    rates(:,20) =  - algebraic(:,70) - algebraic(:,85);
    rates(:,22) =  - algebraic(:,80)+algebraic(:,85);
    algebraic(:,25) = constants(:,88)+states(:,26)+states(:,39);
    algebraic(:,56) =  constants(:,100).*constants(:,101).*log( constants(:,54).*algebraic(:,25));
    algebraic(:,28) = constants(:,89)+states(:,27)+states(:,40);
    algebraic(:,61) =  constants(:,100).*constants(:,101).*log( constants(:,55).*algebraic(:,28));
    algebraic(:,36) = constants(:,91)+states(:,29)+states(:,41);
    algebraic(:,71) =  constants(:,100).*constants(:,101).*log( constants(:,57).*algebraic(:,36));
    algebraic(:,81) =  constants(:,17).*exp((algebraic(:,56)+algebraic(:,61))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,71)./( constants(:,100).*constants(:,101)));
    % oldStim = piecewise({voi>p.stimStart(2)-1 & voi<=p.stimStart(2),  p.stimMag(2).*voi , voi>p.stimStart(2)&voi<p.stimEnd(2), p.stimMag(2)}, p.stimHolding(1)(2));
    m = p.stimMag(2)/p.tRamp(2);
    algebraic(:,2) = piecewise({voi<p.stimStart(2)&voi>p.stimStart(2) - p.tRamp(2), p.stimHolding(2)+ m.*((voi - p.stimStart(2))+p.tRamp(2)) , voi>=p.stimStart(2)&voi<p.stimStart(2)+p.stimDur(2), p.stimMag(2)+p.stimHolding(2) , voi<p.stimStart(2)+p.tRamp(2)+p.stimDur(2)&voi>=p.stimStart(2)+p.stimDur(2), p.stimHolding(2)+  - m.*(((voi - p.stimStart(2)) - p.tRamp(2)) - p.stimDur(2)) }, p.stimHolding(2));

    algebraic(:,23) = algebraic(:,2)+states(:,25);
    algebraic(:,50) =  constants(:,100).*constants(:,101).*log( constants(:,53).*algebraic(:,23));
    algebraic(:,42) = constants(:,92)+states(:,30)+states(:,43);
    algebraic(:,76) =  constants(:,100).*constants(:,101).*log( constants(:,58).*algebraic(:,42));
    algebraic(:,86) =  constants(:,18).*exp((algebraic(:,71)+algebraic(:,50))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,76)./( constants(:,100).*constants(:,101)));
    rates(:,29) = algebraic(:,81) - algebraic(:,86);
    algebraic(:,87) =  constants(:,24).*(exp((algebraic(:,67)+algebraic(:,62))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,45)./( constants(:,100).*constants(:,101))));
    rates(:,32) = ( - algebraic(:,72) - algebraic(:,77))+algebraic(:,87);
    rates(:,37) = (algebraic(:,72)+algebraic(:,77)) - algebraic(:,87);
    rates(:,38) = algebraic(:,82) - algebraic(:,87);
    algebraic(:,32) = constants(:,90)+states(:,28)+states(:,42);
    algebraic(:,66) =  constants(:,100).*constants(:,101).*log( constants(:,56).*algebraic(:,32));
    algebraic(:,90) =  constants(:,19).*exp((algebraic(:,66)+algebraic(:,61))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,76)./( constants(:,100).*constants(:,101)));
    rates(:,27) =  - algebraic(:,81) - algebraic(:,90);
    rates(:,30) = algebraic(:,86)+algebraic(:,90);
    algebraic(:,52) =  constants(:,100).*constants(:,101).*log( constants(:,54).*algebraic(:,25));
    algebraic(:,58) =  constants(:,100).*constants(:,101).*log( constants(:,55).*algebraic(:,28));
    algebraic(:,19) = constants(:,80)+states(:,16)+states(:,44);
    algebraic(:,68) =  constants(:,100).*constants(:,101).*log( constants(:,44).*algebraic(:,19));
    algebraic(:,37) = constants(:,96)+states(:,45);
    algebraic(:,73) =  constants(:,100).*constants(:,101).*log( constants(:,62).*algebraic(:,37));
    algebraic(:,83) =  constants(:,25).*(exp((algebraic(:,52)+algebraic(:,58)+constants(:,108))./( constants(:,100).*constants(:,101))) - exp((algebraic(:,68)+algebraic(:,73)+algebraic(:,52))./( constants(:,100).*constants(:,101))));
    algebraic(:,63) =  constants(:,100).*constants(:,101).*log( constants(:,56).*algebraic(:,32));
    algebraic(:,88) =  constants(:,26).*(exp((algebraic(:,63)+algebraic(:,58)+constants(:,108))./( constants(:,100).*constants(:,101))) - exp((algebraic(:,68)+algebraic(:,73)+algebraic(:,63))./( constants(:,100).*constants(:,101))));
    algebraic(:,46) = constants(:,97)+states(:,46);
    algebraic(:,78) =  constants(:,100).*constants(:,101).*log( constants(:,63).*algebraic(:,46));
    algebraic(:,91) =  constants(:,27).*(exp(algebraic(:,68)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,78)+constants(:,108))./( constants(:,100).*constants(:,101))));
    rates(:,44) = (algebraic(:,83)+algebraic(:,88)) - algebraic(:,91);
    algebraic(:,93) =  constants(:,20).*exp((algebraic(:,56)+algebraic(:,50))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,66)./( constants(:,100).*constants(:,101)));
    rates(:,25) =  - algebraic(:,86) - algebraic(:,93);
    rates(:,26) =  - algebraic(:,81) - algebraic(:,93);
    rates(:,28) =  - algebraic(:,90)+algebraic(:,93);
    algebraic(:,94) =  constants(:,28).*(exp((algebraic(:,78)+algebraic(:,73))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,58)./( constants(:,100).*constants(:,101))));
    rates(:,40) = ( - algebraic(:,83) - algebraic(:,88))+algebraic(:,94);
    rates(:,45) = (algebraic(:,83)+algebraic(:,88)) - algebraic(:,94);
    rates(:,46) = algebraic(:,91) - algebraic(:,94);
    algebraic(:,3) = constants(:,65)+states(:,1);
    algebraic(:,30) =  constants(:,100).*constants(:,101).*log( constants(:,29).*algebraic(:,3));
    algebraic(:,5) = constants(:,66)+states(:,3);
    algebraic(:,38) =  constants(:,100).*constants(:,101).*log( constants(:,31).*algebraic(:,5));
    algebraic(:,6) = constants(:,68)+states(:,4);
    algebraic(:,43) =  constants(:,100).*constants(:,101).*log( constants(:,32).*algebraic(:,6));
    algebraic(:,99) =  constants(:,1).*(exp((algebraic(:,38)+algebraic(:,30))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,43)./( constants(:,100).*constants(:,101))));
    algebraic(:,4) = constants(:,67)+states(:,2);
    algebraic(:,33) =  constants(:,100).*constants(:,101).*log( constants(:,30).*algebraic(:,4));
    algebraic(:,24) = constants(:,82)+states(:,18);
    algebraic(:,98) =  constants(:,100).*constants(:,101).*log( constants(:,46).*algebraic(:,24));
    algebraic(:,100) =  constants(:,2).*(exp(algebraic(:,43)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,38)+algebraic(:,33)+algebraic(:,98))./( constants(:,100).*constants(:,101))));
    rates(:,4) = algebraic(:,99) - algebraic(:,100);
    algebraic(:,7) = constants(:,73)+states(:,5);
    algebraic(:,48) =  constants(:,100).*constants(:,101).*log( constants(:,33).*algebraic(:,7));
    algebraic(:,8) = constants(:,74)+states(:,6);
    algebraic(:,54) =  constants(:,100).*constants(:,101).*log( constants(:,34).*algebraic(:,8));
    algebraic(:,101) =  constants(:,3).*(exp((algebraic(:,48)+algebraic(:,30))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,54)./( constants(:,100).*constants(:,101))));
    algebraic(:,102) =  constants(:,4).*(exp(algebraic(:,54)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,48)+algebraic(:,33)+algebraic(:,98))./( constants(:,100).*constants(:,101))));
    rates(:,6) = algebraic(:,101) - algebraic(:,102);
    algebraic(:,9) = constants(:,70)+states(:,7);
    algebraic(:,59) =  constants(:,100).*constants(:,101).*log( constants(:,35).*algebraic(:,9));
    algebraic(:,10) = constants(:,71)+states(:,8);
    algebraic(:,64) =  constants(:,100).*constants(:,101).*log( constants(:,36).*algebraic(:,10));
    algebraic(:,103) =  constants(:,5).*(exp((algebraic(:,59)+algebraic(:,30))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,64)./( constants(:,100).*constants(:,101))));
    rates(:,1) = ( - algebraic(:,99) - algebraic(:,103)) - algebraic(:,101);
    algebraic(:,104) =  constants(:,6).*(exp(algebraic(:,64)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,59)+algebraic(:,33)+algebraic(:,98))./( constants(:,100).*constants(:,101))));
    rates(:,8) = algebraic(:,103) - algebraic(:,104);
    algebraic(:,11) = constants(:,75)+states(:,9);
    algebraic(:,69) =  constants(:,100).*constants(:,101).*log( constants(:,37).*algebraic(:,11));
    algebraic(:,12) = constants(:,77)+states(:,10);
    algebraic(:,74) =  constants(:,100).*constants(:,101).*log( constants(:,38).*algebraic(:,12));
    algebraic(:,105) =  constants(:,7).*(exp((algebraic(:,69)+algebraic(:,33))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,74)./( constants(:,100).*constants(:,101))));
    rates(:,2) = (algebraic(:,100)+algebraic(:,104)+algebraic(:,102)) - algebraic(:,105);
    algebraic(:,92) =  constants(:,100).*constants(:,101).*log( constants(:,42).*algebraic(:,16));
    algebraic(:,106) =  constants(:,10).*(exp((algebraic(:,38)+algebraic(:,92))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,48)./( constants(:,100).*constants(:,101))));
    rates(:,14) =  - algebraic(:,106);
    rates(:,5) = (algebraic(:,106) - algebraic(:,101))+algebraic(:,102);
    algebraic(:,17) = constants(:,69)+states(:,15);
    algebraic(:,95) =  constants(:,100).*constants(:,101).*log( constants(:,43).*algebraic(:,17));
    algebraic(:,108) =  constants(:,11).*(exp((algebraic(:,95)+algebraic(:,38))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,59)./( constants(:,100).*constants(:,101))));
    rates(:,15) =  - algebraic(:,108);
    rates(:,7) = (algebraic(:,108)+algebraic(:,104)) - algebraic(:,103);
    algebraic(:,13) = constants(:,79)+states(:,11);
    algebraic(:,79) =  constants(:,100).*constants(:,101).*log( constants(:,39).*algebraic(:,13));
    algebraic(:,107) =  constants(:,8).*(exp(algebraic(:,74)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,69)+algebraic(:,79))./( constants(:,100).*constants(:,101))));
    rates(:,10) = algebraic(:,105) - algebraic(:,107);
    rates(:,11) = algebraic(:,107);
    algebraic(:,96) =  constants(:,100).*constants(:,101).*log( constants(:,44).*algebraic(:,19));
    algebraic(:,21) = constants(:,81)+states(:,17);
    algebraic(:,97) =  constants(:,100).*constants(:,101).*log( constants(:,45).*algebraic(:,21));
    algebraic(:,110) =  constants(:,12).*(exp((algebraic(:,38)+algebraic(:,96))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,97)./( constants(:,100).*constants(:,101))));
    rates(:,3) = (((algebraic(:,100) - algebraic(:,99)) - algebraic(:,106)) - algebraic(:,108)) - algebraic(:,110);
    algebraic(:,15) = constants(:,76)+states(:,13);
    algebraic(:,89) =  constants(:,100).*constants(:,101).*log( constants(:,41).*algebraic(:,15));
    algebraic(:,14) = constants(:,78)+states(:,12);
    algebraic(:,84) =  constants(:,100).*constants(:,101).*log( constants(:,40).*algebraic(:,14));
    algebraic(:,109) =  constants(:,9).*(exp((algebraic(:,69)+algebraic(:,84))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,89)./( constants(:,100).*constants(:,101))));
    rates(:,9) = (algebraic(:,107) - algebraic(:,105)) - algebraic(:,109);
    rates(:,12) =  - algebraic(:,109);
    rates(:,13) = algebraic(:,109);
    rates(:,16) =  - algebraic(:,110);
    rates(:,17) = algebraic(:,110);
    rates(:,18) = algebraic(:,110);
   rates = rates';
end

% Calculate algebraic variables
function algebraic = computealgebraic(algebraic, constants, states, voi, p)
    statesSize = size(states);
    statesColumnCount = statesSize(2);
    if ( statesColumnCount == 1)
        states = states';
        utilOnes = 1;
    else
        statesRowCount = statesSize(1);
        utilOnes = ones(statesRowCount, 1);
    end
    algebraic(:,20) = constants(:,83)+states(:,20)+states(:,31);
    algebraic(:,44) =  constants(:,100).*constants(:,101).*log( constants(:,48).*algebraic(:,20));
    algebraic(:,22) = constants(:,84)+states(:,21)+states(:,32);
    algebraic(:,49) =  constants(:,100).*constants(:,101).*log( constants(:,49).*algebraic(:,22));
    algebraic(:,27) = constants(:,86)+states(:,23)+states(:,33);
    algebraic(:,60) =  constants(:,100).*constants(:,101).*log( constants(:,51).*algebraic(:,27));
    algebraic(:,70) =  constants(:,13).*exp((algebraic(:,44)+algebraic(:,49))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,60)./( constants(:,100).*constants(:,101)));
    % oldStim = piecewise({voi>p.stimStart(1)-1 & voi<=p.stimStart(1),  p.stimMag(1).*voi , voi>p.stimStart(1)&voi<p.stimEnd(1), p.stimMag(1)}, p.stimHolding(1)(1));
    m = p.stimMag(1)/p.tRamp(1);
    algebraic(:,1) = piecewise({voi<p.stimStart(1)&voi>p.stimStart(1) - p.tRamp(1), p.stimHolding(1)+ m.*((voi - p.stimStart(1))+p.tRamp(1)) , voi>=p.stimStart(1)&voi<p.stimStart(1)+p.stimDur(1), p.stimMag(1)+p.stimHolding(1) , voi<p.stimStart(1)+p.tRamp(1)+p.stimDur(1)&voi>=p.stimStart(1)+p.stimDur(1), p.stimHolding(1)+  - m.*(((voi - p.stimStart(1)) - p.tRamp(1)) - p.stimDur(1)) }, p.stimHolding(1));

    algebraic(:,18) = algebraic(:,1)+states(:,19);
    algebraic(:,39) =  constants(:,100).*constants(:,101).*log( constants(:,47).*algebraic(:,18));
    algebraic(:,31) = constants(:,87)+states(:,24)+states(:,35);
    algebraic(:,65) =  constants(:,100).*constants(:,101).*log( constants(:,52).*algebraic(:,31));
    algebraic(:,75) =  constants(:,14).*exp((algebraic(:,60)+algebraic(:,39))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,65)./( constants(:,100).*constants(:,101)));
    algebraic(:,26) = constants(:,85)+states(:,22)+states(:,34);
    algebraic(:,55) =  constants(:,100).*constants(:,101).*log( constants(:,50).*algebraic(:,26));
    algebraic(:,80) =  constants(:,15).*exp((algebraic(:,55)+algebraic(:,49))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,65)./( constants(:,100).*constants(:,101)));
    algebraic(:,40) =  constants(:,100).*constants(:,101).*log( constants(:,48).*algebraic(:,20));
    algebraic(:,45) =  constants(:,100).*constants(:,101).*log( constants(:,49).*algebraic(:,22));
    algebraic(:,16) = constants(:,72)+states(:,14)+states(:,36);
    algebraic(:,57) =  constants(:,100).*constants(:,101).*log( constants(:,42).*algebraic(:,16));
    algebraic(:,29) = constants(:,93)+states(:,37);
    algebraic(:,62) =  constants(:,100).*constants(:,101).*log( constants(:,59).*algebraic(:,29));
    algebraic(:,72) =  constants(:,21).*(exp((algebraic(:,40)+algebraic(:,45)+constants(:,107))./( constants(:,100).*constants(:,101))) - exp((algebraic(:,57)+algebraic(:,62)+algebraic(:,40))./( constants(:,100).*constants(:,101))));
    algebraic(:,51) =  constants(:,100).*constants(:,101).*log( constants(:,50).*algebraic(:,26));
    algebraic(:,77) =  constants(:,22).*(exp((algebraic(:,51)+algebraic(:,45)+constants(:,107))./( constants(:,100).*constants(:,101))) - exp((algebraic(:,57)+algebraic(:,62)+algebraic(:,51))./( constants(:,100).*constants(:,101))));
    algebraic(:,34) = constants(:,94)+states(:,38);
    algebraic(:,67) =  constants(:,100).*constants(:,101).*log( constants(:,60).*algebraic(:,34));
    algebraic(:,82) =  constants(:,23).*(exp(algebraic(:,57)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,67)+constants(:,107))./( constants(:,100).*constants(:,101))));
    algebraic(:,85) =  constants(:,16).*exp((algebraic(:,44)+algebraic(:,39))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,55)./( constants(:,100).*constants(:,101)));
    algebraic(:,25) = constants(:,88)+states(:,26)+states(:,39);
    algebraic(:,56) =  constants(:,100).*constants(:,101).*log( constants(:,54).*algebraic(:,25));
    algebraic(:,28) = constants(:,89)+states(:,27)+states(:,40);
    algebraic(:,61) =  constants(:,100).*constants(:,101).*log( constants(:,55).*algebraic(:,28));
    algebraic(:,36) = constants(:,91)+states(:,29)+states(:,41);
    algebraic(:,71) =  constants(:,100).*constants(:,101).*log( constants(:,57).*algebraic(:,36));
    algebraic(:,81) =  constants(:,17).*exp((algebraic(:,56)+algebraic(:,61))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,71)./( constants(:,100).*constants(:,101)));
    % oldStim2 = piecewise({voi>p.stimStart(2)-1 & voi<=p.stimStart(2),  p.stimMag(2).*voi , voi>p.stimStart(2)&voi<p.stimEnd(2), p.stimMag(2)}, p.stimHolding(1)(2));
    m = p.stimMag(2)/p.tRamp(2);
    algebraic(:,2) = piecewise({voi<p.stimStart(2)&voi>p.stimStart(2) - p.tRamp(2), p.stimHolding(2)+ m.*((voi - p.stimStart(2))+p.tRamp(2)) , voi>=p.stimStart(2)&voi<p.stimStart(2)+p.stimDur(2), p.stimMag(2)+p.stimHolding(2) , voi<p.stimStart(2)+p.tRamp(2)+p.stimDur(2)&voi>=p.stimStart(2)+p.stimDur(2), p.stimHolding(2)+  - m.*(((voi - p.stimStart(2)) - p.tRamp(2)) - p.stimDur(2)) }, p.stimHolding(2));

    algebraic(:,23) = algebraic(:,2)+states(:,25);
    algebraic(:,50) =  constants(:,100).*constants(:,101).*log( constants(:,53).*algebraic(:,23));
    algebraic(:,42) = constants(:,92)+states(:,30)+states(:,43);
    algebraic(:,76) =  constants(:,100).*constants(:,101).*log( constants(:,58).*algebraic(:,42));
    algebraic(:,86) =  constants(:,18).*exp((algebraic(:,71)+algebraic(:,50))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,76)./( constants(:,100).*constants(:,101)));
    algebraic(:,87) =  constants(:,24).*(exp((algebraic(:,67)+algebraic(:,62))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,45)./( constants(:,100).*constants(:,101))));
    algebraic(:,32) = constants(:,90)+states(:,28)+states(:,42);
    algebraic(:,66) =  constants(:,100).*constants(:,101).*log( constants(:,56).*algebraic(:,32));
    algebraic(:,90) =  constants(:,19).*exp((algebraic(:,66)+algebraic(:,61))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,76)./( constants(:,100).*constants(:,101)));
    algebraic(:,52) =  constants(:,100).*constants(:,101).*log( constants(:,54).*algebraic(:,25));
    algebraic(:,58) =  constants(:,100).*constants(:,101).*log( constants(:,55).*algebraic(:,28));
    algebraic(:,19) = constants(:,80)+states(:,16)+states(:,44);
    algebraic(:,68) =  constants(:,100).*constants(:,101).*log( constants(:,44).*algebraic(:,19));
    algebraic(:,37) = constants(:,96)+states(:,45);
    algebraic(:,73) =  constants(:,100).*constants(:,101).*log( constants(:,62).*algebraic(:,37));
    algebraic(:,83) =  constants(:,25).*(exp((algebraic(:,52)+algebraic(:,58)+constants(:,108))./( constants(:,100).*constants(:,101))) - exp((algebraic(:,68)+algebraic(:,73)+algebraic(:,52))./( constants(:,100).*constants(:,101))));
    algebraic(:,63) =  constants(:,100).*constants(:,101).*log( constants(:,56).*algebraic(:,32));
    algebraic(:,88) =  constants(:,26).*(exp((algebraic(:,63)+algebraic(:,58)+constants(:,108))./( constants(:,100).*constants(:,101))) - exp((algebraic(:,68)+algebraic(:,73)+algebraic(:,63))./( constants(:,100).*constants(:,101))));
    algebraic(:,46) = constants(:,97)+states(:,46);
    algebraic(:,78) =  constants(:,100).*constants(:,101).*log( constants(:,63).*algebraic(:,46));
    algebraic(:,91) =  constants(:,27).*(exp(algebraic(:,68)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,78)+constants(:,108))./( constants(:,100).*constants(:,101))));
    algebraic(:,93) =  constants(:,20).*exp((algebraic(:,56)+algebraic(:,50))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,66)./( constants(:,100).*constants(:,101)));
    algebraic(:,94) =  constants(:,28).*(exp((algebraic(:,78)+algebraic(:,73))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,58)./( constants(:,100).*constants(:,101))));
    algebraic(:,3) = constants(:,65)+states(:,1);
    algebraic(:,30) =  constants(:,100).*constants(:,101).*log( constants(:,29).*algebraic(:,3));
    algebraic(:,5) = constants(:,66)+states(:,3);
    algebraic(:,38) =  constants(:,100).*constants(:,101).*log( constants(:,31).*algebraic(:,5));
    algebraic(:,6) = constants(:,68)+states(:,4);
    algebraic(:,43) =  constants(:,100).*constants(:,101).*log( constants(:,32).*algebraic(:,6));
    algebraic(:,99) =  constants(:,1).*(exp((algebraic(:,38)+algebraic(:,30))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,43)./( constants(:,100).*constants(:,101))));
    algebraic(:,4) = constants(:,67)+states(:,2);
    algebraic(:,33) =  constants(:,100).*constants(:,101).*log( constants(:,30).*algebraic(:,4));
    algebraic(:,24) = constants(:,82)+states(:,18);
    algebraic(:,98) =  constants(:,100).*constants(:,101).*log( constants(:,46).*algebraic(:,24));
    algebraic(:,100) =  constants(:,2).*(exp(algebraic(:,43)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,38)+algebraic(:,33)+algebraic(:,98))./( constants(:,100).*constants(:,101))));
    algebraic(:,7) = constants(:,73)+states(:,5);
    algebraic(:,48) =  constants(:,100).*constants(:,101).*log( constants(:,33).*algebraic(:,7));
    algebraic(:,8) = constants(:,74)+states(:,6);
    algebraic(:,54) =  constants(:,100).*constants(:,101).*log( constants(:,34).*algebraic(:,8));
    algebraic(:,101) =  constants(:,3).*(exp((algebraic(:,48)+algebraic(:,30))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,54)./( constants(:,100).*constants(:,101))));
    algebraic(:,102) =  constants(:,4).*(exp(algebraic(:,54)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,48)+algebraic(:,33)+algebraic(:,98))./( constants(:,100).*constants(:,101))));
    algebraic(:,9) = constants(:,70)+states(:,7);
    algebraic(:,59) =  constants(:,100).*constants(:,101).*log( constants(:,35).*algebraic(:,9));
    algebraic(:,10) = constants(:,71)+states(:,8);
    algebraic(:,64) =  constants(:,100).*constants(:,101).*log( constants(:,36).*algebraic(:,10));
    algebraic(:,103) =  constants(:,5).*(exp((algebraic(:,59)+algebraic(:,30))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,64)./( constants(:,100).*constants(:,101))));
    algebraic(:,104) =  constants(:,6).*(exp(algebraic(:,64)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,59)+algebraic(:,33)+algebraic(:,98))./( constants(:,100).*constants(:,101))));
    algebraic(:,11) = constants(:,75)+states(:,9);
    algebraic(:,69) =  constants(:,100).*constants(:,101).*log( constants(:,37).*algebraic(:,11));
    algebraic(:,12) = constants(:,77)+states(:,10);
    algebraic(:,74) =  constants(:,100).*constants(:,101).*log( constants(:,38).*algebraic(:,12));
    algebraic(:,105) =  constants(:,7).*(exp((algebraic(:,69)+algebraic(:,33))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,74)./( constants(:,100).*constants(:,101))));
    algebraic(:,92) =  constants(:,100).*constants(:,101).*log( constants(:,42).*algebraic(:,16));
    algebraic(:,106) =  constants(:,10).*(exp((algebraic(:,38)+algebraic(:,92))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,48)./( constants(:,100).*constants(:,101))));
    algebraic(:,17) = constants(:,69)+states(:,15);
    algebraic(:,95) =  constants(:,100).*constants(:,101).*log( constants(:,43).*algebraic(:,17));
    algebraic(:,108) =  constants(:,11).*(exp((algebraic(:,95)+algebraic(:,38))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,59)./( constants(:,100).*constants(:,101))));
    algebraic(:,13) = constants(:,79)+states(:,11);
    algebraic(:,79) =  constants(:,100).*constants(:,101).*log( constants(:,39).*algebraic(:,13));
    algebraic(:,107) =  constants(:,8).*(exp(algebraic(:,74)./( constants(:,100).*constants(:,101))) - exp((algebraic(:,69)+algebraic(:,79))./( constants(:,100).*constants(:,101))));
    algebraic(:,96) =  constants(:,100).*constants(:,101).*log( constants(:,44).*algebraic(:,19));
    algebraic(:,21) = constants(:,81)+states(:,17);
    algebraic(:,97) =  constants(:,100).*constants(:,101).*log( constants(:,45).*algebraic(:,21));
    algebraic(:,110) =  constants(:,12).*(exp((algebraic(:,38)+algebraic(:,96))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,97)./( constants(:,100).*constants(:,101))));
    algebraic(:,15) = constants(:,76)+states(:,13);
    algebraic(:,89) =  constants(:,100).*constants(:,101).*log( constants(:,41).*algebraic(:,15));
    algebraic(:,14) = constants(:,78)+states(:,12);
    algebraic(:,84) =  constants(:,100).*constants(:,101).*log( constants(:,40).*algebraic(:,14));
    algebraic(:,109) =  constants(:,9).*(exp((algebraic(:,69)+algebraic(:,84))./( constants(:,100).*constants(:,101))) - exp(algebraic(:,89)./( constants(:,100).*constants(:,101))));
    algebraic(:,35) =  constants(:,100).*constants(:,101).*log( constants(:,51).*algebraic(:,27));
    algebraic(:,41) =  constants(:,100).*constants(:,101).*log( constants(:,52).*algebraic(:,31));
    algebraic(:,47) =  constants(:,100).*constants(:,101).*log( constants(:,57).*algebraic(:,36));
    algebraic(:,53) =  constants(:,100).*constants(:,101).*log( constants(:,58).*algebraic(:,42));
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

