function [a,b,c,fa,fb,fc]=mnbrak(a,b,f,param);
%MNBRAK    Bracket a minimum of a scalar-valued function.
%          [a,b,c,fa,fb,fc]=mnbrak(a,b,f,param) searches in the 
%          downhill direction (defined by the function f at 
%          initial points a and b, where f(a)>f(b), b>a) to find 
%          points a, b, and c that bracket the minimum of the 
%          function (i.e. f(a)>f(b), f(c)>f(b)). Also returned
%          are the values fa=f(a), fb=f(b), fc=f(c), c > b > a.
%          The input f is a string giving the name of the function
%          file for y=f(x,param) that will be evaluated at 
%          x=a, b, and c, where param is some parameter vector.
%          [a,b,c,fa,fb,fc]=mnbrak(a,b,f) assumes that the function
%          in `f'.m has only one argument.
%

%          Ellen R. McGrattan, 8-22-87

%          References:
%
%          [1]  William H. Press et al. (1986), Numerical Recipes, New York: 
%               Cambridge Univ. Press, pp. 280-282.
%
gold=1.618034;
glimit=100;
tiny=1e-15;
maxits=100;
if nargin==4;
  pstr=',param';
else;
  pstr='';
end;
eval(['fa=',f,'(a',pstr,');']);
eval(['fb=',f,'(b',pstr,');']);
if fb > fa;                         %  Switch roles of a and b so that we go
  dum=a;                            %  downhill in the direction from a to b.
  a=b;
  b=dum;
  dum=fb;
  fb=fa;
  fa=dum;
end;
c=b+gold*(b-a);                     %  First guess for c.
eval(['fc=',f,'(c',pstr,');']);
iter=1;
while fb >= fc & iter<=maxits;      %  Loop until we bracket or hit iteration
  flag=0;                           %  limit.
  r=(b-a)*(fb-fc);                  %  Compute U by parabolic extrapolation 
  q=(b-c)*(fb-fa);                  %  from a,b,c.  Tiny is used to prevent
  if abs(q-r)>0;                    %  division by 0.
    dum=sign(q-r)*max(abs(q-r),tiny); 
  else;
    dum=-tiny;
  end;
  u=b-((b-c)*q-(b-a)*r)/(2*dum);
  ulim=b+glimit*(c-b);              %  We won't go farther than this.  
  if (b-u)*(u-c)>0;                 %  Now to test various possibilities:
    eval(['fu=',f,'(u',pstr,');']); %  1. Parabolic u between b and c
    if fu < fc;                     
      a=b;                          %     got minimum between b and c
      fa=fb;
      b=u;
      fb=fu;
      flag=1;
    end;
    if fu > fb; 
      c=u;                          %     got minimum between a and u
      fc=fu;
      flag=1;
    end;
    if flag==0;
      u=c+gold*(c-b);               %     parabolic fit was no use; use
      eval(['fu=',f,'(u',pstr,');']); %   default magnification
    end;
  elseif (c-u)*(u-ulim)>0;          %  2. Parabolic fit between c and 
    eval(['fu=',f,'(u',pstr,');']); %     allowed limit
    if fu < fc;
      b=c;
      c=u;
      u=c+gold*(c-b);
      fb=fc;
      fc=fu;
    end;
  elseif (u-ulim)*(ulim-c)>0;       %  3. Limit parabolic u to max allowed
    u=ulim;        
    eval(['fu=',f,'(u',pstr,');']);
  else;                             %  4. Reject parabolic u; use 
    u=c+gold*(c-b);                 %     default magnification
    eval(['fu=',f,'(u',pstr,');']);
  end;
  if flag==0;                       %  Eliminate oldest point and continue
    a=b;
    b=c; 
    c=u;
    fa=fb;
    fb=fc;
    fc=fu;
  end;
  iter=iter+1;
end;
if iter>=maxits; disp('Iteration limit of 100 reached in MNBRAK'); end;
