function N2O5_kinetics_kmb()

time = [0; 600; 1200; 1800; 2400; 3000; 3600; 4200; 4800; 6000]; 
        %time in seconds
data = [0.31;0.254;0.208;0.172;0.141;0.116;0.0964;0.0812;0.0669;0.0464]; 
        %N2O5 conc in mol/L

%--------------------------------------------------------------------------
%ZEROTH-ORDER:
%--------------------------------------------------------------------------

%If the reaction is zero-order, then this plot will look linear:
    figure(1)
    plot(time, data,'b*','Markersize',8,'LineWidth',3)
    xlabel('time, in seconds')
    ylabel('concentration, in mol/L')
    title('Concentration of N2O5 over Time')
    grid on
    
%If we think it is a zero-order reaction, we try to find k:
    coeffs0 = polyfit(time, 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 = coeffs0(1);
    yint_of_line = coeffs0(2);
    
    initial_k_zero_order = -slope_of_line; %initial k is the negative of the 
                                           %slope of the best-fit line FOR
                                            %ZERO-ORDER model. 
%Now look for the best parameters y0 and k:
    %We need a starting point for our search:
        init_params_zero_order = [data(1),initial_k_zero_order]; 
                             %= [y(0), k]
                             
    %Now find the best-fit parameters:                                    
        bestfit_params_zero_order ...
            = fminsearch(@(params) sum_of_squared_errors(params, ...
                                                      @zero_order_model, ...
                                                      data, time), ...
                                                    init_params_zero_order);
    %bestfit_params_zero_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_zero_order_output ...
                = zero_order_model(bestfit_params_zero_order, tvec);
        
    %Plot:
        figure(2)
        plot(time, data, 'b*', tvec, bestfit_model_zero_order_output,'k', ...
                'Markersize',8,'LineWidth',3)
        xlabel('time, in seconds')
        ylabel('concentration, in mol/L')
        title('Concentration of N2O5 over Time')
        legend('data','best-fit zero-order model')
        grid on
    

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

%If the reaction is first-order, then this plot will look linear:
    log_of_data = log(data);

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

%Get an initial guess for the 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:
    %We need a starting point for our search:
        init_params_1st_order = [data(1),initial_k_1st_order]; 
                             %= [y(0), k]
                                        
                                        
    %Now find the best-fit parameters:                                    
        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?                            
        bestfit_model_1st_order_output ...
                    = first_order_model(bestfit_params_1st_order, tvec);
                
    %Plot:
        figure(4)
        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 N2O5')
        legend('data','best-fit first-order model')
        grid on

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

%If this was a second-order reaction, this plot would be linear:
    one_over_data = 1./data;
    
    figure(5)
    plot(time, one_over_data, 'b*','Markersize',8,'LineWidth',3)
    xlabel('time, in seconds')
    ylabel('1/concentration')
    title('1/Concentration of N2O5')
    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(6)
        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 N2O5')
        legend('data','best-fit second-order model')
        grid on

%--------------------------------------------------------------------------
%COMPARE SSEs:
%--------------------------------------------------------------------------
SSE_zero_order_model = sum_of_squared_errors(bestfit_params_zero_order, ...
                                                @zero_order_model, ...
                                                data, time)
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 = zero_order_model(params, time)
    y0 = params(1);
    k = params(2);
    out = y0 - k*time;
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function out = sum_of_squared_errors(params, model, data, time)
    error = (model(params,time) - data);
    squared_errors = error.^2;
    out = sum(squared_errors);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%