planck.c 5.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209
  1. /* planck.c
  2. *
  3. * Integral of Planck's black body radiation formula
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * double lambda, T, y, plancki();
  10. *
  11. * y = plancki( lambda, T );
  12. *
  13. *
  14. *
  15. * DESCRIPTION:
  16. *
  17. * Evaluates the definite integral, from wavelength 0 to lambda,
  18. * of Planck's radiation formula
  19. * -5
  20. * c1 lambda
  21. * E = ------------------
  22. * c2/(lambda T)
  23. * e - 1
  24. *
  25. * Physical constants c1 = 3.7417749e-16 and c2 = 0.01438769 are built in
  26. * to the function program. They are scaled to provide a result
  27. * in watts per square meter. Argument T represents temperature in degrees
  28. * Kelvin; lambda is wavelength in meters.
  29. *
  30. * The integral is expressed in closed form, in terms of polylogarithms
  31. * (see polylog.c).
  32. *
  33. * The total area under the curve is
  34. * (-1/8) (42 zeta(4) - 12 pi^2 zeta(2) + pi^4 ) c1 (T/c2)^4
  35. * = (pi^4 / 15) c1 (T/c2)^4
  36. * = 5.6705032e-8 T^4
  37. * where sigma = 5.6705032e-8 W m^2 K^-4 is the Stefan-Boltzmann constant.
  38. *
  39. *
  40. * ACCURACY:
  41. *
  42. * The left tail of the function experiences some relative error
  43. * amplification in computing the dominant term exp(-c2/(lambda T)).
  44. * For the right-hand tail see planckc, below.
  45. *
  46. * Relative error.
  47. * The domain refers to lambda T / c2.
  48. * arithmetic domain # trials peak rms
  49. * IEEE 0.1, 10 50000 7.1e-15 5.4e-16
  50. *
  51. */
  52. /*
  53. Cephes Math Library Release 2.8: July, 1999
  54. Copyright 1999 by Stephen L. Moshier
  55. */
  56. #include "mconf.h"
  57. #ifdef ANSIPROT
  58. extern double polylog(int, double);
  59. extern double exp(double);
  60. extern double log1p(double); /* log(1+x) */
  61. extern double expm1(double); /* exp(x) - 1 */
  62. double planckc(double, double);
  63. double plancki(double, double);
  64. #else
  65. double polylog(), exp(), log1p(), expm1();
  66. double planckc(), plancki();
  67. #endif
  68. /* NIST value (1999): 2 pi h c^2 = 3.741 7749(22) �× 10-16 W m2 */
  69. double planck_c1 = 3.7417749e-16;
  70. /* NIST value (1999): h c / k = 0.014 387 69 m K */
  71. double planck_c2 = 0.01438769;
  72. double plancki(w, T) double w, T;
  73. {
  74. double b, h, y, bw;
  75. b = T / planck_c2;
  76. bw = b * w;
  77. if (bw > 0.59375) {
  78. y = b * b;
  79. h = y * y;
  80. /* Right tail. */
  81. y = planckc(w, T);
  82. /* pi^4 / 15 */
  83. y = 6.493939402266829149096 * planck_c1 * h - y;
  84. return y;
  85. }
  86. h = exp(-planck_c2 / (w * T));
  87. y = 6. * polylog(4, h) * bw;
  88. y = (y + 6. * polylog(3, h)) * bw;
  89. y = (y + 3. * polylog(2, h)) * bw;
  90. y = (y - log1p(-h)) * bw;
  91. h = w * w;
  92. h = h * h;
  93. y = y * (planck_c1 / h);
  94. return y;
  95. }
  96. /* planckc
  97. *
  98. * Complemented Planck radiation integral
  99. *
  100. *
  101. *
  102. * SYNOPSIS:
  103. *
  104. * double lambda, T, y, planckc();
  105. *
  106. * y = planckc( lambda, T );
  107. *
  108. *
  109. *
  110. * DESCRIPTION:
  111. *
  112. * Integral from w to infinity (area under right hand tail)
  113. * of Planck's radiation formula.
  114. *
  115. * The program for large lambda uses an asymptotic series in inverse
  116. * powers of the wavelength.
  117. *
  118. * ACCURACY:
  119. *
  120. * Relative error.
  121. * The domain refers to lambda T / c2.
  122. * arithmetic domain # trials peak rms
  123. * IEEE 0.6, 10 50000 1.1e-15 2.2e-16
  124. *
  125. */
  126. double planckc(w, T) double w;
  127. double T;
  128. {
  129. double b, d, p, u, y;
  130. b = T / planck_c2;
  131. d = b * w;
  132. if (d <= 0.59375) {
  133. y = 6.493939402266829149096 * planck_c1 * b * b * b * b;
  134. return (y - plancki(w, T));
  135. }
  136. u = 1.0 / d;
  137. p = u * u;
  138. #if 0
  139. y = 236364091.*p/365866013534056632601804800000.;
  140. y = (y - 15458917./475677107995483570176000000.)*p;
  141. y = (y + 174611./123104841613737984000000.)*p;
  142. y = (y - 43867./643745871363538944000.)*p;
  143. y = ((y + 3617./1081289781411840000.)*p - 1./5928123801600.)*p;
  144. y = ((y + 691./78460462080000.)*p - 1./2075673600.)*p;
  145. y = ((((y + 1./35481600.)*p - 1.0/544320.)*p + 1.0/6720.)*p - 1./40.)*p;
  146. y = y + log(d * expm1(u));
  147. y = y - 5.*u/8. + 1./3.;
  148. #else
  149. y = -236364091. * p / 45733251691757079075225600000.;
  150. y = (y + 77683. / 352527500984795136000000.) * p;
  151. y = (y - 174611. / 18465726242060697600000.) * p;
  152. y = (y + 43867. / 107290978560589824000.) * p;
  153. y = ((y - 3617. / 202741834014720000.) * p + 1. / 1270312243200.) * p;
  154. y = ((y - 691. / 19615115520000.) * p + 1. / 622702080.) * p;
  155. y = ((((y - 1. / 13305600.) * p + 1. / 272160.) * p - 1. / 5040.) * p +
  156. 1. / 60.) *
  157. p;
  158. y = y - 0.125 * u + 1. / 3.;
  159. #endif
  160. y = y * planck_c1 * b / (w * w * w);
  161. return y;
  162. }
  163. /* planckd
  164. *
  165. * Planck's black body radiation formula
  166. *
  167. *
  168. *
  169. * SYNOPSIS:
  170. *
  171. * double lambda, T, y, planckd();
  172. *
  173. * y = planckd( lambda, T );
  174. *
  175. *
  176. *
  177. * DESCRIPTION:
  178. *
  179. * Evaluates Planck's radiation formula
  180. * -5
  181. * c1 lambda
  182. * E = ------------------
  183. * c2/(lambda T)
  184. * e - 1
  185. *
  186. */
  187. double planckd(w, T) double w, T;
  188. {
  189. return (planck_c2 / ((w * w * w * w * w) * (exp(planck_c2 / (w * T)) - 1.0)));
  190. }
  191. /* Wavelength, w, of maximum radiation at given temperature T.
  192. c2/wT = constant
  193. Wein displacement law.
  194. */
  195. double planckw(T) double T;
  196. { return (planck_c2 / (4.96511423174427630 * T)); }