''' Solar spectrum models at the Earth orbit
.. codeauthor:: Yoshifumi Futaana
.. autosummary::
BlackBody
BlackBodyModFuv
AM0
Solar spectrum models at the Earth orbit (1 AU).
See also :mod:`irfpy.util.planck`.
'''
import os
import logging
_logger = logging.getLogger(__name__)
from pkg_resources import resource_filename
import numpy as np
import scipy.interpolate
from irfpy.util import planck
[docs]class BlackBody:
r''' Black body assumption of the Sun.
Sun is proxied by a black body with tempereature 5778 K.
Using the planck's formula (see :class:`irfpy.util.planck.BlackBodyMksa`)
the spectrum at the Sun, i.e. power by area by wavelength [W / m2 / sr/ m],
is immediately known as :math:`B(\lambda, T)`.
This is implemented in :mod:`irfpy.util.planck`.
The total spectral power is the integral over the Sun surface area and
the emission angle (in solid angle).
As isotropic emission is assumed, the solid angle integral can be done
.. math::
\int{\cos\theta d\Omega} = \int_0^{2\pi}d\phi\int_0^{\pi/2}\cos\theta\sin\theta d\theta = \pi
The :math:`\cos\theta` appears from the inner product of the flux and
the surface area. See also how to derive Stefan-Boltzmann law.
Same calculation is done.
Anyway, the total spectral power, in the unit of [W / m], is
.. math::
F(\lambda, T) = \pi B(\lambda, T) \times 4\pi Rs^2
Here ``Rs`` is the radius of the Sun.
Consider a unit area at the Earth's orbit. The spectral power through
the area is, as the emission is isotropic, calculated by multiplying the
fraction of the area to the sphare surface at 1AU, namely,
.. math::
u(\lambda, T) = F(\lambda, T) \times \frac{1}{4\pi a^2}
where ``a`` is 1AU. The unit of ``u`` is [W / m2 / m].
In summary, the spectrum power at the Earth orbit, ``u`` [W / m2 / m] is
.. math::
u(\lambda, T) = \pi \frac{Rs^2}{a^2} B(\lambda, T)
Substituting Rs=695500 km and a=149597871 km,
.. math::
u(\lambda, T) = 6.8687\times10^5 B(\lambda, T)
'''
def __init__(self, temperature=5778.):
self.sun = planck.SunProxy(temperature=temperature)
self.scale = 6.8687e-5 # Scaling from the spectrum at the sun (W / m2 / sr / m)
# to the spectrum at the Earth (W / m2 / m).
# sr is considered into the scale.
self.name = 'BlackBody(%gK)' % temperature
''' Name of the model '''
[docs] def spectral_irradiance(self, wavelength):
''' Return the spectral irradiance (W / m2 / m).
Return the spectral irradiance, in the unit of W / m2 / m.
:param wavelength: Wavelength in m.
:returns: Spectral irradiance in W / m2 / m
>>> sun = BlackBody(5800) # 5800 K blackbody approx.
>>> irr500nm = sun.spectral_irradiance(500e-9)
>>> print('%.3e' % irr500nm)
1.846e+09
'''
spectral_radiance_at_sun = self.sun.Blambda(wavelength) # in W/m2 sr m
spectral_irradiance_at_earth = spectral_radiance_at_sun * self.scale
return spectral_irradiance_at_earth
[docs] def photon_flux(self, wavelength):
r''' Return the photon's differential flux (# / m2 / s / m).
The spectral irradiance, u [W / m2 / m = J / m2 / s / m],
can be converted to the photon flux.
A photon with wavelength, :math:`\lambda`, will have the energy
:math:`E = hc/\lambda`.
Thus, the number of photons through a unit area per unit time is
:math:`F = u / E`.
>>> sun = BlackBody(5800)
>>> fph = sun.photon_flux(500e-9) # in /m2 /s /m
>>> fph = fph / 1e9 # in m2 /s /nm
>>> print('%.3e' % fph)
4.648e+18
'''
u = self.spectral_irradiance(wavelength)
ene = planck.wavlen2energy(wavelength)
f = u / ene
return f
[docs]class BlackBodyModFuv(BlackBody):
''' Black body assumption with constant FUV (<~100 nm) modification
'''
def __init__(self, *args, **kwds):
BlackBody.__init__(self, *args, **kwds)
self.name = self.name.replace('BlackBody', 'BlackBody+constFUV')
''' Name of the model '''
[docs] def spectral_irradiance(self, wavelength):
'''
>>> sun = BlackBodyModFuv(5800)
>>> print('%.3e' % (sun.spectral_irradiance(500e-9)))
1.846e+09
'''
irr = BlackBody.spectral_irradiance(self, wavelength)
irrarr = np.array(irr)
irrarr = np.where(np.logical_and(irrarr < 1e5, wavelength < 500e-9), 1e5, irrarr)
# if irr < 1e5 and wavelength < 500e-9:
# return 1e5
# else:
# return irr
return irrarr
[docs]class AM0:
''' Air mass 0 model (at extraterrestorial) under ASTM E-490.
>>> am0 = AM0()
>>> irr = am0.spectral_irradiance(500e-9)
>>> print('%.3e' % irr)
1.914e+09
>>> fph = am0.photon_flux(np.array([500e-9, 1000e-9]))
>>> print('%.3e %.3e' % (fph[0] / 1e9, fph[1] / 1e9)) # in m2 /nm /s
4.816e+18 3.765e+18
'''
logger = logging.getLogger(__name__ + '.AM0')
def __init__(self):
fn = resource_filename(__name__,
os.path.join('data', 'astm_e490_am0.txt.gz'))
self.logger.debug('Filename %s, exists: %s' % (fn,
os.path.exists(fn)))
dat = np.loadtxt(fn)
self.logger.debug('Data size = %s' % str(dat.shape))
wavlen = dat[:, 0] # Data file in nm
irrad = dat[:, 1] # Data file in W/m2/nm
self.irradfunc = scipy.interpolate.interp1d(wavlen, irrad)
''' The irradiance function from wavelength (nm) to irradiance (W/m2/nm)
'''
self.name = 'ASTM/E490/AM0'
[docs] def spectral_irradiance(self, wavelength):
''' Return the irradiance with interface of MKSA
'''
wavlen_nm = np.array(wavelength) * 1e9 # m -> nm.
irrad_nm = self.irradfunc(wavlen_nm)
self.logger.debug('%s' % str(irrad_nm))
irrad_mksa = irrad_nm * 1e9 # W/m2/nm -> W/m2/m
return irrad_mksa
[docs] def photon_flux(self, wavelength):
u = self.spectral_irradiance(wavelength)
ene = planck.wavlen2energy(wavelength)
f = u / ene
return f
import unittest
import doctest
[docs]def doctests():
return unittest.TestSuite((
doctest.DocTestSuite(),
))
if __name__ == '__main__':
unittest.main(defaultTest='doctests')