% Modified version of the MATLAB template file FitNumDifSR.m from "A First
% Course on Kinetics and Reaction Engineering" used in the solution of
% Activity 16.2 of "A First Course on Kinetics and Reaction Engineering."
%
function Activity_16_2(p_guess)
    % Known quantities and constants (in consistent units)
    V = 50; % mL
    VFR = 5; % mL per min
    CS0=[12.6 11.2 9.0 8.1 6.3 5.6 4.3 3.6 2.3 1.0]'; % mmol per L
    CS0 = CS0/1000; % mmol per mL
    CP=[1.01 0.98 0.92 0.90 0.83 0.79 0.71 0.65 0.52 0.29]'; % mmol per L
    CP = CP/1000; % mmol per mL

    % Define set variables and response variable
    x = CS0;
    y_hat = CP;
    
    % Function that evaluates the model equations
	function g = nlaeqns(u)
        % u(1) = nS
        % u(2) = nP
        % u(3) = nW
        Vmax = p(1);
        Km = p(2);
        nS0 = x_set(1)*VFR;
        CS = u(1)/VFR;
        rate = Vmax*CS/(Km + CS);
		g = [
			nS0 - u(1) - V*rate;
            -u(2) + V*rate;
            -u(3) + V*rate;
		];
    end % of internal function nlaeqns

    % Function that calculates the responses using the model solution
    function y = nlmodel(p_current,x)
        % Make the current parameters available to the model equations
        p = p_current;
        y = zeros(n_data,1);
        
        for i = 1:n_data
            % Make the set variables available to the model equations
            for j = 1:n_set
                x_set(j) = x(i,j);
            end
            
            % Guesses
            % The measured response is available here as y_hat(i)
            nS0 = x_set(1)*VFR;
            nP_guess = y_hat(i)*VFR;
            u_guess = [
                nS0 - nP_guess
                nP_guess
                nP_guess
            ];
        
            % Solve the set of nonlinear modle equations
            % u = fsolve(@nlaeqns, u_guess);
            
            % Once you are sure fsolve is obtaining converged solutions,
            % you can suppress the "equations solved" messages by
            % commenting out the line above this comment and uncommenting
            % the next two lines.
            options = optimoptions('fsolve','Display','off');
            u = fsolve(@nlaeqns, u_guess, options);
            
            % Calculate the response
            y(i) = u(2)/VFR;
        end
    end % of internal function nlmodel

    % All the necessary internal functions are now defined

    % get the number of data points and set variables
    n_data = size(x,1);
    n_set = size(x,2);
    
    % create global variables to make parameters and set variables
    % available to the internal function that evaluates the model equations
    p = zeros(length(p_guess),1);
    x_set = zeros(size(x,2),1);
    
    % Perform the regression analysis
    [pf,resid,J,sigma,mse] = nlinfit(x,y_hat,@nlmodel,p_guess);
    confint = nlparci(pf,resid,'covar',sigma);
    pf_u = (confint(:,2) - confint(:,1))/2.0;
    
    % Calculate the model-predicted responses and the correlation
    % coefficient
    y = zeros(n_data,1);
    y = nlmodel(pf,x);

    % Generate the graphical output
    n_set = size(x,2);
    if n_set == 1
        % Model plot
        % Use many data points for the model so the curve is smooth
        xx = linspace(min(x),max(x))';
        yy = nlmodel(pf,xx);
        figure;
        plot(x,y,'r',x,y_hat,'o')
        xlabel('X')
        ylabel('Y')
        title('Model Plot')
        legend('Model','Data')
    else
        % Parity plot
        figure;
        plot(y_hat,y_hat,'r',y_hat,y,'o')
        xlabel('Measured Response')
        ylabel('Predicted Response')
        title('Parity Plot')

        % Residuals plots for each set variable
        for n = 1:n_set
            figure;
            plot(x(:,n),zeros(length(y),1),'r',x(:,n),resid,'o')
            xlabel(['x',num2str(n)])
            ylabel('Residual')
            title(['Residuals Plot vs. x',num2str(n)])
        end
    end
    
    % Display the correlation coefficient, parameters and uncertainties
    format shortE;
    r_squared = 1.0 - (sum(resid.^2))/(sum((y_hat-mean(y_hat)).^2))
    display('Best Values for the Parameters:')
    pf
    display('95% Confidence Intervals for the Parameters:')
    pf_u
end % of Activity_16_2.m