/* planck.c
*
* Integral of Planck's black body radiation formula
*
*
*
* SYNOPSIS:
*
* double lambda, T, y, plancki();
*
* y = plancki( lambda, T );
*
*
*
* DESCRIPTION:
*
* Evaluates the definite integral, from wavelength 0 to lambda,
* of Planck's radiation formula
* -5
* c1 lambda
* E = ------------------
* c2/(lambda T)
* e - 1
*
* Physical constants c1 = 3.7417749e-16 and c2 = 0.01438769 are built in
* to the function program. They are scaled to provide a result
* in watts per square meter. Argument T represents temperature in degrees
* Kelvin; lambda is wavelength in meters.
*
* The integral is expressed in closed form, in terms of polylogarithms
* (see polylog.c).
*
* The total area under the curve is
* (-1/8) (42 zeta(4) - 12 pi^2 zeta(2) + pi^4 ) c1 (T/c2)^4
* = (pi^4 / 15) c1 (T/c2)^4
* = 5.6705032e-8 T^4
* where sigma = 5.6705032e-8 W m^2 K^-4 is the Stefan-Boltzmann constant.
*
*
* ACCURACY:
*
* The left tail of the function experiences some relative error
* amplification in computing the dominant term md_exp(-c2/(lambda T)).
* For the right-hand tail see planckc, below.
*
* Relative error.
* The domain refers to lambda T / c2.
* arithmetic domain # trials peak rms
* IEEE 0.1, 10 50000 7.1e-15 5.4e-16
*
*/
/*
Cephes Math Library Release 2.8: July, 1999
Copyright 1999 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef ANSIPROT
extern double polylog (int, double);
extern double md_exp (double);
extern double md_log1p (double); /* md_log(1+x) */
extern double expm1 (double); /* md_exp(x) - 1 */
double planckc(double, double);
double plancki(double, double);
#else
double polylog(), md_exp(), md_log1p(), expm1();
double planckc(), plancki();
#endif
/* NIST value (1999): 2 pi h c^2 = 3.741 7749(22) × 10-16 W m2 */
double planck_c1 = 3.7417749e-16;
/* NIST value (1999): h c / k = 0.014 387 69 m K */
double planck_c2 = 0.01438769;
double
plancki(w, T)
double w, T;
{
double b, h, y, bw;
b = T / planck_c2;
bw = b * w;
if (bw > 0.59375)
{
y = b * b;
h = y * y;
/* Right tail. */
y = planckc (w, T);
/* pi^4 / 15 */
y = 6.493939402266829149096 * planck_c1 * h - y;
return y;
}
h = md_exp(-planck_c2/(w*T));
y = 6. * polylog (4, h) * bw;
y = (y + 6. * polylog (3, h)) * bw;
y = (y + 3. * polylog (2, h)) * bw;
y = (y - md_log1p (-h)) * bw;
h = w * w;
h = h * h;
y = y * (planck_c1 / h);
return y;
}
/* planckc
*
* Complemented Planck radiation integral
*
*
*
* SYNOPSIS:
*
* double lambda, T, y, planckc();
*
* y = planckc( lambda, T );
*
*
*
* DESCRIPTION:
*
* Integral from w to infinity (area under right hand tail)
* of Planck's radiation formula.
*
* The program for large lambda uses an asymptotic series in inverse
* powers of the wavelength.
*
* ACCURACY:
*
* Relative error.
* The domain refers to lambda T / c2.
* arithmetic domain # trials peak rms
* IEEE 0.6, 10 50000 1.1e-15 2.2e-16
*
*/
double
planckc (w, T)
double w;
double T;
{
double b, d, p, u, y;
b = T / planck_c2;
d = b*w;
if (d <= 0.59375)
{
y = 6.493939402266829149096 * planck_c1 * b*b*b*b;
return (y - plancki(w,T));
}
u = 1.0/d;
p = u * u;
#if 0
y = 236364091.*p/365866013534056632601804800000.;
y = (y - 15458917./475677107995483570176000000.)*p;
y = (y + 174611./123104841613737984000000.)*p;
y = (y - 43867./643745871363538944000.)*p;
y = ((y + 3617./1081289781411840000.)*p - 1./5928123801600.)*p;
y = ((y + 691./78460462080000.)*p - 1./2075673600.)*p;
y = ((((y + 1./35481600.)*p - 1.0/544320.)*p + 1.0/6720.)*p - 1./40.)*p;
y = y + md_log(d * expm1(u));
y = y - 5.*u/8. + 1./3.;
#else
y = -236364091.*p/45733251691757079075225600000.;
y = (y + 77683./352527500984795136000000.)*p;
y = (y - 174611./18465726242060697600000.)*p;
y = (y + 43867./107290978560589824000.)*p;
y = ((y - 3617./202741834014720000.)*p + 1./1270312243200.)*p;
y = ((y - 691./19615115520000.)*p + 1./622702080.)*p;
y = ((((y - 1./13305600.)*p + 1./272160.)*p - 1./5040.)*p + 1./60.)*p;
y = y - 0.125*u + 1./3.;
#endif
y = y * planck_c1 * b / (w*w*w);
return y;
}
/* planckd
*
* Planck's black body radiation formula
*
*
*
* SYNOPSIS:
*
* double lambda, T, y, planckd();
*
* y = planckd( lambda, T );
*
*
*
* DESCRIPTION:
*
* Evaluates Planck's radiation formula
* -5
* c1 lambda
* E = ------------------
* c2/(lambda T)
* e - 1
*
*/
double
planckd(w, T)
double w, T;
{
return (planck_c2 / ((w*w*w*w*w) * (md_exp(planck_c2/(w*T)) - 1.0)));
}
/* Wavelength, w, of maximum radiation at given temperature T.
c2/wT = constant
Wein displacement law.
*/
double
planckw(T)
double T;
{
return (planck_c2 / (4.96511423174427630 * T));
}