function integral = rect_integ(start_x, end_x, eps, func)
//
/// Integrate the given function over the given interval, using a simple
//    first-order rectangle method.
//
// The input arguments of this function are
//
//              start_x             x-value at which to start integrating
//
//              end_x               x-value at which to stop integrating
//
//              eps                 termination criterion: quit when the
//                                         fractional change in integral
//                                         drops to this value
//
//              func                name of a function which takes a single
//                                         argument
//
// How many slices (rectangles) do we use?  We start at two slices, 
//    then iterate by doubling the number of slices each time.  
//    Quit when either the fractional change in integral is less than 
//    "eps", or when we have reached some maximum number of 
//    iterations (see code).
//
// MWR 4/8/2007

  verbose = 1;

  // sanity checks
  if (end_x <= start_x)
    mprintf("rect_integ: end_x %f must be > start_x %f \n", ...
                   end_x, start_x);
    error("rect_integ: quitting now \n");
  end
  if (eps <= 0) 
    mprintf("rect_integ: eps %f must be greater than zero \n", eps);
    error("rect_integ: quitting now \n");
  end

  // quit if we take more than this many iterations
  max_iter = 100;

  // prepare for the main loop
  prev_total_area = 0;
  num_slices = 2;
  iter = 0;
  done = 0;
  while (done == 0)
    
    if (verbose > 0)
      mprintf("  iter %6d:  num_slices %6d \n", iter, num_slices);
    end

    // compute the total light using the current number of slices
    dx = (end_x - start_x) / num_slices;
    x = start_x;
    total_area = 0;
    for i = 1 : num_slices
      height = feval(x, func);
      area = height * dx; 
      total_area = total_area + area;
      if (verbose > 1) 
        mprintf("  x %9.4e height %9.4e area %9.4e = %9.4e \n", ...
                      x, height, area, total_area);    
      end
      x = x + dx;
    end
    
    // check to see if change is small enough to declare victory
    frac_change = abs((total_area - prev_total_area) / total_area);
    if (verbose > 0) 
      mprintf("   total %9.4e  frac_change %9.4e  eps %9.4e   ", ...
                  total_area, frac_change, eps);
    end
    if (frac_change < eps) 
      if (verbose > 0) 
        mprintf("  done! \n");
      end
      done = 1;
    else
      if (verbose > 0) 
        mprintf("  keep going \n");
      end
    end

    // quit if it has taken too long
    if (iter > max_iter)
      mprintf("rect_integ: reached limit %5d iter \n", max_iter);
      error("rect_integ: quitting now \n");
    end

    // prepare for the next time through the loop
    iter = iter + 1;
    num_slices = num_slices * 2;
    prev_total_area = total_area;

  end

  integral = total_area;

endfunction

