# Source code for astropy.modeling.blackbody

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
Model and functions related to blackbody radiation.

.. _blackbody-planck-law:

-------------------

Blackbody flux is calculated with Planck law
(:ref:Rybicki & Lightman 1979 <ref-rybicki1979>):

.. math::

B_{\\lambda}(T) = \\frac{2 h c^{2} / \\lambda^{5}}{exp(h c / \\lambda k T) - 1}

B_{\\nu}(T) = \\frac{2 h \\nu^{3} / c^{2}}{exp(h \\nu / k T) - 1}

where the unit of :math:B_{\\lambda}(T) is
:math:erg \\; s^{-1} cm^{-2} \\mathring{A}^{-1} sr^{-1}, and
:math:B_{\\nu}(T) is :math:erg \\; s^{-1} cm^{-2} Hz^{-1} sr^{-1}.
:func:~astropy.modeling.blackbody.blackbody_lambda and
:func:~astropy.modeling.blackbody.blackbody_nu calculate the
blackbody flux for :math:B_{\\lambda}(T) and :math:B_{\\nu}(T),
respectively.

For blackbody representation as a model, see :class:BlackBody1D.

.. _blackbody-examples:

Examples
^^^^^^^^

>>> import numpy as np
>>> from astropy import units as u
>>> from astropy.modeling.blackbody import blackbody_lambda, blackbody_nu

Calculate blackbody flux for 5000 K at 100 and 10000 Angstrom while suppressing
any Numpy warnings:

>>> wavelengths = [100, 10000] * u.AA
>>> temperature = 5000 * u.K
>>> with np.errstate(all='ignore'):
...     flux_lam = blackbody_lambda(wavelengths, temperature)
...     flux_nu = blackbody_nu(wavelengths, temperature)
>>> flux_lam  # doctest: +FLOAT_CMP
<Quantity [  1.27452545e-108,  7.10190526e+005] erg / (Angstrom cm2 s sr)>
>>> flux_nu  # doctest: +FLOAT_CMP
<Quantity [  4.25135927e-123,  2.36894060e-005] erg / (cm2 Hz s sr)>

Alternatively, the same results for flux_nu can be computed using
:class:BlackBody1D with blackbody representation as a model. The difference between
this and the former approach is in one additional step outlined as follows:

>>> from astropy import constants as const
>>> from astropy.modeling import models
>>> temperature = 5000 * u.K
>>> bolometric_flux = const.sigma_sb * temperature ** 4 / np.pi
>>> bolometric_flux.to(u.erg / (u.cm * u.cm * u.s))  # doctest: +FLOAT_CMP
<Quantity 1.12808367e+10 erg / (cm2 s)>
>>> wavelengths = [100, 10000] * u.AA
>>> bb_astro = models.BlackBody1D(temperature, bolometric_flux=bolometric_flux)
>>> bb_astro(wavelengths).to(u.erg / (u.cm * u.cm * u.Hz * u.s)) / u.sr  # doctest: +FLOAT_CMP
<Quantity [4.25102471e-123, 2.36893879e-005] erg / (cm2 Hz s sr)>

where bb_astro(wavelengths) computes the equivalent result as flux_nu above.

Plot a blackbody spectrum for 5000 K:

.. plot::

import matplotlib.pyplot as plt
import numpy as np
from astropy import constants as const
from astropy import units as u
from astropy.modeling.blackbody import blackbody_lambda

temperature = 5000 * u.K
wavemax = (const.b_wien / temperature).to(u.AA)  # Wien's displacement law
waveset = np.logspace(
0, np.log10(wavemax.value + 10 * wavemax.value), num=1000) * u.AA
with np.errstate(all='ignore'):
flux = blackbody_lambda(waveset, temperature)

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(waveset.value, flux.value)
ax.axvline(wavemax.value, ls='--')
ax.get_yaxis().get_major_formatter().set_powerlimits((0, 1))
ax.set_xlabel(r'$\\lambda$ ({0})'.format(waveset.unit))
ax.set_ylabel(r'$B_{\\lambda}(T)$')
ax.set_title('Blackbody, T = {0}'.format(temperature))

Note that an array of temperatures can also be given instead of a single
temperature. In this case, the Numpy broadcasting rules apply: for instance, if
the frequency and temperature have the same shape, the output will have this
shape too, while if the frequency is a 2-d array with shape (n, m) and the
temperature is an array with shape (m,), the output will have a shape
(n, m).

^^^^^^^^

.. _ref-rybicki1979:

Rybicki, G. B., & Lightman, A. P. 1979, Radiative Processes in Astrophysics (New York, NY: Wiley)

"""

import warnings
from collections import OrderedDict

import numpy as np

from .core import Fittable1DModel
from .parameters import Parameter
from astropy import constants as const
from astropy import units as u
from astropy.utils.exceptions import AstropyUserWarning

__all__ = ['BlackBody1D', 'blackbody_nu', 'blackbody_lambda']

# Units
FNU = u.erg / (u.cm**2 * u.s * u.Hz)
FLAM = u.erg / (u.cm**2 * u.s * u.AA)

# Some platform implementations of expm1() are buggy and Numpy uses
# them anyways--the bug is that on certain large inputs it returns
# NaN instead of INF like it should (it should only return NaN on a
# NaN input
# See https://github.com/astropy/astropy/issues/4171
with warnings.catch_warnings():
warnings.simplefilter('ignore', RuntimeWarning)
_has_buggy_expm1 = np.isnan(np.expm1(1000)) or np.isnan(np.expm1(1e10))

[docs]class BlackBody1D(Fittable1DModel):
"""
One dimensional blackbody model.

Parameters
----------
temperature : :class:~astropy.units.Quantity
Blackbody temperature.
bolometric_flux : :class:~astropy.units.Quantity
The bolometric flux of the blackbody (i.e., the integral over the
spectral axis).

Notes
-----

Model formula:

.. math:: f(x) = \\pi B_{\\nu} f_{\\text{bolometric}} / (\\sigma  T^{4})

Examples
--------
>>> from astropy.modeling import models
>>> from astropy import units as u
>>> bb = models.BlackBody1D()
>>> bb(6000 * u.AA)  # doctest: +FLOAT_CMP
<Quantity 1.3585381201978953e-15 erg / (cm2 Hz s)>

.. plot::
:include-source:

import numpy as np
import matplotlib.pyplot as plt

from astropy.modeling.models import BlackBody1D
from astropy.modeling.blackbody import FLAM
from astropy import units as u
from astropy.visualization import quantity_support

bb = BlackBody1D(temperature=5778*u.K)
wav = np.arange(1000, 110000) * u.AA
flux = bb(wav).to(FLAM, u.spectral_density(wav))

with quantity_support():
plt.figure()
plt.semilogx(wav, flux)
plt.axvline(bb.lambda_max.to(u.AA).value, ls='--')
plt.show()

"""

# We parametrize this model with a temperature and a bolometric flux. The
# bolometric flux is the integral of the model over the spectral axis. This
# is more useful than simply having an amplitude parameter.
temperature = Parameter(default=5000, min=0, unit=u.K)
bolometric_flux = Parameter(default=1, min=0, unit=u.erg / u.cm ** 2 / u.s)

# We allow values without units to be passed when evaluating the model, and
# in this case the input x values are assumed to be frequencies in Hz.
_input_units_allow_dimensionless = True

# We enable the spectral equivalency by default for the spectral axis
input_units_equivalencies = {'x': u.spectral()}

[docs]    def evaluate(self, x, temperature, bolometric_flux):
"""Evaluate the model.

Parameters
----------
x : float, ~numpy.ndarray, or ~astropy.units.Quantity
Frequency at which to compute the blackbody. If no units are given,
this defaults to Hz.

temperature : float, ~numpy.ndarray, or ~astropy.units.Quantity
Temperature of the blackbody. If no units are given, this defaults
to Kelvin.

bolometric_flux : float, ~numpy.ndarray, or ~astropy.units.Quantity
Desired integral for the blackbody.

Returns
-------
y : number or ndarray
Blackbody spectrum. The units are determined from the units of
bolometric_flux.
"""

# We need to make sure that we attach units to the temperature if it
# doesn't have any units. We do this because even though blackbody_nu
# can take temperature values without units, the / temperature ** 4
# factor needs units to be defined.
if isinstance(temperature, u.Quantity):
temperature = temperature.to(u.K, equivalencies=u.temperature())
else:
temperature = u.Quantity(temperature, u.K)

# We normalize the returned blackbody so that the integral would be
# unity, and we then multiply by the bolometric flux. A normalized
# blackbody has f_nu = pi * B_nu / (sigma * T^4), which is what we
# calculate here. We convert to 1/Hz to make sure the units are
# simplified as much as possible, then we multiply by the bolometric
# flux to get the normalization right.
fnu = ((np.pi * u.sr * blackbody_nu(x, temperature) /
const.sigma_sb / temperature ** 4).to(1 / u.Hz) *
bolometric_flux)

# If the bolometric_flux parameter has no unit, we should drop the /Hz
# and return a unitless value. This occurs for instance during fitting,
# since we drop the units temporarily.
if hasattr(bolometric_flux, 'unit'):
return fnu
else:
return fnu.value

@property
def input_units(self):
# The input units are those of the 'x' value, which should always be
# Hz. Because we do this, and because input_units_allow_dimensionless
# is set to True, dimensionless values are assumed to be in Hz.
return {'x': u.Hz}

def _parameter_units_for_data_units(self, inputs_unit, outputs_unit):
return OrderedDict([('temperature', u.K),
('bolometric_flux', outputs_unit['y'] * u.Hz)])

@property
def lambda_max(self):
"""Peak wavelength when the curve is expressed as power density."""
return const.b_wien / self.temperature

[docs]def blackbody_nu(in_x, temperature):
"""Calculate blackbody flux per steradian, :math:B_{\\nu}(T).

.. note::

Use numpy.errstate to suppress Numpy warnings, if desired.

.. warning::

Output values might contain nan and inf.

Parameters
----------
in_x : number, array-like, or ~astropy.units.Quantity
Frequency, wavelength, or wave number.
If not a Quantity, it is assumed to be in Hz.

temperature : number, array-like, or ~astropy.units.Quantity
Blackbody temperature.
If not a Quantity, it is assumed to be in Kelvin.

Returns
-------
flux : ~astropy.units.Quantity
Blackbody monochromatic flux in
:math:erg \\; cm^{-2} s^{-1} Hz^{-1} sr^{-1}.

Raises
------
ValueError
Invalid temperature.

ZeroDivisionError
Wavelength is zero (when converting to frequency).

"""
# Convert to units for calculations, also force double precision
freq = u.Quantity(in_x, u.Hz, dtype=np.float64)
temp = u.Quantity(temperature, u.K, dtype=np.float64)

# Check if input values are physically possible
if np.any(temp < 0):
raise ValueError('Temperature should be positive: {0}'.format(temp))
if not np.all(np.isfinite(freq)) or np.any(freq <= 0):
warnings.warn('Input contains invalid wavelength/frequency value(s)',
AstropyUserWarning)

log_boltz = const.h * freq / (const.k_B * temp)
boltzm1 = np.expm1(log_boltz)

if _has_buggy_expm1:
# Replace incorrect nan results with infs--any result of 'nan' is
# incorrect unless the input (in log_boltz) happened to be nan to begin
# with.  (As noted in #4393 ideally this would be replaced by a version
# of expm1 that doesn't have this bug, rather than fixing incorrect
# results after the fact...)
boltzm1_nans = np.isnan(boltzm1)
if np.any(boltzm1_nans):
if boltzm1.isscalar and not np.isnan(log_boltz):
boltzm1 = np.inf
else:
boltzm1[np.where(~np.isnan(log_boltz) & boltzm1_nans)] = np.inf

# Calculate blackbody flux
bb_nu = (2.0 * const.h * freq ** 3 / (const.c ** 2 * boltzm1))
flux = bb_nu.to(FNU, u.spectral_density(freq))

return flux / u.sr  # Add per steradian to output flux unit

[docs]def blackbody_lambda(in_x, temperature):
"""Like :func:blackbody_nu but for :math:B_{\\lambda}(T).

Parameters
----------
in_x : number, array-like, or ~astropy.units.Quantity
Frequency, wavelength, or wave number.
If not a Quantity, it is assumed to be in Angstrom.

temperature : number, array-like, or ~astropy.units.Quantity
Blackbody temperature.
If not a Quantity, it is assumed to be in Kelvin.

Returns
-------
flux : ~astropy.units.Quantity
Blackbody monochromatic flux in
:math:erg \\; cm^{-2} s^{-1} \\mathring{A}^{-1} sr^{-1}.

"""
if getattr(in_x, 'unit', None) is None:
in_x = u.Quantity(in_x, u.AA)

bb_nu = blackbody_nu(in_x, temperature) * u.sr  # Remove sr for conversion
flux = bb_nu.to(FLAM, u.spectral_density(in_x))

return flux / u.sr  # Add per steradian to output flux unit