function peak_lambda = planck_peak(start_lambda, end_lambda, temp, epsilon)
//
// Given an interval in wavelength and a temperature, find the peak wavelength
// of the Planck function within the interval.  Iterate until the fractional
// change in peak wavelength is less than epsilon.
//
// 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 
//                                              peak wavelength at which we
//                                              stop looking
//
//             peak_lambda     (output)     peak wavelength of Planck function
//                                              within the interval (m)
//
// We'll use a bisection approach to find the peak.  At each iteration,
// we use the value of the slope at the current wavelength in order to
// decide if the peak lies in the lower or upper half of the current 
// interval.  In order to find the slope, we evaluate the Planck function
// at the current wavelength, and at another neighboring wavelength.
//
// MWR 3/31/2002
//

// check input values for validity

// this is the name of the Planck function we'll call

// in order to find the slope near the current wavelength, we'll check
//   the function's value at the current wavelength, AND at a wavelength
//   which is this fraction of the current interval away from it.

// prepare to enter main loop in which we search for the peak

// enter the main loop

  // use the middle of the current interval as the current guess

  // pick a neighboring wavelength

  // evaluate the Planck function at the current wavelength, and at
  // the neighboring wavelength
  
  // decide which half of the current interval contains the peak
  //   and move the appropriate boundary
  //  (and check to see if we're at the peak right now)
  
  
  // check to see if the fractional change since last time is small
  //   enough that we can quit
  
// end of main loop
  
// set the output argument "planck_peak" to the current wavelength