%% Finds the root xr where f(xr) = 0, given:
% @(x) f(x) function handle
% a and b are initial trial points which have opposite signs
% delta is the numerical tolerance on the root (in Newtons) -- not zero
% Obtained from the pseudo-code in https://en.wikipedia.org/wiki/Brent%27s_method
% Retrieved April 24, 2019
function [root,err_flag,num_evals] = Brents_Method(f,a,b,delta)

    %% Setup
    num_evals = 0;
    orig_a = a;
    orig_b = b;
    fa = f(a); num_evals = num_evals + 1;
    fb = f(b); num_evals = num_evals + 1;

    if fa*fb >= 0
        err_flag = -160; % Error 160 - f(a) and f(b) have the same sign
        root = -160;
        return;
    end
    
    if abs(fa) < abs(fb) % Swap 
        temp = a;
        a = b;      fa = f(a); num_evals = num_evals + 1;
        b = temp;   fb = f(b); num_evals = num_evals + 1;
    end
    
    c = a; fc = f(c); num_evals = num_evals + 1;
    d = [];
    m_flag = 1;
    converge = 0;
    num_iters = 0;
    max_iters = 10000;
    while converge == 0 && num_iters < max_iters
        
        %% Select advanced root finding method step
        if fa ~= fc && fb ~= fc
            s = a*fb*fc/((fa-fb)*(fa-fc)) +...
                b*fa*fc/((fb-fa)*(fb-fc)) +...
                c*fa*fb/((fc-fa)*(fc-fb)); % Inverse quadratic interpolation
        else
            s = b - fb*(b-a)/(fb-fa); % Secant method
        end
        
        %% Assess flag conditions
        condition = [0 0 0 0 0];
        % Condition 1
        if (s > (3*a+b)/4 && s < b) || (s > b && s < (3*a+b)/4)
        else
            condition(1) = 1; % s is outside of the two bounds
        end
        
        % Conditionn 2
        if m_flag == 1 && abs(s-b) >= abs(b-c)/2
            condition(2) = 1;
        end
        
        % Condition 3
        if isempty(d)
        else
            if m_flag == 0 && abs(s-b) >= abs(c-d)/2
                condition(3) = 1;
            end
        end
        
        % Condition 4
        if m_flag == 1 && abs(b-c) < abs(delta)
            condition(4) = 1;
        end
        
        % Condition 5
        if isempty(d)
        else
            if m_flag == 0 && abs(c-d) < abs(delta)
                condition(5) = 1;
            end
        end
        
        if sum(condition) > 0
            s = (a+b)/2;
            m_flag = 1;
        else
            m_flag = 0;
        end
        
        %% End of iteration processing
        fs = f(s);  num_evals = num_evals + 1;
        d = c;
        c = b;
        if fa*fs < 0
            b = s; fb = f(b); num_evals = num_evals + 1;
        else
            a = s; fa = f(a); num_evals = num_evals + 1;
        end
        
        if abs(fa) < abs(fb) % Swap
            temp = a;
            a = b;      fa = f(a); num_evals = num_evals + 1;
            b = temp;   fb = f(b); num_evals = num_evals + 1;
        end
        
        %% Convergence check
        if fb == 0 || abs(b-a) < abs(delta)
            converge = 1;
        end
        num_iters = num_iters + 1;
        
    end
    
    %% Output
    if converge == 1
        err_flag = 0;
        root = b;
    else
        err_flag = -161;
        root = -161;
    end
end