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 (epsilon <= 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(lower_bound, func);
y_high = feval(upper_bound, func);
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(lower_bound, func);
   y_high = feval(upper_bound, func);
   
   // 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(mid, func);
   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;

endfunction