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
//

verbose = 0;

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

  
// 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.
neighbor_fraction = 0.01;

// 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
bot = start_lambda;
top = end_lambda;
prev_mid = bot;
num_iter = 0;

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

  if (verbose > 0)
    mprintf(" iter %5d:  bot %9.5e    top %9.5e \n", num_iter, bot, top);
  end

  // use the middle of the current interval as the current guess
  mid = (bot + top)/2.0;

  // pick a neighboring wavelength
  neighbor = mid + (top - bot)*neighbor_fraction;

  // evaluate the Planck function at the current wavelength, and at
  // the neighboring wavelength
  mid_planck = feval(temp, mid, planck);
  neighbor_planck = feval(temp, neighbor, planck);
  if (verbose > 0)
    mprintf("    mid %9.4e = %9.4e   neighbor %9.4e = %9.4e ", ...
             mid, mid_planck, neighbor, neighbor_planck);
  end
  
  // 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)
  if (mid_planck > neighbor_planck)
    if (verbose > 0)
      mprintf(" move top \n");
    end
    top = mid;
  else
    if (verbose > 0)
      mprintf(" move bot \n");
    end
    bot = mid;
  end
  
  // check to see if the fractional change since last time is small
  //   enough that we can quit
  frac_change = abs((mid - prev_mid)/mid);
  if (verbose > 0)
    mprintf("  new %9.4e  old %9.4e  frac %9.4e ", ...
               mid, prev_mid, 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_peak: reached limit of %d iterations \n", max_iter);
    end
    done = 1;
  end

  prev_mid = mid;
  num_iter = num_iter + 1;
  
// end of main loop
end
  
// set the output argument "planck_peak" to the current wavelength
peak_lambda = mid;


endfunction
