"""Parametric isothermal transmission spectrum from the shared zeta grid.
zeta(lambda) = -ln P_photo(lambda)/bar is tabulated on a (T, [M/H], C/O) grid
for two fiducial planets (1 R_J at g = 10 m/s^2; 0.3 R_J at g = 30 m/s^2),
i.e. at two nodes x = 0 and -1.15 of the slant-column axis x = 0.5 ln[(R_p H)/(R_p H)_0]
(precompute/make_zeta_grid.py). For any other planet, with
x_p = 0.5 ln[(R_p H)/(R_p,0 H_0)] - ln(T_eq/T_chem)::
s(lam) = (zeta_1 - zeta_0) / (x_1 - x_0) slant-column exponent, ~1/(1+a)
zeta_p = zeta_0 + s (x_p - x_0) linear in x, extrapolated
+ 1.5 s [ln(r/R_p) - ln(r_0/R_fid)] local column at r(lam), iterated
zeta_g1 = zeta_p - s gamma_E gas tau_chord = 1 level
zeta_h1 = -ln P_haze + 1.5 ln(r/R_p) haze tau_chord = 1 level
tau(z) = exp((zeta_g1 - z)/s) + exp(zeta_h1 - z)
zeta = zeta_d + int_0^inf [1 - exp(-tau(zeta_d + u))] du
R(lam) = R_p / (1 - (H / R_p) (zeta - median_band(zeta))), depth = (R / R_s)^2
s < 1 where the opacity grows with pressure (tau ~ P^(1+a): CIA, Lorentz
wings), so both the column scaling and the vertical tau scale height (s H)
shrink. The grid zeta is an effective-radius level (POSEIDON depth -> R ->
r(P)), which for an exponential tau profile sits s gamma_E above the
tau_chord = 1 level; the haze level is a tau_chord = 1 pressure with exponent
1. Both are therefore taken back to tau = 1 and the thin-atmosphere effective
radius of the summed profile above an opaque deck at zeta_d = -ln P_cloud
(no deck: zeta_d deep enough that tau(zeta_d) >= e^12) is integrated
numerically. For pure gas this is zeta_d + s Ein(exp((zeta_p - zeta_d)/s
- gamma_E)) (Ein the entire exponential integral; the no-haze fast path),
which tends to max(zeta_p, zeta_d) in both limits but keeps the partial
absorption of gas whose tau = 1 level lies just below the deck; the integral
also gives the haze its own tau scale height (H, not s H) where it dominates.
so R_p is the band-median transit radius and H = k T_eq / (mu m_H g(R_p)).
The radius mapping is the exact isothermal hydrostatic solution with
g ~ 1/r^2 (1/r is linear in ln P), which POSEIDON also integrates; its first-
order form R_p + H dzeta under-predicts puffy atmospheres (H/R_p ~ 0.03 gives
~30 % too little amplitude). The chord column sqrt(2 pi r H(r)) ~ r^1.5 is
likewise evaluated at each wavelength's own radius r(lam) (fixed point in
ln r, 2-4 steps), net of the curvature the fiducial already carries at its
own r_0(lam); without it a 1 R_J, g = 1.7 m/s^2 planet (H/R_p = 0.033) came
out 7 % low in amplitude. The residual 0.1-0.2 H there is the second-order
sphericity term (tau sampled across r(z), not at r_eff). Fidelity: isothermal,
equilibrium chemistry, gray deck, power-law haze; no pressure-broadening
correction.
"""
import numpy as np
from . import tables
K_B = 1.380649e-23
M_H = 1.6735575e-27
G = 6.674e-11
R_J = 7.1492e7
M_J = 1.898e27
R_SUN = 6.957e8
SIGMA_H2_RAYLEIGH_350NM = 5.31e-31 # m^2, POSEIDON MacMad17 haze anchor
HAZE_LAMBDA0_UM = 0.35
BAND_UM = (0.40, 1.70)
EULER_GAMMA = 0.5772156649015329
S_RANGE = (0.25, 1.5) # physical range of the slant-column exponent
# Ein(a) = int_0^a (1 - e^-s)/s ds, tabulated once on a log grid (numpy only);
# Ein(a) = gamma + ln a + E1(a) with E1(50) ~ 5e-24, so a > 50 is analytic.
_EIN_LOG_A = np.linspace(np.log(1e-8), np.log(50.0), 6001)
_EIN_A = np.exp(_EIN_LOG_A)
_EIN_TABLE = np.concatenate([[_EIN_A[0]], _EIN_A[0] + np.cumsum(
0.5 * (-np.expm1(-_EIN_A[1:]) / _EIN_A[1:] + -np.expm1(-_EIN_A[:-1]) / _EIN_A[:-1])
* np.diff(_EIN_A))])
def _ein(a):
"""Entire exponential integral Ein(a) for a >= 0 (vectorised)."""
a = np.asarray(a, float)
out = np.interp(np.log(np.clip(a, 1e-300, None)), _EIN_LOG_A, _EIN_TABLE)
out = np.where(a < _EIN_A[0], a, out)
return np.where(a > _EIN_A[-1], EULER_GAMMA + np.log(np.clip(a, 1e-300, None)), out)
def _zeta_eff(zeta_g1, s, zeta_h1, zeta_d=None, du=0.02, pad=16.0, chunk=2048):
"""Effective-radius zeta of tau(z) = exp((zeta_g1 - z)/s) + exp(zeta_h1 - z)
above an opaque deck at zeta_d: zeta_d + int_0^inf [1 - exp(-tau)] du.
Trapezoid in u with the Euler-Maclaurin end correction (O(du^4); 1e-6
agreement with s Ein(...) for pure gas), up to where both components have
tau < e^-pad. Without a deck the bottom is placed where either component
has tau = e^12, which the result does not depend on.
"""
zeta_g1, s, zeta_h1 = (np.asarray(np.atleast_1d(v), float) for v in (zeta_g1, s, zeta_h1))
if zeta_d is None:
zeta_d = np.maximum(zeta_g1 - 12.0 * s, zeta_h1 - 12.0)
zeta_g1, s, zeta_h1, zeta_d = np.broadcast_arrays(zeta_g1, s, zeta_h1, zeta_d)
top = np.maximum(zeta_g1 + pad * s, zeta_h1 + pad)
u = np.arange(0.0, max(float((top - zeta_d).max()), 0.0) + du, du)
out = np.empty(zeta_g1.shape)
for i in range(0, len(out), chunk): # chunked: ~30 MB per temporary
sl = slice(i, i + chunk)
z = zeta_d[sl, None] + u[None, :]
tau_g = np.exp(np.minimum((zeta_g1[sl, None] - z) / s[sl, None], 350.0))
tau_h = np.exp(np.minimum(zeta_h1[sl, None] - z, 350.0))
f = -np.expm1(-(tau_g + tau_h))
df0 = -np.exp(-(tau_g[:, 0] + tau_h[:, 0])) * (tau_g[:, 0] / s[sl] + tau_h[:, 0])
out[sl] = (zeta_d[sl] + du * (f.sum(axis=1) - 0.5 * (f[:, 0] + f[:, -1]))
+ du ** 2 / 12.0 * df0)
return out
[docs]
def grid_wavelengths(tables_dir=None):
return tables.load_zeta_grid(tables_dir)["wavelength_um"]
def _interp_weights(axis, x):
"""Lower index and fractional weight for linear interpolation on `axis`."""
axis = np.asarray(axis, float)
if len(axis) == 1:
return 0, 0.0
i = int(np.clip(np.searchsorted(axis, x, side="right") - 1, 0, len(axis) - 2))
return i, float((x - axis[i]) / (axis[i + 1] - axis[i]))
def _trilinear(values, axes, point):
"""Trilinear interpolation of values[a, b, c, ...] at point = (x, y, z)."""
(i0, t0), (i1, t1), (i2, t2) = [_interp_weights(a, p)
for a, p in zip(axes, point)]
out = 0.0
for di, w0 in ((0, 1 - t0), (1, t0)):
for dj, w1 in ((0, 1 - t1), (1, t1)):
for dk, w2 in ((0, 1 - t2), (1, t2)):
w = w0 * w1 * w2
if w:
corner = values[min(i0 + di, values.shape[0] - 1),
min(i1 + dj, values.shape[1] - 1),
min(i2 + dk, values.shape[2] - 1)]
out = out + w * np.asarray(corner, dtype=float)
return out
def _check_range(name, x, axis):
lo, hi = float(np.min(axis)), float(np.max(axis))
if not (lo <= x <= hi):
raise ValueError(f"{name}={x} outside the grid range [{lo}, {hi}]")
[docs]
def transmission_depth(wl_um, *, T_eq_K, R_p_RJ, M_p_MJ, R_s_RSun,
log_met=0.0, c_to_o=0.55, mu=None,
log_P_cloud_bar=None, log_a_haze=None, haze_slope=-4.0,
tables_dir=None):
"""Transit depth on `wl_um` plus a dict of derived quantities.
R_p_RJ is the band-median (0.40-1.70 um) transit radius. T_eq below the
grid floor clamps the chemistry lookup only (reported as T_clamped); the
scale height always uses the true T_eq. mu=None takes the equilibrium-
chemistry mean molecular weight from the grid.
"""
g = tables.load_zeta_grid(tables_dir)
wl_grid = g["wavelength_um"]
wl_um = np.atleast_1d(np.asarray(wl_um, float))
for name, x in (("T_eq_K", T_eq_K), ("R_p_RJ", R_p_RJ),
("M_p_MJ", M_p_MJ), ("R_s_RSun", R_s_RSun)):
if not (np.isfinite(x) and x > 0):
raise ValueError(f"{name} must be > 0, got {x}")
if mu is not None and not (np.isfinite(mu) and mu > 0):
raise ValueError(f"mu must be > 0, got {mu}")
# Stored zeta are bin averages at bin centres and the last partial bin of
# the nominal span is dropped by the binning, so accept requests up to
# 1.5 bin widths beyond the outer centres (half a bin past the outer bin
# edges, which covers the nominal 0.35-1.85 um span); np.interp holds the
# edge value there.
lo = wl_grid[0] - 1.5 * (wl_grid[1] - wl_grid[0])
hi = wl_grid[-1] + 1.5 * (wl_grid[-1] - wl_grid[-2])
if wl_um.min() < lo or wl_um.max() > hi:
raise ValueError(
f"wavelength range [{wl_um.min():.3f}, {wl_um.max():.3f}] um "
f"outside the grid [{wl_grid.min():.3f}, {wl_grid.max():.3f}] um")
_check_range("log_met", log_met, g["log_met"])
_check_range("c_to_o", c_to_o, g["c_to_o"])
T_chem = float(np.clip(T_eq_K, g["T_K"].min(), g["T_K"].max()))
T_clamped = T_chem != float(T_eq_K)
axes = (g["T_K"], g["log_met"], g["c_to_o"])
point = (T_chem, log_met, c_to_o)
mu_chem = float(_trilinear(g["mu"], axes, point))
mu_used = mu_chem if mu is None else float(mu)
fids = g["meta"]["fiducials"]
fid = fids[0]
R_p = R_p_RJ * R_J
R_s = R_s_RSun * R_SUN
g_p = G * M_p_MJ * M_J / R_p ** 2
H = K_B * T_eq_K / (mu_used * M_H * g_p)
H_fid = K_B * T_chem / (mu_chem * M_H * fid["g_ms2"])
shift = (0.5 * np.log((R_p * H) / (fid["R_p_m"] * H_fid))
- np.log(T_eq_K / T_chem))
# Linear in x through the two fiducial nodes (extrapolated outside them);
# the slope is the slant-column exponent, also used for the deck tail.
x_nodes = np.array([f["x_column"] for f in fids], float)
zeta_0 = _trilinear(g["zeta"][..., 0, :], axes, point)
zeta_1 = _trilinear(g["zeta"][..., 1, :], axes, point)
slope = np.clip((zeta_1 - zeta_0) / (x_nodes[1] - x_nodes[0]), *S_RANGE)
zeta_p = zeta_0 + slope * (shift - x_nodes[0])
zeta_d = None if log_P_cloud_bar is None else -log_P_cloud_bar * np.log(10.0)
zeta_h1 = None
if log_a_haze is not None:
kappa = 10.0 ** log_a_haze * SIGMA_H2_RAYLEIGH_350NM * (
wl_grid / HAZE_LAMBDA0_UM) ** haze_slope # m^2 per molecule
P_haze_bar = K_B * T_eq_K / (kappa * np.sqrt(2 * np.pi * R_p * H)) / 1e5
zeta_h1 = -np.log(P_haze_bar)
def zeta_eff(dz_gas, dz_haze):
"""Effective-radius zeta with the local-column corrections applied to the
gas (effective-radius level zeta_p) and haze (tau_chord = 1 level)."""
zg = zeta_p + dz_gas
if zeta_h1 is not None:
# grid zeta is an effective-radius level (tau_chord = 1 + s gamma_E);
# the haze level is tau_chord = 1 with exponent 1: sum the two tau
# profiles at their tau = 1 levels and integrate the effective radius
return _zeta_eff(zg - slope * EULER_GAMMA, slope, zeta_h1 + dz_haze, zeta_d)
if zeta_d is not None:
with np.errstate(over="ignore"): # _ein is analytic for large a
return zeta_d + slope * _ein(np.exp((zg - zeta_d) / slope - EULER_GAMMA))
return zg
band = (wl_grid >= BAND_UM[0]) & (wl_grid <= BAND_UM[1])
# Local slant column: the chord column at wavelength lambda is
# sqrt(2 pi r H(r)) ~ r^1.5 at that wavelength's own radius r(lambda), not
# at R_p, so zeta_g1 gains 1.5 s ln(r/R_p) and zeta_h1 gains 1.5 ln(r/R_p).
# The grid already carries the fiducial's own curvature, 1.5 s ln(r_0/R_fid)
# with r_0 = r(zeta_0) through the same mapping, which is subtracted so it
# is not counted twice (a planet on the node gets no correction; constant
# offsets drop out in the band median). Fixed point in ln r: the
# contraction factor is ~1.5 s H/R_p, so thin planets converge in 2 steps.
ln_r0 = -np.log1p(-(H_fid / fid["R_p_m"]) * (zeta_0 - np.median(zeta_0[band])))
ln_r = ln_r0 # start from the fiducial's shape
for n_iter in range(1, 9):
dz_gas, dz_haze = 1.5 * slope * (ln_r - ln_r0), 1.5 * ln_r
zeta = zeta_eff(dz_gas, dz_haze)
zeta_ref = float(np.median(zeta[band]))
# 1/r linear in ln P (isothermal, g ~ r^-2); the denominator only
# approaches 0 for an unbound atmosphere (H/R_p * dzeta -> 1).
denom = 1.0 - (H / R_p) * (zeta - zeta_ref)
if denom.min() <= 0.1:
raise ValueError(
f"atmosphere is unbound for this (T_eq, mu, g): H/R_p = {H / R_p:.3f}")
ln_r, ln_r_prev = -np.log(denom), ln_r
if np.abs(ln_r - ln_r_prev).max() < 1e-4:
break
radius = R_p / denom
depth_grid = (radius / R_s) ** 2
depth = np.interp(wl_um, wl_grid, depth_grid)
info = {"g_ms2": float(g_p), "H_m": float(H), "H_over_Rp": float(H / R_p),
"mu_used": mu_used,
"mu_source": "chemistry" if mu is None else "user",
"T_chem_K": T_chem, "T_clamped": bool(T_clamped),
"zeta_ref": zeta_ref, "zeta_shift": float(shift),
"column_slope_median": float(np.median(slope[band])),
"column_iterations": n_iter,
"column_correction_max_H": float(max(
np.abs(dz_gas).max(), 0.0 if zeta_h1 is None else np.abs(dz_haze).max())),
"feature_amplitude_ppm": float(np.ptp(depth_grid[band]) * 1e6),
"band_um": list(BAND_UM)}
return depth, info