function root = root_bisect(lower_bound, upper_bound, epsilon, func)
%
% Find the root of the given function within the given interval,
% using the bisection method.
%
% Arguments:
%
%      lower_bound     (input)    lower bound of interval
%                                     on which to search for root
%
%      upper_bound     (input)    upper bound of interval
%
%      epsilon         (input)    when fractional change in root
%                                     reaches this limit, stop
%
%      func            (input)    name of function whose root
%                                     we seek -- should take
%                                     a single argument
%
%      root           (output)    argument at which the given
%                                     function has value zero
%
% MWR 3/25/2002
%

% sanity checks
if (lower_bound >= upper_bound)
   error('lower_bound must be less than upper_bound');
end
if (eps <= 0)
   error('epsilon must be greater than zero');
end

% if there is no guaranteed root between the given bounds, 
% then quit with an error message
y_low = feval(func, lower_bound);
y_high = feval(func, upper_bound);
if (sign(y_low) == sign(y_high))
   error('function has same sign at both bounds -- quitting');
end

done = 0;
iterations = 0;
old_mid = lower_bound;
while (done == 0) 
   y_low = feval(func, lower_bound);
   y_high = feval(func, upper_bound);
   
   % Pick the half of the interval in which the root must exist
   % 'mid' is going to be next root candidate
   mid = (upper_bound - lower_bound)/2.0;
   y_mid = feval(func, mid);
   if (sign(y_mid) == sign(y_low))
      % root is in upper half
      lower_bound = mid;
   else
      % root is in lower half
      upper_bound = mid;
   end
   
   % now, check to see if we are close enough
   if (mid == 0)
      % we can't calculate fractional change here, so use criterion
      % that absolute value of the change in candidates is less than epsilon
      if (abs(mid - old_mid) < epsilon)
         done = 1;
      end
   else
      frac_change = abs((mid - old_mid)/mid);
      if (frac_change < epsilon)
         done = 1;
      end
   end
   
   % we need to save a copy of the current root candidate, so we
   % can calculate the fractional change since the previous one
   old_mid = mid;
   
   iterations = iterations + 1;
end

% when we get here, the latest root candidate is in 'mid'
root = mid;
