1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
|
// Copyright Matt Borland, 2022
// Use, modification and distribution are subject to the
// Boost Software License, Version 1.0.
// (See accompanying file LICENSE_1_0.txt
// or copy at http://www.boost.org/LICENSE_1_0.txt)
#include <cmath>
#include <limits>
#include <boost/math/quadrature/tanh_sinh.hpp>
#include <boost/math/constants/constants.hpp>
#include "math_unit_test.hpp"
// numerically evaluate the integral for Stefan-Boltzmann Law from Planck's Law
// https://en.wikipedia.org/wiki/Stefan%E2%80%93Boltzmann_law#Derivation_from_Planck's_law
template<typename Real>
Real StefanBoltzmann()
{
constexpr auto h = Real(6.62607015e-34);
constexpr auto c = Real(299792458.0);
constexpr auto kB = Real(1.380649e-23);
constexpr auto c1 = Real(2) * boost::math::constants::pi<Real>() * h * c * c;
constexpr auto c2 = h * c / kB;
constexpr auto T = 1000;
auto integrand = [&](const Real l)
{
return c1 / (std::pow(l, Real(5)) * (std::exp(c2 / (T * l)) - Real(1)));
};
auto integrator = boost::math::quadrature::tanh_sinh<Real>(15);
return static_cast<Real>(integrator.integrate(integrand, Real(0.0), std::numeric_limits<Real>::infinity()));
}
int main() {
constexpr auto sigma = 56703.74419;
// Double Precision
CHECK_MOLLIFIED_CLOSE(sigma, StefanBoltzmann<double>(), 1e-5);
// integate using single precision. It gets the right answer, but there is integer underflow
// if this line is commented out, there is no underflow
CHECK_MOLLIFIED_CLOSE(sigma, StefanBoltzmann<float>(), 1e-5f);
return boost::math::test::report_errors();
}
|
