Appendix A: Algorithms for Computing In-band Radiance

Below are C++ functions for computing the integrated spectral radiance (W·m^-2·sr^-1) and integrated spectral photon radiance (photon s^-1·m^-2·sr^-1). The functions compute the integral from the specified wavenumber to infinity for a blackbody at the input temperature. Finite spectral regions can be computed by using this function twice-once with each end point of the spectral region. The difference of the two gives the radiance for the spectral region.

#include  <math.h> // for “exp” function

double planck_integral (double sigma, double temperature) {

//  integral of spectral radiance from sigma (cm-1) to infinity.

//  result is W/m2/sr.

//  follows Widger and Woodall, Bulletin of the American Meteorological

//  Society, Vol. 57, No. 10, pp. 1217

//  constants

      double Planck =  6.6260693e-34 ;    

      double  Boltzmann = 1.380658e-23 ;

      double  Speed_of_light = 299792458.0 ;

      double  Speed_of_light_sq = Speed_of_light * Speed_of_light ;

//  compute powers of x, the dimensionless spectral coordinate

      double c1 =  (Planck*Speed_of_light/Boltzmann) ;

      double x =  c1 * 100 * sigma / temperature ;

      double x2 = x *  x  ;

      double x3 = x *  x2 ;

//  decide how many terms of sum are needed

      double iterations = 2.0 + 20.0/x ;

      iterations = (iterations<512) ? iterations : 512 ;

      int iter = int(iterations) ;

//  add up terms of sum

      double sum = 0  ;

      for (int n=1;  n<iter; n++) {

            double  dn = 1.0/n ;

            sum  += exp(-n*x)*(x3 + (3.0 * x2 + 6.0*(x+dn)*dn)*dn)*dn;

            }

//  return result, in units of W/m2/sr

      double c2 =  (2.0*Planck*Speed_of_light_sq) ;

      return  c2*pow(temperature/c1,4)*sum ;

}

#include  <math.h>  // for “exp”  function

double planck_photon_integral (double sigma, double temperature) {

//  integral of spectral photon radiance from sigma (cm-1) to infinity.

//  result is photons/s/m2/sr.

//  follows Widger and Woodall, Bulletin of the American Meteorological

//  Society, Vol. 57, No. 10, pp. 1217

//  constants

      double Planck =  6.6260693e-34 ;    

      double Boltzmann = 1.380658e-23 ;

      double Speed_of_light = 299792458.0 ;

//  compute powers of x, the dimensionless spectral coordinate

      double c1 =  Planck*Speed_of_light/Boltzmann ;

      double x =  c1*100*sigma/temperature ;

      double x2 = x *  x  ;

//  decide how many terms of sum are needed

      double iterations = 2.0 + 20.0/x ;

      iterations = (iterations<512) ? iterations : 512 ;

      int iter = int(iterations) ;

//  add up terms of sum

      double sum = 0  ;

      for (int n=1;  n<iter; n++) {

            double  dn = 1.0/n ;

            sum  += exp(-n*x) * (x2 + 2.0*(x + dn)*dn)*dn ;

            }

//  return result, in units of photons/s/m2/sr

      double kTohc =  Boltzmann*temperature/(Planck*Speed_of_light) ;

      double c2 =  2.0* pow(kTohc,3)*Speed_of_light ;

      return c2 *sum  ;

}