function approx = derivative(f, x, varargin)
%derivate computes a numerical derivative using the chosen method: forward
% difference, backward difference, centered difference, and Richardson
% extrapolation.
%
% Syntax:  approx = derivative(f, x, 'Method', 'forward', 'StepSize', 1e-4)
%
% Inputs:
% f             function handle 
% x             array of values to compute the numerical derivative
%
%   ---Optional Arguments---
% field name    name of target parameter 'Method', 'StepSize'
% field         parameter value 'forward', 'backward', 'centered',
%               'Richardson', for field name 'Method'. The default method
%               is forward difference.
%
%               floating point step size h for field name 'StepSize'. The
%               default step size depends on the difference method.
%
% Outputs:
% approx        the numerical approximation to the derivative of function
%
% Example: 
%   zero = derivative(@cos, pi/3)
%   zero = derivative(@sin, pi/2, 'Method', 'richardson')
%   zero = derivative(@(x) exp(x)-1, -1, 'stepsize', 1e-6)
%
% Other m-files required: none
% Subfunctions: 
% MAT-files required: 
%
% See also: 
%
% Author: Matlab Instructor
% 2022-Apr-04; Last revision: 2022-Apr-04
%
%------------- BEGIN CODE --------------

    % data structure
    data.f = f;
    data.x = x;
    data.h = sqrt(eps);
    method = @forwarddiff;
    
    % we don't know the order of the arguments
    stepSizeSet   = false;

    % varargin has some optional arguments!
    if nargin > 2
        
        % loop through the field names (odd indexes)
        for k=1:2:length(varargin)

            % compare the field name to constants
            switch(lower(varargin{k}))
                
                case char('method')
                    % use the even (following) element of varargin for the
                    % name of the method
                    
                    switch(lower(varargin{k+1}))
                        case char('backward')
                            method = @forwarddiff;
                            if ~stepSizeSet
                                data.h = -sqrt(eps);
                            end
                            if data.h > 0
                                data.h = -data.h;
                            end

                        case char('centered')
                            method = @centered;
                            if ~stepSizeSet
                                data.h = eps^(1/3);
                            end

                        case char('richardson')
                            method = @richardson;
                            if ~stepSizeSet
                                data.h = eps^(1/5);
                            end
                    end
                    
                case char('stepsize')
                    % use the next value in varargin for the floating point
                    % step size h
                    data.h      = varargin{k+1};
                    stepSizeSet = true;
            end
        end
    end
    
    approx = method(data);
end

function approx = forwarddiff(data)
%derivate computes a numerical derivative using forward difference, or
% backward if data.h < 0.
%
% Syntax:  approx = fowarddiff(data)
%
% Inputs:
% data  structure with the following fields:
%           f   function handle 
%           x   array of values to compute the numerical derivative
%           h   step size for difference methods
%
% Outputs:
% approx   the numerical approximation to the derivative of function
%
% Example: 
%   zero = forwarddiff(data)
%
% Other m-files required: none
% Subfunctions: 
% MAT-files required: 
%
% See also: 
%
% Author: Matlab Instructor
% 2022-Apr-04; Last revision: 2022-Apr-04
%
%------------- BEGIN CODE --------------

    % f'(x) = (f(x+h)-f(x))/h
    % if h<0, then the sign of the denominator will reverse the difference
    % in the numerator. 
    
    approx = (data.f(data.x + data.h) - data.f(data.x))./(data.h);

end

function approx = centered(data)
end

function approx = richardson(data)
end

%------------- END OF CODE --------------
%Please send suggestions for improvement of the above template header 
%to Denis Gilbert at this email address: gilbertd@dfo-mpo.gc.ca.
%Your contribution towards improving this template will be acknowledged in
%the "Changes" section of the TEMPLATE_HEADER web page on the Matlab
%Central File Exchange


