diff options
Diffstat (limited to 'libm/double/planck.c')
| -rw-r--r-- | libm/double/planck.c | 223 |
1 files changed, 223 insertions, 0 deletions
diff --git a/libm/double/planck.c b/libm/double/planck.c new file mode 100644 index 000000000..834c85dff --- /dev/null +++ b/libm/double/planck.c @@ -0,0 +1,223 @@ +/* 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 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 <math.h> +#ifdef ANSIPROT +extern double polylog (int, double); +extern double exp (double); +extern double log1p (double); /* log(1+x) */ +extern double expm1 (double); /* exp(x) - 1 */ +double planckc(double, double); +double plancki(double, double); +#else +double polylog(), exp(), 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 = 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 - 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 + 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) * (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)); +} |