%
% *** function [ir,r] = secant(user_func,a,b,tol,itemax,iplot) ***
%
% Author: SJT
% Date: 1/19/2005
%
% Input
%   user_func = an inline function defined by the user
%   a         = one estimate of the root
%   b         = a second estimate of the root
%   tol       = desired tolerance
%   itemax    = the maximum number of iterations allowed
%   iplot     = 1 for plotting
%
% Output
%   ir = 0 if secant method converges
%      = 1 if maximum number of iterations exceeded
%      = 2 if divide by zero
%   r  = an estimate of the root
%
function [ir,r] = secant(user_func,a,b,tol,itemax,iplot)

% Check input arguments
if nargin < 6, iplot = 0; end
if nargin < 5, itemax = 8; end
if nargin < 4, tol = 1E-6; end

% Ensure that the two estimates are different
if abs(a-b) < tol
    fprintf('Error: Two estimates are the same \n')
    return
end

% Evaluate the function at a and b
fa = feval(user_func,a);
fb = feval(user_func,b);

% Ensure that the initial bracket spans a root (warning only)
if fa*fb > 0
    fprintf('WARNING: Initial range does not necessarily span a root \n')
end

ir = 0;
cnt = 0;
while abs(a-b) > tol | abs(fa) > tol
    cnt = cnt + 1;
    
    % Relabel a and b so that the (absolute) function value at a is the smaller
    if abs(fa) > abs(fb)
        tmp = a;
        a = b;
        b = tmp;
        tmp = fa;
        fa = fb;
        fb = tmp;
    end
    
    % Calculate slope of secant line
    slope = (fb-fa)/(b-a);  
    
    % Plot if desired
    if iplot == 1
        secant_plot(user_func,a,b,fa/slope)
    end
    
    % Old a        -> New b
    % New estimate -> a
    b  = a;
    fb = fa;
    if abs(slope) > 1E-13
        da = fa/slope;
    else
        fprintf('ERROR: Divide by zero in function secant \n')
        ir = 2;
        return
    end
    a  = a - da;
    fa = feval(user_func, a);
    
    fprintf('Iteration %3i, a = %12.6f, fa = %12.6f, b = %12.6f, fb = %12.6f \n',...
             cnt, a, fa, b, fb)
                  
    % Store convergence data
    xstep(cnt) = da;

    % Check for maximum number of iterations
    if cnt >= itemax
        fprintf('ERROR: Maximum number of iterations exceeded \n')
        ir = 1;
        break
    end
    
end

% Set best estimate of root
r = a;

fprintf('Convergence data for secant method (superlinear) \n')
fprintf('xstep(i)/xstep(i-1)  xstep(i)/xstep(i-1)^2 \n') 
for i = 1:cnt-1
    xratio(i)  = xstep(i+1)/xstep(i);
    xratio2(i) = xstep(i+1)/xstep(i)^2;
    fprintf('%16.5e  %18.5e \n', xratio(i), xratio2(i))
end