function [xmin,vof]=golden(a,b,c,f,tol,param);
%GOLDEN section search in one dimension.  This program searches for a 
%  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. 280-282.
%
if nargin==6;
  pstr=',param';
else;
  pstr='';
end;
q=.61803399;                          %  Golden ratios
r=1-q;
maxits=500;
x0=a;                                 %  At any time, we keep track of 4
x3=c;                                 %  points: x0, x1, x2, x3.
if abs(c-b) > abs(b-a);               %  Make x0 to x1 the smaller segment
  x1=b;
  x2=b+r*(c-b);                       %  and fill in the new point to be tried
else;
  x2=b;
  x1=b-r*(b-a);
end;
eval(['f1=',f,'(x1',pstr,');']);      %  The initial function evaluations.
eval(['f2=',f,'(x2',pstr,');']);
iter=1;                               %  Loop until minimum found or maximum
while abs(x3-x0)>tol*(abs(x1)+abs(x2)) & iter<=maxits;          % iterations 
  if f2<f1;                           %  One possible outcome
    x0=x1;
    x1=x2;
    x2=q*x1+r*x2;
    f0=f1;
    f1=f2;
    eval(['f2=',f,'(x2',pstr,');']);  %  and a new function evaluation
  else;                               %  The other outcome
    x3=x2;
    x2=x1;
    x1=q*x2+r*x0;
    f3=f2;
    f2=f1;
    eval(['f1=',f,'(x1',pstr,');']);  %  and a new function evaluation
  end;
  iter=iter+1;
end;
if f1<f2;                             %  Output the best of the two current
  vof=f1;                             %  points
  xmin=x1;
else;
  vof=f2;
  xmin=x2;
end;
if iter>=maxits; 
  fprintf('Iteration limit of %g reached in GOLDEN\n',maxits); 
end;
