% Modified version of the MATLAB template file FitNumDifSR.m from "A First
% Course on Kinetics and Reaction Engineering" used in the solution of
% Example 16.1 of "A First Course on Kinetics and Reaction Engineering."
%
function Example_16_1(p_guess)
    % Known quantities and constants (in consistent units)
    Diam=0.25; % cm
    L=10.; % cm
    T=1500; % K
    P=1.0; % atm
    Rgas=82.06; % (atm-cm^3)/(mol-K)

    x = [10 20 30 40 60 75 100]'; % inlet volume flow (cm^3/min) 
    y_hat = [0.99 0.93 0.86 0.78 0.66 0.60 0.51]'; % COS conversion

	% Function that evaluates the ODEs
    function dnCOSdz = odeqns(z,nCOS)
        % Current parameter values are available in column vector p
        k = p(1);
        % Current set variable values are available in column vector x_set
        VFR0 = x_set(1);
        % Calculate quantities needed to evaluate the derivative
        n0COS = (P/Rgas/T)*VFR0;
        extent = n0COS - nCOS;
        ntot = n0COS + extent;
        PCOS = P*nCOS/ntot;
        % Evaluate the derivative
        dnCOSdz = [
            -pi()*Diam^2/4*k*PCOS
        ];
    end % of internal function odeqns
    
    % Function that solves the model equations and calculates the model-
    % predicted responses
    function y = nlmodel(p_current,x)
        % Declare y
        y = zeros(n_data,1);
        % Make the current parameters available to the model equations
        p = p_current;
        
        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
            VFR0 = x_set(1);
            
            % Initial and final values
            z0 = 0; 
            n0COS = (P/Rgas/T)*VFR0;
            zf = L;
            
            % Solve the set of ordinary differential model equations
            [z,uu] = ode45(@odeqns,[z0 zf],n0COS);
% NOTE: For stiff equations replace ode45 with ode15s
            last_value = length(z);
            nCOS = uu(last_value,:);
            
            % calculate the response
            y(i) = (n0COS-nCOS)/n0COS;
        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 and report the results
    [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
    y = zeros(n_data,1);
    y = nlmodel(pf,x);
    
    % Display the correlation coefficient, parameters and uncertainties
    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

    % 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

end % of Example_16_1.m