''' Collection of Planck's black body radiation formulae.
.. codeauthor:: Yoshifumi Futaana
.. todo::
Move freq2wavlen, wavlen2freq, energy2wavelen, wavelen2energy to :mod:`irfpy.util.physics` module.
'''
import math
import numpy as np
from irfpy.util import unitq as u
from irfpy.util import constant as k
[docs]def freq2wavlen(freq):
''' Frequency (Hz) to wavelength (m) conversion.
:param freq: Frequency
:type freq: Scalar or np.array.
:returns: Wavelength in m
The following command confirms that 1 MHz ~ 300 m.
>>> print('%.3f' % freq2wavlen(1e6)) # 1MHz
299.792
Unitful version is also defined.
>>> freq = u.Hz * np.array([1e6, 1e7, 1e8]) # 1MHz, 10 and 100MHz
>>> wavlen = freq2wavlen_u(freq)
>>> print('%.5f %.5f %.5f' % tuple(wavlen.rescale(u.km)))
0.29979 0.02998 0.00300
'''
wavlen = k.c / freq
return wavlen
freq2wavlen_u = u.make_unitful_function(freq2wavlen, [u.Hz], u.m)
''' Unitful version of :func:`freq2wavlen`
'''
[docs]def wavlen2freq(wavlen):
''' Wavelength (m) to frequency (Hz) conversion.
:param wavlen: Wavelength in m.
:type wavlen: Sacalar or np.array.
:returns: Frequency in Hz.
>>> print('%.3e' % wavlen2freq(1e-9)) # 1 nm
2.998e+17
It confirms 1 nm ~ 300 000 000 GHz = 3e17 Hz.
Unitful version :func:`wavlen2freq_u` is also defined.
>>> wavlen = np.array([5, 10, 30])
>>> wavlen = u.km * wavlen # 5km, 10km, 30 km
>>> freq = wavlen2freq_u(wavlen)
>>> print('%.1f %.1f %.1f' % tuple(freq.rescale(u.Hz)))
59958.5 29979.2 9993.1
'''
freq = k.c / wavlen
return freq
wavlen2freq_u = u.make_unitful_function(wavlen2freq, [u.m], u.Hz)
''' Unitful version of :func:`wavlen2freq`
'''
[docs]def energy2wavlen(energy):
r''' Convert the energy (J) to wavelength (m)
>>> print('%.2e m' % (energy2wavlen(3.973e-19)))
5.00e-07 m
>>> print('3 [eV] = %.2f [nm]' % energy2wavlen_u(3 * u.eV).rescale(u.nm))
3 [eV] = 413.28 [nm]
'''
return k.h * k.c / energy
energy2wavlen_u = u.make_unitful_function(energy2wavlen, [u.J], u.m)
''' Unitful version of :func:`energy2wavlen`.
'''
[docs]def wavlen2energy(wavlen):
r''' Convert the wavelength (m) to energy (J)
.. math::
E = h\nu = \frac{hc}{\lambda}
>>> print('%.3e J' % wavlen2energy(500e-9)) # 500 nm
3.973e-19 J
>>> print('%.3f eV' % (wavlen2energy(500e-9) / k.qe))
2.480 eV
Unitful version also defined.
>>> wavlen = 500 * u.nm
>>> energy = wavlen2energy_u(wavlen)
>>> print('500 nm = %.3e J' % energy.rescale(u.J))
500 nm = 3.973e-19 J
>>> print('500 nm = %.3f eV' % energy.rescale(u.eV))
500 nm = 2.480 eV
'''
return k.h * k.c / wavlen
wavlen2energy_u = u.make_unitful_function(wavlen2energy, [u.m], u.J)
''' Unitful version of wavlen2energy
'''
[docs]class BlackBodyMksa:
''' Black body class with interfaces to MKSA unit system.
'''
lh = math.log10(k.h)
lc = math.log10(k.c)
def __init__(self, temperature):
''' A black body is instanced.
:param temperature: Temperature in K.
'''
self.tempK = temperature
self.kT = temperature * k.kB # In the unit of J.
[docs] def spectral_radiance_per_frequency(self, frequency):
r''' Return the spectral radiance in frequency domain
:param frequency: Frequency in Hz (/s).
:type frequency: Scalar or np.array
:returns: The spectral radiance, in the unit of :math:`W/m^2 sr Hz`
The spectral radiance is calcualted as
.. math::
B_\nu(T) = \frac{2h\nu^3}{c^2}\frac{1}{\exp{\frac{h\nu}{k_BT} - 1}}
>>> sun = BlackBodyMksa(6000) # Sun is 6000 K black body.
>>> print('%.3e' % sun.spectral_radiance_per_frequency(1.5e14)) # 2000 nm
2.145e-08
'''
lh, lc = self.lh, self.lc
nu = frequency
lnu = np.log10(nu)
t1 = lh + 3 * lnu - 2 * lc # h\nu^3 / c^2 in log for numerics
e1 = lh + lnu - math.log10(self.kT)
eterm = 1 / (np.exp(10 ** e1) - 1)
b = 2 * (10 ** t1) * eterm
return b
[docs] def Bnu(self, *args, **kwds):
''' Return :meth:`spectral_radiance_per_frequency`
An alias to :meth:`spectral_radiance_per_frequency`.
>>> sun = BlackBodyMksa(6000)
>>> print('%.3e' % sun.Bnu(1.5e14))
2.145e-08
'''
return self.spectral_radiance_per_frequency(*args, **kwds)
[docs] def spectral_radiance_per_wavelength(self, wavelength):
r''' Return the spectral radiance in wavelength domain
:param wavelength: Wavelength in ``m``.
:returns: The spectral radiance, in the unit of :math:`W/m^2 sr m`
.. math::
B_\lambda(T) = \frac{2hc^2}{\lambda^5}{\frac{1}{\exp{\frac{hc}{\lambda k_BT}}-1}}
>>> sun = BlackBodyMksa(6000)
>>> print('%.3e' % sun.spectral_radiance_per_wavelength(500e-9))
3.176e+13
'''
lh, lc = self.lh, self.lc
llmd = np.log10(wavelength)
t1 = lh + lc * 2 - 5 * llmd
e1 = lh + lc - llmd - math.log10(self.kT)
eterm = 1 / (np.exp(10 ** e1) - 1)
b = 2 * (10 ** t1) * eterm
return b
[docs] def Blambda(self, *args, **kwds):
''' An alias to :meth:`spectral_radiance_per_wavelength`.
'''
return self.spectral_radiance_per_wavelength(*args, **kwds)
[docs]class BlackBody_u(BlackBodyMksa):
''' Unitful version of :class:`BlackBodyMksa`.
>>> earth = BlackBody_u(300 * u.K)
>>> b = earth.Bnu(500 * u.kHz)
>>> print('{0:.3e} {1}'.format(float(b.n), b.u))
2.304e-26 W/(m**2*sr*Hz)
>>> b = earth.Blambda(500 * u.nm)
>>> print('{0:.3e} {1}'.format(float(b.n), b.u))
8.399e-27 W/(m**3*sr)
>>> sun = BlackBody_u(0.5 * u.eVT) # ~5800K
>>> b = sun.Bnu(500 * u.kHz)
>>> print('{0:.3e} {1}'.format(float(b.n), b.u))
4.457e-25 W/(m**2*sr*Hz)
>>> sun = BlackBody_u(5800 * u.K) # ~5800K
>>> b = sun.Bnu(500 * u.kHz)
>>> print('{0:.3e} {1}'.format(float(b.n), b.u))
4.455e-25 W/(m**2*sr*Hz)
'''
def __init__(self, temperature):
tempK = float(temperature.rescale(u.K))
BlackBodyMksa.__init__(self, tempK)
spectral_radiance_per_frequency = u.make_unitful_function(
BlackBodyMksa.spectral_radiance_per_frequency,
[u.Hz], u.W / (u.m **2 * u.ster * u.Hz), offset=1)
spectral_radiance_per_wavelength = u.make_unitful_function(
BlackBodyMksa.spectral_radiance_per_wavelength,
[u.m], u.W / (u.m ** 2 * u.ster * u.m), offset=1)
[docs]def SunProxy(temperature=5778.):
''' The Sun proxied by a black body.
The default temperature is 5778K, which is consdered as a average
surface temperatuer of the Sun.
:keyword temperature: Temperature in K.
>>> sun = SunProxy()
>>> print('%.3e' % sun.Blambda(500e-9)) # At 200 nm
2.638e+13
>>> lamrng = np.array([200, 500, 1000, 2000]) * 1e-9 # 200-2000 nm
>>> b = sun.Blambda(lamrng) / 1e9 # in W/m2 sr nm
>>> print('%.2e %.2e %.2e %.2e' % (b[0], b[1], b[2], b[3]))
1.46e+03 2.64e+04 1.08e+04 1.50e+03
'''
return BlackBodyMksa(temperature)