function planck_integ = planck_integ(start_lambda, end_lambda, temp, epsilon)
//
// Given an interval in wavelength and a temperature, integrate the Planck
// function for the given temperature over the interval.  We use a simple
// first-order method: break the interval into equal slices, compute the
// area of each slice.
//
// We iterate by doubling the number of slices each time through the
// main loop.  When the new area differs from the previous area by
// a fractional amount less than "epsilon", we declare victory and return
// the result.
//
// Arguments:    
//             start_lambda     (input)     starting point of interval (m)
//
//             end_lambda       (input)     ending point of interval (m)
//
//             temp             (input)     temperature of Planck function (K)
//
//             epsilon          (input)     fractional change in guessed 
//                                              integrated area at which we
//                                              stop looking
//
//             planck_integ     (output)     integrated energy per unit area
//                                              per second per steradian
//                                              radiated over the interval 
//
// MWR 4/6/2007
//

verbose = 0;

// check input values for validity
if (start_lambda >= end_lambda) 
  mprintf("planck_integ: start_lambda must be < end_lambda \n");
  error("planck_integ: quitting \n");
end
if ((start_lambda < 0) | (end_lambda < 0)) 
  mprintf("planck_integ: both wavelength limits must be positive \n");
  error("planck_integ: quitting \n");
end
if (temp <= 0)
  mprintf("planck_integ: temp must be positive \n");
  error("planck_integ: quitting \n");
end
if (epsilon <= 0)
  mprintf("planck_integ: epsilon (termination criterion) must be positive \n");
  error("planck_integ: quitting \n");
end

// if we iterate more than this many times, something is wrong --
//   stop and return the current best guess
max_iter = 100;

// prepare to enter main loop in which we search for the peak
prev_area = 0;
num_slices = 1;
num_iter = 0;

// enter the main loop
done = 0;
while (done == 0)

  if (verbose > 0)
    mprintf(" iter %5d:  num_slices %5d  \n", num_iter, num_slices);
  end

  // how wide is each slice?
  slice_width = (end_lambda - start_lambda)/num_slices;

  // break the interval into rectangles; compute the height and area
  //   of each rectangle, accumulate the area in the "total_area" variable.
  //   We will use the mid-point of each rectangle to find its height.
  total_area = 0;
  lambda = start_lambda + 0.5*slice_width;
  for i = 1 : num_slices
    height = feval(temp, lambda, planck);
    // one could be more efficient here by adding up the heights,
    //    then multiplying by width just once at the end ...
    area = height*slice_width;
    if (verbose > 1)
      mprintf("  slice %6d  lam %9.4e  height %9.4e  area %9.4e = %9.4e \n", ...
                 i, lambda, height, area, total_area);
    end
    total_area = total_area + area;
    lambda = lambda + slice_width;
  end
  if (verbose > 0)
    mprintf("    total area %9.4e \n", total_area);
  end
  
  // check to see if the fractional change since last time is small
  //   enough that we can quit
  if (total_area == 0) 
    frac_change = abs(total_area - prev_area);
  else
    frac_change = abs((total_area - prev_area)/total_area);
  end
  if (verbose > 0)
    mprintf("  new %9.4e   prev %9.4e  frac %9.4e ", ...
               total_area, prev_area, frac_change);
  end
  if (frac_change <= epsilon)
    if (verbose > 0)
      mprintf("  done! \n");
    end
    done = 1;
  else
    if (verbose > 0)
      mprintf("  keep going \n");
    end
  end

  // check to see if we should quit because it has taken too long
  if (num_iter > max_iter) 
    if (verbose > 0)
      mprintf("  planck_integ: reached limit of %d iterations \n", max_iter);
    end
    done = 1;
  end

  num_slices = num_slices*2;
  prev_area = total_area;
  num_iter = num_iter + 1;
  
  
// end of main loop
end
  
// set the output argument "planck_integ" to the current integrated area
planck_integ = total_area;


endfunction
