% Modified version of the MATLAB template file FitNumDifMR.m from "A First
% Course on Kinetics and Reaction Engineering" used in the solution of
% Example 16.2 of "A First Course on Kinetics and Reaction Engineering."
%
function Example_16_2(p_guess)
    % Data from the problem statement (in consistent units)
    CG0 = [0.4809 0.4809 0.4809 0.4809 0.4725 0.4725 0.4725 0.4725]';
    t = [36.0 143.0 296.0 561.0 25.0 185.0 335.5 462.0]';
    CG = [0.4188 0.2882 0.2188 0.0868 0.4163 0.2629 0.1598 0.1074]';
    CH = [0.0063 0.0046 0.0027 0.0011 0.0067 0.0043 0.0021 0.0014]';
    CL = [0.0298 0.1200 0.1682 0.2578 0.0279 0.1334 0.2056 0.2326]';
    
    x = [CG0 t];
    y_hat = [CG CH CL];
    
	% Function that evaluates the ODEs
    function dudv = odeqns(v,u)
        % Current parameter values are available in column vector p
        % Current set variable values are available in column vector x_set
        dudv = [
            -p(1)*u(1)
			p(1)*u(1) - (p(2) + p(3))*u(2)
            p(2)*u(2)
        ];
    end % of internal function odeqns

    % Function that solves the model equations and calculates the model-
    % predicted responses
    function y = mrmodel(p_current)
        % Declare y
        y = zeros(n_data,n_resp);
        % 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
            
            % Initial and final values
            v0 = 0;
            u0 = [
                x_set(1)
                0
                0
            ];
            vf = x_set(2);
            
            % Solve the set of ordinary differential model equations
            [v,uu] = ode45(@odeqns,[v0 vf],u0);
% NOTE: For stiff equations replace ode45 with ode15s
            last_value = length(v);
            u = uu(last_value,:);
            
            % Use the solution to calculate the responses, y(1), y(2), ...
            % by evaluating f1(u(1),...,u(n),x_set(1),...,x_set(n_set)),
            % f2(u(1),...,u(n),x_set(1),...,x_set(n_set)), etc.
            y(i,:) = [
                u(1)
                u(2)
                u(3)
            ];
        end
    end % of internal function mrmodel

    % Function that evaluates the objective function
    function phi = objfcn(p_curr)
        % calculate the model responses using the current parameter values
        y = mrmodel(p_curr);
        
        % calculate the errors
        err = y_hat - y;
        
        % form the error matrix
        Emat = zeros(n_resp);
        for i = 1:n_resp
            % diagonal element
            for k = 1:n_data
                Emat(i,i) = Emat(i,i) + err(k,i)^2;
            end
            % below diagonal elements
            for j = i+1:n_resp
                for k = 1:n_data
                    Emat(j,i) = Emat(j,i) + err(k,j)*err(k,i);
                end
            end
            % above diagonal elements
            for j = i+1:n_resp
                Emat(i,j) = Emat(j,i);
            end
        end
        
        % evaluate the objective function
        phi = det(Emat);
    end % of internal function objfcn

    % All the necessary internal functions are now defined
    
    % get the number of data, parameters response and set variables
    n_data = size(y_hat,1);
    n_param = length(p_guess);
    n_resp = size(y_hat,2);
    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(n_param,1);
    x_set = zeros(n_set,1);

    % minimize the objective function and compute the responses
    y = zeros(n_data,n_resp);
    pf = fminsearch(@objfcn,p_guess);
    y = mrmodel(pf);
    resid = y - y_hat;
    
    % create a parity plot for each response variable
    for i = 1:n_resp
        figure;
        plot(y_hat(:,i),y_hat(:,i),'r',y_hat(:,i),y(:,i),'o')
        xlabel('Measured Response')
        ylabel('Predicted Response')
        title(['Parity Plot for Response variable ',num2str(i)])
        % create a residuals plot for this response variable versus each
        % set variable
        for n = 1:n_set
            figure;
            plot(x(:,n),zeros(size(y,1),1),'r',x(:,n),resid(:,i),'o')
            xlabel(['x',num2str(n)])
            ylabel('Residual')
            title(['y',num2str(i),' Residuals Plot vs. x',num2str(n)])
        end
    end

    % Display the parameters 
    display('Best Values for the Parameters:')
    pf
end % of Example_16_2.m