% [a,b,c]=dil8(fname)
% Function to bracket a minimum for a scalar function of one real variable.
% Inputs:
%  fname = string containing the name of a function m-file
% Returns:
%  a = left end of bracketing interval
%  b = point interior to bracketing interval with low function value
%  c = right end of bracketing interval

function [a,b,c] = dil8(fname)

% Philip D. Loewen
% 30 December 1997

gold = (3.0+sqrt(5.0))/2.0;   % Golden ratio expansion factor
maxswaps = 10;                % Maximum number of dilations allowed

x0 = 0.0; y0 = eval([fname,'(x0)']);   % Starting interval has endpoints
x1 = 1.0; y1 = eval([fname,'(x1)']);   % x0=0 and x1=1.

% Relabel the endpoints according to which has the greater y-value.
if y0<y1,
  xlo = x0; xhi = x1;
else
  xlo = x1; xhi = x0;
end;
ylo = min([y0,y1]); yhi = max([y0,y1]);

% Do a sample dilation to get the ball rolling.
xnew = xhi + gold*(xlo-xhi);
ynew = eval([fname,'(xnew)']);

% disp(['dil8:  Table of intermediate results follows:']);
% disp(['   xlo        ylo        xhi        yhi        xnew       ynew']);
% prtmtx([xlo ylo xhi yhi xnew ynew]);

% Loop through these steps until the minimum is bracketed.
nswaps = 0;
while ynew <= min(ylo,yhi),
  nswaps = nswaps+1;
  if nswaps > maxswaps, break, end;
  xlo = xnew; ylo = ynew;
  xnew = xhi + gold*(xlo-xhi);
  ynew = eval([fname,'(xnew)']);
%   prtmtx([xlo ylo xhi yhi xnew ynew]);
end;

% Report the results in increasing order.
a = min([xhi,xnew]);
b = xlo;
c = max([xhi,xnew]);

if nswaps > maxswaps,
  a = min([xlo,xhi,xnew]);   c = max([xlo,xhi,xnew]);
  disp(['*** dil8:  Function ',fname,' appears to be monotonic on the interval']);
  disp(['*** [',num2str(a),',',num2str(c),']:  try ', ...
                 'rescaling the independent variable.']);
end;

% disp(['dilate:  Finished; returning values ', ...
%       'a=',num2str(a),', b=',num2str(b),', c=',num2str(c),'.']);
return;
