Below are C++ functions for computing the integrated spectral radiance (W m-2 sr-1) and integrated spectral photon radiance (photon s-1m-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 ;
} |