function H2O2_kinetics_kmb()

time = [0;120;300;600;1200;1800;2400;3000;3600]; %time in seconds
data = [1;0.91;0.78;0.59;0.37;0.22;0.13;0.08;0.05]; %H2O2 conc in mol/L

%--------------------------------------------------------------------------
%FIRST-ORDER:
%--------------------------------------------------------------------------

%If this is a first-order reaction, this plot looks linear:
    log_of_data = log(data);

    figure(1)
    plot(time, log_of_data, 'b*','Markersize',8,'LineWidth',3)
    xlabel('time, in seconds')
    ylabel('log of concentration')
    title('Log of Concentration of H2O2')
    grid on

%Get an initial guess for the first-order rate constant k:
    coeffs1 = polyfit(time, log_of_data, 1);  %This gets us the best-fit
                                              %line.  The 1 indicates a
                                              %first-degree polynomial.
                                              %The variable named coeffs
                                              %will be a vector that
                                              %contains the slope and
                                              %y-intercept of the
                                              %best-fit line.
    slope_of_line = coeffs1(1);
    yint_of_line = coeffs1(2);
    
    initial_k_1st_order = -slope_of_line; %initial k is the negative of the 
                                          %slope of the best-fit line FOR
                                          %FIRST-ORDER model. See
                                          %Equation (2).
    
%Now look for the best parameters y0 and k:
    %We need a starting point for our search:
        init_params_1st_order = [data(1),initial_k_1st_order]; 
                             %= [y(0), k]
                                        
    %Now the big event!  The next command starts with our initial guess at
    %the parameters for the model, and returns the parameters that produce
    %the best fit for the data (best fit here means that it minimizes the
    %sum of squared errors).
        bestfit_params_1st_order ...
            = fminsearch(@(params) sum_of_squared_errors(params, ...
                                                      @first_order_model, ...
                                                      data, time), ...
                                                    init_params_1st_order);
                                                
    %bestfit_params_1st_order   %Uncomment this line if you want to
                                %print the parameters to the screen
                                
    %What does our best-fit model predict for the given times in the table?
        tvec = time(1):0.1:time(end);
        bestfit_model_1st_order_output ...
                    = first_order_model(bestfit_params_1st_order, tvec);
        
    %Plot:
        figure(2)
        plot(time, data, 'b*', tvec, bestfit_model_1st_order_output,'k', ...
                'Markersize',8,'LineWidth',3)
        xlabel('time, in seconds')
        ylabel('concentration, in mol/L')
        title('Concentration of H2O2')
        legend('data','best-fit first-order model')
        grid on

%--------------------------------------------------------------------------
%SECOND-ORDER:
%--------------------------------------------------------------------------

%If this is a second-order reaction, this plot looks linear:
    one_over_data = 1./data;
    
    figure(3)
    plot(time, one_over_data, 'b*','Markersize',8,'LineWidth',3)
    xlabel('time, in seconds')
    ylabel('1/concentration')
    title('1/Concentration of H2O2')
    grid on
    
%If we try for a linear fit for that plot:
    coeffs2 = polyfit(time, one_over_data, 1);
    slope_of_line = coeffs2(1);
    yint_of_line = coeffs2(2);    
    
    initial_k_2nd_order = slope_of_line; %initial k is equal to the slope 
                                         %of the best-fit line FOR
                                         %SECOND-ORDER model.  See
                                         %Equation (4).    
%Now look for the best parameters:
    %We need a starting point for our search:
        init_params_2nd_order = [data(1), initial_k_2nd_order]; 
                             %= [y(0), k]
                             
    %Now find the best-fit parameters:                                    
        bestfit_params_2nd_order ...
            = fminsearch(@(params) sum_of_squared_errors(params, ...
                                                     @second_order_model, ...
                                                     data, time), ...
                                                    init_params_2nd_order);
                                                
    %bestfit_params_2nd_order   %Uncomment this line if you want to
                                %print the parameters to the screen
                                
    %What does our best-fit model predict for the given times in the table?                            
        bestfit_model_2nd_order_output = ...
                    second_order_model(bestfit_params_2nd_order, tvec);
        
    %Plot:
        figure(4)
        plot(time, data, 'b*', tvec, bestfit_model_2nd_order_output,'k', ...
                'Markersize',8,'LineWidth',3)
        xlabel('time, in seconds')
        ylabel('concentration, in mol/L')
        title('Concentration of H2O2')
        legend('data','best-fit second-order model')
        grid on

%--------------------------------------------------------------------------
%COMPARE SSEs (smallest SSE is an indicator of best model)
%--------------------------------------------------------------------------
SSE_first_order_model = sum_of_squared_errors(bestfit_params_1st_order, ...
                                                @first_order_model, ...
                                                data, time)
SSE_second_order_model = sum_of_squared_errors(bestfit_params_2nd_order, ...
                                                @second_order_model, ...
                                                data, time)

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Functions:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = sum_of_squared_errors(params, model, data, time)
    error = (model(params,time) - data);
    squared_errors = error.^2;
    out = sum(squared_errors);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = first_order_model(params,time)
    y0 = params(1);
    k = params(2);
    out = y0 * exp(-k * time);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = second_order_model(params,time)
    y0 = params(1);
    k = params(2);
    out = y0 ./(y0*k*time + 1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%