Radiometry

Radiometric source functions, the lowest layer every emitting or measuring stage reaches (#246): the Planck radiance B(ν̃, T) per unit wavenumber in W m⁻² sr⁻¹ (cm⁻¹)⁻¹, its log, their closed-form derivatives, the radiation constants, and the per-µm and per-GHz views. Wavenumber is the one internal axis; wavelength and frequency are views. Brightness temperature inverts Planck: in closed form per wavenumber, and per band of an Instrument by Newton with implicit differentiation.

The Planck function is checked against independent oracles (gate TH12): the constants against the SI 2019 exact h, c and k, the Stefan–Boltzmann law from the spectral integral, the Wien peaks, and pyradi's ryplanck (radiance, derivatives and Widger–Woodall band integrals, rtol 1e-12). pyradi is never imported: scripts/make_planck_oracle.py evaluates ryplanck.py at commit c7cfe9df96fdf0e6c258ee6fd624f62f19ebab06 (MIT licence at the root, BSD-3 per file) in an isolated interpreter and stores the values, with the commit, the NumPy and SciPy versions and SciPy's constants, in tests/data/planck_pyradi_c7cfe9d.npz.

radiometry

Radiometric source functions, below every stage that emits or measures: the Planck radiance in wavenumber with its derivatives, the radiation constants and the per-µm and per-GHz views (radtrax#247), and brightness temperature, monochromatic and through an instrument's bands (radtrax#248).

surface/ (emission), solvers/ (layer sources), instrument/ (NEΔT) and the retrievals all reach it; it depends only on atmosphere.gas and exceptions, and names Instrument for type checkers only (radtrax#246).

planck

Planck spectral radiance in wavenumber, its derivatives and its unit views.

The black-body spectral radiance per unit wavenumber is

\[ B(\tilde\nu, T) = \frac{c_1 \tilde\nu^3}{e^{x} - 1} = \frac{c_1 \tilde\nu^3 e^{-x}}{1 - e^{-x}}, \qquad x = \frac{c_2 \tilde\nu}{T}, \]

in W m⁻² sr⁻¹ (cm⁻¹)⁻¹, with \(\tilde\nu\) in cm⁻¹ and \(T\) in K. Radiance, not exitance (\(M = \pi B\)), crosses every public boundary. Wavenumber is the one internal axis; wavelength and frequency enter only through the views to_per_um, from_per_um, to_per_ghz and from_per_ghz (#246, D8).

The radiation constants are computed at import from the SI 2019 exact \(h\), \(c\) and \(k\) (BOLTZMANN_J_PER_K lives in radtrax.atmosphere.gas):

• \(c_1 = 2hc^2 = 1.1910429723971885 \times 10^{-8}\) W m⁻² sr⁻¹ cm⁴ (C1_W_M2_SR_CM4),
• \(c_2 = hc/k = 1.4387768775039338\) cm K (C2_CM_K),
• \(\sigma = 2\pi^5 k^4 / (15 h^3 c^2) = 5.670374419 \times 10^{-8}\) W m⁻² K⁻⁴ (STEFAN_BOLTZMANN_W_M2_K4), so that \(\pi \int_0^\infty B \, d\tilde\nu = \sigma T^4\).

Numerics. The form with \(e^{-x}\) and \(1 - e^{-x} = -\mathrm{expm1}(-x)\) cannot overflow: over 0.4–15 µm and 150–2500 K, \(x\) runs from 0.384 to 239.8, and \(e^{x}\) overflows float32 above \(x \approx 88\). \(e^{-x}\) underflows gracefully to 0 (in float32, \(B\) at 0.4 µm is 0 below about 360 K), where log_planck_radiance stays finite. Both functions carry their closed-form derivatives through jax.custom_jvp, written in jnp operations so that second derivatives work:

\[ \frac{\partial B}{\partial T} = B \frac{x}{T (1 - e^{-x})}, \qquad \frac{\partial B}{\partial \tilde\nu} = B \left( \frac{3}{\tilde\nu} - \frac{c_2}{T (1 - e^{-x})} \right), \]

and the same brackets for \(\ln B\). Limits: \(B \to c_1 \tilde\nu^2 T / c_2\) as \(x \to 0\) (Rayleigh–Jeans, the passive-microwave end) and \(B \to c_1 \tilde\nu^3 e^{-x}\) as \(x \to \infty\) (Wien).

Inputs must satisfy \(\tilde\nu > 0\) and \(T > 0\) (at \(\tilde\nu = 0\) the form is \(0 \cdot \infty\)); nothing checks it, so the functions stay usable inside jit over a swath. Every function keeps the caller's floating dtype: the constants are Python floats (weakly typed), so float32 inputs stay float32 under jax_enable_x64.

planck_radiance

planck_radiance(nu_cm: Float[ArrayLike, '*#shape'], t_k: Float[ArrayLike, '*#shape']) -> Float[Array, '*#shape']

Planck spectral radiance per unit wavenumber [W m⁻² sr⁻¹ (cm⁻¹)⁻¹].

\[ B(\tilde\nu, T) = \frac{c_1 \tilde\nu^3 e^{-x}}{1 - e^{-x}}, \qquad x = \frac{c_2 \tilde\nu}{T}, \]

evaluated with expm1 so that nothing overflows; \(B\) underflows to 0 where \(e^{-x}\) does. The derivatives are the closed forms \(\partial B / \partial T = B x / (T (1 - e^{-x}))\) and \(\partial B / \partial \tilde\nu = B (3 / \tilde\nu - c_2 / (T (1 - e^{-x})))\) (jax.custom_jvp), differentiable again.

Parameters:

Name	Type	Description	Default
nu_cm	Float[ArrayLike, '*#shape']	Wavenumber \(\tilde\nu\) [cm⁻¹], > 0.	required
t_k	Float[ArrayLike, '*#shape']	Temperature \(T\) [K], > 0; broadcasts against nu_cm, so a (G,) grid and (L, 1) layer temperatures give (L, G).	required

Returns:

Type	Description
Float[Array, '*#shape']	\(B\), the broadcast shape of the inputs, in their floating dtype.

Examples:

>>> import jax
>>> nu = 1e4 / 10.9  # 10.9 µm
>>> f"{float(planck_radiance(nu, 300.0)):.6e}"
'1.143269e-01'
>>> f"{float(jax.grad(planck_radiance, argnums=1)(nu, 300.0)):.6e}"
'1.697611e-03'

Source code in src/radtrax/radiometry/planck.py

def planck_radiance(
    nu_cm: Float[ArrayLike, "*#shape"], t_k: Float[ArrayLike, "*#shape"]
) -> Float[Array, "*#shape"]:
    """Planck spectral radiance per unit wavenumber [W m⁻² sr⁻¹ (cm⁻¹)⁻¹].

    $$
    B(\\tilde\\nu, T) = \\frac{c_1 \\tilde\\nu^3 e^{-x}}{1 - e^{-x}},
    \\qquad x = \\frac{c_2 \\tilde\\nu}{T},
    $$

    evaluated with `expm1` so that nothing overflows; $B$ underflows to 0
    where $e^{-x}$ does. The derivatives are the closed forms
    $\\partial B / \\partial T = B x / (T (1 - e^{-x}))$ and
    $\\partial B / \\partial \\tilde\\nu = B (3 / \\tilde\\nu - c_2 / (T (1 - e^{-x})))$
    (`jax.custom_jvp`), differentiable again.

    Args:
        nu_cm: Wavenumber $\\tilde\\nu$ [cm⁻¹], > 0.
        t_k: Temperature $T$ [K], > 0; broadcasts against ``nu_cm``, so a
            ``(G,)`` grid and ``(L, 1)`` layer temperatures give ``(L, G)``.

    Returns:
        $B$, the broadcast shape of the inputs, in their floating dtype.

    Examples:
        >>> import jax
        >>> nu = 1e4 / 10.9  # 10.9 µm
        >>> f"{float(planck_radiance(nu, 300.0)):.6e}"
        '1.143269e-01'
        >>> f"{float(jax.grad(planck_radiance, argnums=1)(nu, 300.0)):.6e}"
        '1.697611e-03'
    """
    return _planck(*_float(nu_cm, t_k))

log_planck_radiance

log_planck_radiance(nu_cm: Float[ArrayLike, '*#shape'], t_k: Float[ArrayLike, '*#shape']) -> Float[Array, '*#shape']

Natural log of the Planck radiance, finite wherever \(\tilde\nu, T > 0\).

\[ \ln B = \ln c_1 + 3 \ln \tilde\nu - x - \ln(1 - e^{-x}), \]

with \(B\) in W m⁻² sr⁻¹ (cm⁻¹)⁻¹. It stays finite where \(B\) underflows (the visible at low temperature, especially in float32), and so do its derivatives \(\partial \ln B / \partial T = x / (T (1 - e^{-x}))\) and \(\partial \ln B / \partial \tilde\nu = 3/\tilde\nu - c_2 / (T (1 - e^{-x}))\) (jax.custom_jvp).

Parameters:

Name	Type	Description	Default
nu_cm	Float[ArrayLike, '*#shape']	Wavenumber \(\tilde\nu\) [cm⁻¹], > 0.	required
t_k	Float[ArrayLike, '*#shape']	Temperature \(T\) [K], > 0; broadcasts against nu_cm.	required

Returns:

Type	Description
Float[Array, '*#shape']	\(\ln B\), the broadcast shape of the inputs, in their floating dtype.

Examples:

>>> import jax.numpy as jnp
>>> nu = jnp.float32(1e4 / 0.4)  # 0.4 µm, where B underflows float32
>>> float(planck_radiance(nu, jnp.float32(150.0)))
0.0
>>> round(float(log_planck_radiance(nu, jnp.float32(150.0))), 1)
-227.7

Source code in src/radtrax/radiometry/planck.py

def log_planck_radiance(
    nu_cm: Float[ArrayLike, "*#shape"], t_k: Float[ArrayLike, "*#shape"]
) -> Float[Array, "*#shape"]:
    """Natural log of the Planck radiance, finite wherever $\\tilde\\nu, T > 0$.

    $$
    \\ln B = \\ln c_1 + 3 \\ln \\tilde\\nu - x - \\ln(1 - e^{-x}),
    $$

    with $B$ in W m⁻² sr⁻¹ (cm⁻¹)⁻¹. It stays finite where $B$ underflows
    (the visible at low temperature, especially in float32), and so do its
    derivatives $\\partial \\ln B / \\partial T = x / (T (1 - e^{-x}))$ and
    $\\partial \\ln B / \\partial \\tilde\\nu = 3/\\tilde\\nu - c_2 / (T (1 - e^{-x}))$
    (`jax.custom_jvp`).

    Args:
        nu_cm: Wavenumber $\\tilde\\nu$ [cm⁻¹], > 0.
        t_k: Temperature $T$ [K], > 0; broadcasts against ``nu_cm``.

    Returns:
        $\\ln B$, the broadcast shape of the inputs, in their floating dtype.

    Examples:
        >>> import jax.numpy as jnp
        >>> nu = jnp.float32(1e4 / 0.4)  # 0.4 µm, where B underflows float32
        >>> float(planck_radiance(nu, jnp.float32(150.0)))
        0.0
        >>> round(float(log_planck_radiance(nu, jnp.float32(150.0))), 1)
        -227.7
    """
    return _log_planck(*_float(nu_cm, t_k))

to_per_um

to_per_um(radiance: Float[ArrayLike, '*#shape'], nu_cm: Float[ArrayLike, '*#shape']) -> Float[Array, '*#shape']

A spectral density per cm⁻¹ as one per µm: \(L_\lambda = L_{\tilde\nu} \tilde\nu^2 / 10^4\).

From \(\lambda = 10^4 / \tilde\nu\) (µm), \(|d\tilde\nu / d\lambda| = \tilde\nu^2 / 10^4\); the factor is positive, so the conversion is the Jacobian's magnitude.

Parameters:

Name	Type	Description	Default
radiance	Float[ArrayLike, '*#shape']	\(L_{\tilde\nu}\) [W m⁻² sr⁻¹ (cm⁻¹)⁻¹].	required
nu_cm	Float[ArrayLike, '*#shape']	Wavenumber \(\tilde\nu\) [cm⁻¹] at which it is given.	required

Returns:

Type	Description
Float[Array, '*#shape']	\(L_\lambda\) [W m⁻² sr⁻¹ µm⁻¹], the broadcast shape of the inputs.

Examples:

>>> nu = 1e4 / 10.9
>>> f"{float(to_per_um(planck_radiance(nu, 300.0), nu)):.4f}"
'9.6227'

Source code in src/radtrax/radiometry/planck.py

def to_per_um(
    radiance: Float[ArrayLike, "*#shape"], nu_cm: Float[ArrayLike, "*#shape"]
) -> Float[Array, "*#shape"]:
    """A spectral density per cm⁻¹ as one per µm: $L_\\lambda = L_{\\tilde\\nu} \\tilde\\nu^2 / 10^4$.

    From $\\lambda = 10^4 / \\tilde\\nu$ (µm), $|d\\tilde\\nu / d\\lambda| = \\tilde\\nu^2 / 10^4$;
    the factor is positive, so the conversion is the Jacobian's magnitude.

    Args:
        radiance: $L_{\\tilde\\nu}$ [W m⁻² sr⁻¹ (cm⁻¹)⁻¹].
        nu_cm: Wavenumber $\\tilde\\nu$ [cm⁻¹] at which it is given.

    Returns:
        $L_\\lambda$ [W m⁻² sr⁻¹ µm⁻¹], the broadcast shape of the inputs.

    Examples:
        >>> nu = 1e4 / 10.9
        >>> f"{float(to_per_um(planck_radiance(nu, 300.0), nu)):.4f}"
        '9.6227'
    """
    return jnp.asarray(radiance) * jnp.asarray(nu_cm) ** 2 / _CM2_PER_UM

from_per_um

from_per_um(radiance_um: Float[ArrayLike, '*#shape'], nu_cm: Float[ArrayLike, '*#shape']) -> Float[Array, '*#shape']

A spectral density per µm as one per cm⁻¹: \(L_{\tilde\nu} = 10^4 L_\lambda / \tilde\nu^2\).

The inverse of to_per_um.

Parameters:

Name	Type	Description	Default
radiance_um	Float[ArrayLike, '*#shape']	\(L_\lambda\) [W m⁻² sr⁻¹ µm⁻¹].	required
nu_cm	Float[ArrayLike, '*#shape']	Wavenumber \(\tilde\nu\) [cm⁻¹] at which it is given, > 0.	required

Returns:

Type	Description
Float[Array, '*#shape']	\(L_{\tilde\nu}\) [W m⁻² sr⁻¹ (cm⁻¹)⁻¹], the broadcast shape of the inputs.

Examples:

>>> f"{float(from_per_um(9.622663, 1e4 / 10.9)):.5e}"
'1.14327e-01'

Source code in src/radtrax/radiometry/planck.py

def from_per_um(
    radiance_um: Float[ArrayLike, "*#shape"], nu_cm: Float[ArrayLike, "*#shape"]
) -> Float[Array, "*#shape"]:
    """A spectral density per µm as one per cm⁻¹: $L_{\\tilde\\nu} = 10^4 L_\\lambda / \\tilde\\nu^2$.

    The inverse of `to_per_um`.

    Args:
        radiance_um: $L_\\lambda$ [W m⁻² sr⁻¹ µm⁻¹].
        nu_cm: Wavenumber $\\tilde\\nu$ [cm⁻¹] at which it is given, > 0.

    Returns:
        $L_{\\tilde\\nu}$ [W m⁻² sr⁻¹ (cm⁻¹)⁻¹], the broadcast shape of the inputs.

    Examples:
        >>> f"{float(from_per_um(9.622663, 1e4 / 10.9)):.5e}"
        '1.14327e-01'
    """
    return jnp.asarray(radiance_um) * _CM2_PER_UM / jnp.asarray(nu_cm) ** 2

to_per_ghz

to_per_ghz(radiance: Float[ArrayLike, '*shape']) -> Float[Array, '*shape']

A spectral density per cm⁻¹ as one per GHz: \(L_f = L_{\tilde\nu} / c\).

With \(f = c \tilde\nu\) and \(c = 29.9792458\) GHz cm, \(d\tilde\nu / df = 1/c\).

Parameters:

Name	Type	Description	Default
radiance	Float[ArrayLike, '*shape']	\(L_{\tilde\nu}\) [W m⁻² sr⁻¹ (cm⁻¹)⁻¹].	required

Returns:

Type	Description
Float[Array, '*shape']	\(L_f\) [W m⁻² sr⁻¹ GHz⁻¹], the shape of radiance.

Examples:

>>> f"{float(to_per_ghz(29.9792458)):.12f}"
'1.000000000000'

Source code in src/radtrax/radiometry/planck.py

def to_per_ghz(radiance: Float[ArrayLike, "*shape"]) -> Float[Array, "*shape"]:
    """A spectral density per cm⁻¹ as one per GHz: $L_f = L_{\\tilde\\nu} / c$.

    With $f = c \\tilde\\nu$ and $c = 29.9792458$ GHz cm, $d\\tilde\\nu / df = 1/c$.

    Args:
        radiance: $L_{\\tilde\\nu}$ [W m⁻² sr⁻¹ (cm⁻¹)⁻¹].

    Returns:
        $L_f$ [W m⁻² sr⁻¹ GHz⁻¹], the shape of ``radiance``.

    Examples:
        >>> f"{float(to_per_ghz(29.9792458)):.12f}"
        '1.000000000000'
    """
    return jnp.asarray(radiance) / _C_GHZ_CM

from_per_ghz

from_per_ghz(radiance_ghz: Float[ArrayLike, '*shape']) -> Float[Array, '*shape']

A spectral density per GHz as one per cm⁻¹: \(L_{\tilde\nu} = c L_f\).

The inverse of to_per_ghz, with \(c = 29.9792458\) GHz cm.

Parameters:

Name	Type	Description	Default
radiance_ghz	Float[ArrayLike, '*shape']	\(L_f\) [W m⁻² sr⁻¹ GHz⁻¹].	required

Returns:

Type	Description
Float[Array, '*shape']	\(L_{\tilde\nu}\) [W m⁻² sr⁻¹ (cm⁻¹)⁻¹], the shape of radiance_ghz.

Examples:

>>> f"{float(from_per_ghz(1.0)):.4f}"
'29.9792'

Source code in src/radtrax/radiometry/planck.py

def from_per_ghz(radiance_ghz: Float[ArrayLike, "*shape"]) -> Float[Array, "*shape"]:
    """A spectral density per GHz as one per cm⁻¹: $L_{\\tilde\\nu} = c L_f$.

    The inverse of `to_per_ghz`, with $c = 29.9792458$ GHz cm.

    Args:
        radiance_ghz: $L_f$ [W m⁻² sr⁻¹ GHz⁻¹].

    Returns:
        $L_{\\tilde\\nu}$ [W m⁻² sr⁻¹ (cm⁻¹)⁻¹], the shape of ``radiance_ghz``.

    Examples:
        >>> f"{float(from_per_ghz(1.0)):.4f}"
        '29.9792'
    """
    return jnp.asarray(radiance_ghz) * _C_GHZ_CM

brightness

Brightness temperature: Planck inverted, per wavenumber and per band.

The monochromatic brightness temperature is the closed-form inverse of planck_radiance,

\[ T_b(\tilde\nu, L) = \frac{c_2 \tilde\nu}{\ln(1 + c_1 \tilde\nu^3 / L)}, \]

evaluated through logs, \(\ln(1 + r) = \mathrm{logaddexp}(0, \ln r)\) with \(\ln r = \ln c_1 + 3 \ln \tilde\nu - \ln L\), so that it stays finite in any precision; autodiff of that form gives \(\partial T_b / \partial L = 1 / (\partial B / \partial T)|_{T_b}\).

An instrument is linear in radiance, \(y = A L\) with \(A\) its response matrix (Instrument.response_matrix, (B, G)), so the band-integrated Planck radiance at one temperature per band is

\[ \bar B_b(T_b) = \sum_g A_{bg} B(\tilde\nu_g, T_b), \]

and the band brightness temperature solves \(\bar B_b(T_b) = y_b\) for every band. The bands decouple (the Jacobian of the residual is diagonal), so band_brightness_temperature runs optimistix.Newton with a diagonal linear solve, from the monochromatic \(T_b\) at each band's centroid \(\bar\nu_b = \sum_g A_{bg} \tilde\nu_g / \sum_g A_{bg}\), and is differentiated by the implicit function theorem: \(\partial T_b / \partial y_b = 1 / (\partial \bar B_b / \partial T_b)\) and \(\partial T_b / \partial \theta = -(\partial \bar B_b / \partial \theta) / (\partial \bar B_b / \partial T_b)\) for instrument parameters \(\theta\).

Every formula holds for any linear \(A\); the centroid divides by the row sum, so an instrument in other radiance units (rows that do not sum to 1) works in its own units. The dense \(A\) is materialised once per call, so band brightness temperature is for multispectral thermal bands (tens of bands); for a hyperspectral instrument use the monochromatic \(T_b\) at the centroid.

Shapes: G wavenumbers, B bands; *B means a scalar or (B,).

brightness_temperature

brightness_temperature(nu_cm: Float[ArrayLike, '*#shape'], radiance: Float[ArrayLike, '*#shape']) -> Float[Array, '*#shape']

Monochromatic brightness temperature [K], the inverse of planck_radiance.

\[ T_b = \frac{c_2 \tilde\nu}{\ln(1 + c_1 \tilde\nu^3 / L)}, \]

with \(\ln(1 + r)\) as \(\mathrm{logaddexp}(0, \ln r)\), finite and accurate from the Rayleigh–Jeans limit (\(r \to 0\)) to the Wien limit (\(r \to \infty\)). A non-positive radiance has no brightness temperature: it gives NaN, through a double where, so the gradient elsewhere stays finite.

Parameters:

Name	Type	Description	Default
nu_cm	Float[ArrayLike, '*#shape']	Wavenumber \(\tilde\nu\) [cm⁻¹], > 0.	required
radiance	Float[ArrayLike, '*#shape']	\(L\) [W m⁻² sr⁻¹ (cm⁻¹)⁻¹]; broadcasts against nu_cm.	required

Returns:

Type	Description
Float[Array, '*#shape']	\(T_b\) [K], the broadcast shape of the inputs; NaN where \(L \le 0\).

Examples:

>>> from radtrax.radiometry import planck_radiance
>>> nu = 1e4 / 10.9
>>> round(float(brightness_temperature(nu, planck_radiance(nu, 300.0))), 9)
300.0

Source code in src/radtrax/radiometry/brightness.py

def brightness_temperature(
    nu_cm: Float[ArrayLike, "*#shape"], radiance: Float[ArrayLike, "*#shape"]
) -> Float[Array, "*#shape"]:
    """Monochromatic brightness temperature [K], the inverse of `planck_radiance`.

    $$
    T_b = \\frac{c_2 \\tilde\\nu}{\\ln(1 + c_1 \\tilde\\nu^3 / L)},
    $$

    with $\\ln(1 + r)$ as $\\mathrm{logaddexp}(0, \\ln r)$, finite and accurate
    from the Rayleigh–Jeans limit ($r \\to 0$) to the Wien limit
    ($r \\to \\infty$). A non-positive radiance has no brightness temperature:
    it gives NaN, through a double where, so the gradient elsewhere stays
    finite.

    Args:
        nu_cm: Wavenumber $\\tilde\\nu$ [cm⁻¹], > 0.
        radiance: $L$ [W m⁻² sr⁻¹ (cm⁻¹)⁻¹]; broadcasts against ``nu_cm``.

    Returns:
        $T_b$ [K], the broadcast shape of the inputs; NaN where $L \\le 0$.

    Examples:
        >>> from radtrax.radiometry import planck_radiance
        >>> nu = 1e4 / 10.9
        >>> round(float(brightness_temperature(nu, planck_radiance(nu, 300.0))), 9)
        300.0
    """
    dtype = jnp.result_type(nu_cm, radiance, float)
    nu = jnp.asarray(nu_cm, dtype)
    radiance = jnp.asarray(radiance, dtype)
    positive = radiance > 0
    safe = jnp.where(positive, radiance, 1.0)
    log_r = _LOG_C1 + 3.0 * jnp.log(nu) - jnp.log(safe)
    t_b = C2_CM_K * nu / jnp.logaddexp(0.0, log_r)
    return jnp.where(positive, t_b, jnp.nan)

band_planck_radiance

band_planck_radiance(instrument: Instrument, t_k: Float[ArrayLike, '*B']) -> Float[Array, ' B']

Band-integrated Planck radiance \(\bar B_b(T_b) = \sum_g A_{bg} B(\tilde\nu_g, T_b)\).

One temperature for every band (a scalar) is one instrument call, instrument(planck_radiance(instrument.nu_cm, t_k)); one temperature per band, (B,), contracts the response matrix with \(B\) on the (B, G) grid. Differentiable in t_k and in the instrument.

Parameters:

Name	Type	Description	Default
instrument	Instrument	The bands, any Instrument (linear in radiance).	required
t_k	Float[ArrayLike, '*B']	Temperature [K], a scalar or one per band, (B,).	required

Returns:

Type	Description
Float[Array, ' B']	\(\bar B\) (B,), in the instrument's radiance units
Float[Array, ' B']	(W m⁻² sr⁻¹ (cm⁻¹)⁻¹ for a response whose rows sum to 1).

Examples:

>>> import jax.numpy as jnp
>>> from radtrax.instrument import GaussianInstrument
>>> nu = jnp.linspace(700.0, 1300.0, 1201)
>>> tir = GaussianInstrument(nu, jnp.array([10900.0]), jnp.array([500.0]))
>>> f"{float(band_planck_radiance(tir, 300.0)[0]):.4f}"
'0.1142'

Source code in src/radtrax/radiometry/brightness.py

def band_planck_radiance(
    instrument: Instrument, t_k: Float[ArrayLike, "*B"]
) -> Float[Array, " B"]:
    """Band-integrated Planck radiance $\\bar B_b(T_b) = \\sum_g A_{bg} B(\\tilde\\nu_g, T_b)$.

    One temperature for every band (a scalar) is one instrument call,
    ``instrument(planck_radiance(instrument.nu_cm, t_k))``; one temperature
    per band, ``(B,)``, contracts the response matrix with $B$ on the
    ``(B, G)`` grid. Differentiable in ``t_k`` and in the instrument.

    Args:
        instrument: The bands, any `Instrument` (linear in radiance).
        t_k: Temperature [K], a scalar or one per band, ``(B,)``.

    Returns:
        $\\bar B$ ``(B,)``, in the instrument's radiance units
        (W m⁻² sr⁻¹ (cm⁻¹)⁻¹ for a response whose rows sum to 1).

    Examples:
        >>> import jax.numpy as jnp
        >>> from radtrax.instrument import GaussianInstrument
        >>> nu = jnp.linspace(700.0, 1300.0, 1201)
        >>> tir = GaussianInstrument(nu, jnp.array([10900.0]), jnp.array([500.0]))
        >>> f"{float(band_planck_radiance(tir, 300.0)[0]):.4f}"
        '0.1142'
    """
    t = jnp.asarray(t_k)
    if t.ndim == 0:
        return instrument(planck_radiance(instrument.nu_cm, t))
    return _band_planck(instrument.response_matrix(), instrument.nu_cm, t)

band_brightness_temperature

band_brightness_temperature(instrument: Instrument, y: Float[ArrayLike, ' B'], *, rtol: float = 1e-12, atol: float = 1e-10, max_steps: int = 32, throw: bool = True) -> Float[Array, ' B']

Band brightness temperature: the \(T_b\) with \(\bar B_b(T_b) = y_b\) [K].

Solves \(\ln \bar B_b(T_b) - \ln y_b = 0\) for every band by optimistix.Newton, with a diagonal Jacobian (one JVP per step) and a diagonal linear solve, from the monochromatic \(T_b\) at each band's centroid. The residual is taken in log radiance so that the tolerances mean the same at 0.4 µm and at 15 µm; it has the same root and the same implicit derivatives as \(\bar B_b - y_b\). The unknown is \(u = \ln T_b\), so no Newton step can leave \(T_b > 0\): in \(T\) a broad band's first step from the centroid guess can overshoot below 0 K. Gradients come from the implicit function theorem (optimistix.ImplicitAdjoint), in y and in the instrument's parameters alike.

A band with \(y_b \le 0\) has no brightness temperature: it is solved on a safe value and returned as NaN (double where), so a masked band does not poison the gradient of the others.

Parameters:

Name	Type	Description	Default
instrument	Instrument	The bands, any Instrument (linear in radiance).	required
y	Float[ArrayLike, ' B']	Band values, (B,), in the instrument's radiance units.	required
rtol	float	Relative tolerance of the Newton solve in \(\ln T\).	1e-12
atol	float	Absolute tolerance of the Newton solve (in \(\ln T\) and in the log residual).	1e-10
max_steps	int	The most Newton steps.	32
throw	bool	Raise when the solve fails to converge (equinox.EquinoxRuntimeError under jit); otherwise return the last iterate.	True

Returns:

Type	Description
Float[Array, ' B']	\(T_b\) [K], (B,); NaN where \(y_b \le 0\).

Examples:

>>> import jax.numpy as jnp
>>> from radtrax.instrument import GaussianInstrument
>>> nu = jnp.linspace(700.0, 1300.0, 1201)
>>> tir = GaussianInstrument(nu, jnp.array([10900.0]), jnp.array([500.0]))
>>> y = band_planck_radiance(tir, 300.0)
>>> round(float(band_brightness_temperature(tir, y)[0]), 9)
300.0

Source code in src/radtrax/radiometry/brightness.py

def band_brightness_temperature(
    instrument: Instrument,
    y: Float[ArrayLike, " B"],
    *,
    rtol: float = 1e-12,
    atol: float = 1e-10,
    max_steps: int = 32,
    throw: bool = True,
) -> Float[Array, " B"]:
    """Band brightness temperature: the $T_b$ with $\\bar B_b(T_b) = y_b$ [K].

    Solves $\\ln \\bar B_b(T_b) - \\ln y_b = 0$ for every band by
    `optimistix.Newton`, with a diagonal Jacobian (one JVP per step) and a
    diagonal linear solve, from the monochromatic $T_b$ at each band's
    centroid. The residual is taken in log radiance so that the tolerances
    mean the same at 0.4 µm and at 15 µm; it has the same root and the same
    implicit derivatives as $\\bar B_b - y_b$. The unknown is $u = \\ln T_b$,
    so no Newton step can leave $T_b > 0$: in $T$ a broad band's first step
    from the centroid guess can overshoot below 0 K. Gradients come from the
    implicit function theorem (`optimistix.ImplicitAdjoint`), in ``y`` and in
    the instrument's parameters alike.

    A band with $y_b \\le 0$ has no brightness temperature: it is solved on a
    safe value and returned as NaN (double where), so a masked band does not
    poison the gradient of the others.

    Args:
        instrument: The bands, any `Instrument` (linear in radiance).
        y: Band values, ``(B,)``, in the instrument's radiance units.
        rtol: Relative tolerance of the Newton solve in $\\ln T$.
        atol: Absolute tolerance of the Newton solve (in $\\ln T$ and in the
            log residual).
        max_steps: The most Newton steps.
        throw: Raise when the solve fails to converge
            (``equinox.EquinoxRuntimeError`` under ``jit``); otherwise
            return the last iterate.

    Returns:
        $T_b$ [K], ``(B,)``; NaN where $y_b \\le 0$.

    Examples:
        >>> import jax.numpy as jnp
        >>> from radtrax.instrument import GaussianInstrument
        >>> nu = jnp.linspace(700.0, 1300.0, 1201)
        >>> tir = GaussianInstrument(nu, jnp.array([10900.0]), jnp.array([500.0]))
        >>> y = band_planck_radiance(tir, 300.0)
        >>> round(float(band_brightness_temperature(tir, y)[0]), 9)
        300.0
    """
    response = instrument.response_matrix()  # (B, G)
    nu = instrument.nu_cm  # (G,)
    y = jnp.asarray(y, dtype=response.dtype)
    row_sum = jnp.sum(response, axis=1)  # (B,)
    centroid = response @ nu / row_sum  # (B,)
    positive = y > 0
    y_safe = jnp.where(positive, y, row_sum * planck_radiance(centroid, 300.0))
    t0 = brightness_temperature(centroid, y_safe / row_sum)

    def residual(
        u: Float[Array, " B"],
        args: tuple[Float[Array, "B G"], Float[Array, " G"], Float[Array, " B"]],
    ) -> Float[Array, " B"]:
        a, grid, log_y = args
        return jnp.log(_band_planck(a, grid, jnp.exp(u))) - log_y

    solver = optx.Newton(rtol=rtol, atol=atol, linear_solver=lx.Diagonal())
    solution = optx.root_find(
        residual,
        solver,
        jnp.log(t0),
        args=(response, nu, jnp.log(y_safe)),
        max_steps=max_steps,
        throw=throw,
        tags=frozenset({lx.diagonal_tag}),
        adjoint=optx.ImplicitAdjoint(linear_solver=lx.Diagonal()),
    )
    return jnp.where(positive, jnp.exp(solution.value), jnp.nan)