%% pwmnbrak.m is my MATLAB version of the mnbrak program from Numerical Recipes
%% (see p393 Section 10.1). For the given inputted function name (ie string input
%% required which gives name of the m file that contains the function of interest) 
%% and initial location guesses it returns three points that bracket the minimum 
%% of the function and the also returns the function values at the three points.
%%

function [ax,bx,cx,fa,fb,fc]=pwmnbrak(func,axi,bxi)

% func = function name (.m file name) (string).
% axi = 1st initial point (scalar).
% bxi = 2nd initial point (scalar).
% ax = 1st final point (scalar).
% bx = 2nd final point (scalar).
% cx = 3rd final point (scalar).
% fa = 1st final point function value (scalar).
% fb = 2nd final point function value (scalar).
% fc = 3rd final point function value (scalar).

%Assign input values to output
ax=axi;
bx=bxi;

%Required constants
gold=1.618034;
glimit=100;
tiny=1e-20;

%Determine initial function values and swap roles of a and b if required
%so go downhill from a to b.
fa=feval(func,ax);
fb=feval(func,bx);
if fb>fa
 dum=ax;
 ax=bx;
 bx=dum;
 dum=fb;
 fb=fa;
 fa=dum;
end

%First guess for c
cx=bx+gold*(bx-ax);
fc=feval(func,cx);

while fb>=fc
 r=(bx-ax)*(fb-fc);
 q=(bx-cx)*(fb-fa);
 ui=max(abs(q-r),tiny);
 if (q-r)<0;
  ui=-ui;
 end
 u=bx-((bx-cx)*q-(bx-ax)*r)/(2*ui);
 ulim=bx+glimit*(cx-bx);

 if ((bx-u)*(u-cx))>0                   %Parabolic u is between b and c: try it
  fu=feval(func,u);                
  if fu<fc                              %Got a minimum between b and c     
   ax=bx;
   fa=fb;
   bx=u;
   fb=fu;
   return
  elseif fu>fb                          %Got a minimum between a and u
   cx=u;
   fc=fu;
   return
  end
  u=cx+gold*(cx-bx);                    %Parabolic fit was not use. Use default magnification
  fu=eval(func,u);
 elseif ((cx-u)*(u-ulim))>0             %Parabolic fit is between c and its allowed limit
  fu=feval(func,u);
  if fu<fc
   bx=cx;
   cx=u;
   u=cx+gold*(cx-bx);
   fb=fc;
   fc=fu;
   fu=feval(func,u);
  end
 elseif ((u-ulim)*(ulim-cx))>=0         %Limit parabolic u to maximum allowed value
  u=ulim;
  fu=feval(func,u);
 else                                   %Reject parabolic u, use default magnification
  u=cx+gold*(cx-bx);
  fu=feval(func,u);              
 end
 ax=bx;                                 %Eliminate oldest point and continue
 bx=cx;
 cx=u;
 fa=fb;
 fb=fc;
 fc=fu;
end




