function [xmin,vof]=brent(a,b,c,f,tol,param);
%BRENT solves one dimensional minimization by parabolic interpolation.  
%  This program searches for the argmin, xmin, of f(x) given a bracketing 
%  triplet abscissas a,b,c.
%
%  INPUTS:
%    (i) a = lower bound of bracketing interval  [a,c]
%   (ii) b = point in bracketing interval where f(a) > f(b), f(c) > f(b), 
%            c > b > a
%  (iii) c = upper bound of bracketing interval  [a,c]
%            where a, b, c can be obtained using MNBRAK.M
%   (iv) f = a string giving the name of the function file for  y=f(x) 
%    (v) tol= convergence criterion
%
%  OUTPUTS:
%    (i) xmin = arg min f(x)
%   (ii) vof  = f(xmin)
%
%  SOURCE: William H. Press et al. (1986), Numerical Recipes, New York: 
%            Cambridge Univ. Press, pp. 283-286.
%
if nargin==6;
  pstr=',param';
else;
  pstr='';
end;
gold=.38196601;
maxits=100;
zeps=1e-10;
am=min(a,c);                   %  am and bm must be in ascending order,
bm=max(a,c);                   %  though the input abscissas need not be.
xm=.5*(am+bm);                 %  Initializations ...
v=b;                           
w=v;
x=v;
tol1=tol*abs(x)+zeps;
tol2=2*tol1;
e=0;                           %  This will be the distance moved on the 
eval(['fx=',f,'(x',pstr,');']); % step before last
fv=fx;
fw=fx;
iter=1;
while abs(x-xm) > tol2-.5*(bm-am) & iter<=maxits;
  flag=0;
  if abs(e) > tol1;            %  Construct a trial parabolic fit.
    r=(x-w)*(fx-fv);
    q=(x-v)*(fx-fw);
    p=(x-v)*q-(x-w)*r;
    q=2*(q-r);
    if q>0; p=-p; end;
    q=abs(q);
    etemp=e;
    e=d;                       
    if abs(p)<abs(.5*q*etemp) & p>q*(am-x) & p<q*(bm-x);
      d=p/q;                   %  Take the parabolic step.
      u=x+d;
      if u-am<tol2 | bm-u<tol2;
        if xm-x>0; 
          d= abs(tol1);
        else; 
          d=-abs(tol1);
        end;
      end;
      flag=1;                  %  Skip over the golden section step below.
    end;
  end;
  if flag~=1;                  %  Golden section step, which we take into
    if x>=xm; e=am-x; else; e=bm-x; end;   % the larger of the 2 segments.
    d=gold*e;
  end;
  if abs(d)>=tol1;             %  d has been computed either from parabolic
    u=x+d;                     %  fit or from golden section.
  else;
    if d>0; u=x+abs(tol1); else; u=x-abs(tol1); end;
  end;
  eval(['fu=',f,'(u',pstr,');']); %  This is the one function evaluation per 
  if fu<=fx;                   %     iteration
    if u>=x; am=x; else; bm=x; end;
    v=w;                       %  Housekeeping ...
    fv=fw;
    w=x;
    fw=fx;
    x=u;
    fx=fu;
  else;
    if u<x; am=u; else; bm=u; end;
    if fu<=fw | w==x;
      v=w;
      fv=fw;
      w=u;
      fw=fu;
    elseif fu<=fv | v==x | v==w;
      v=u;
      fv=fu;
    end;
  end;
  xm=.5*(am+bm);
  tol1=tol*abs(x)+zeps;
  tol2=2*tol1;
  iter=iter+1;
end;         % of while loop
xmin=x;
vof=fx;
if iter>=maxits; disp('Iteration limit of 100 reached in BRENT'); end;
