Spectral Contrast Models#
The default multi-band GP uses the Planck function ratio \(f_{\rm spot}(\lambda) = B_\lambda(T_{\rm spot})/B_\lambda(T_{\rm phot})\) evaluated at a single wavelength per band. This is a good approximation for broadband continuum, but ignores molecular absorption features (e.g., TiO bands below ~4000 K) and the finite width of real filter bandpasses.
The spectral module provides SpectralContrastModel, which precomputes
bandpass-integrated contrast ratios from model-atmosphere spectra,
giving a more realistic mapping from \(T_{\rm spot}\) to \(c(\lambda)\).
What you will learn
Build a
SpectralContrastModelfrom a spectrum provider and bandpass setVisualize the precomputed contrast table and compare to the Planck approximation
Plug the spectral model into
MultiBandGPSolveras a drop-in replacementUse
KorgProviderfor 1D LTE model-atmosphere spectraUse
pyphotfor real instrument filter curves
Prerequisites
This tutorial builds on the Multi-Band GP Tutorial.
Run that notebook first to understand the chromatic kernel and MultiBandGPSolver.
1. Setup: simulate multi-band data#
We quickly reproduce the simulated 3-band dataset from the Multi-Band GP Tutorial.
import numpy as np
import matplotlib.pyplot as plt
import jax.numpy as jnp
from spotgp import (
TrapezoidSymmetricEnvelope,
VisibilityFunction,
SpotEvolutionModel,
LightcurveModel,
)
from spotgp.contrast import contrast_factor
from spotgp.multiband import MultiBandData, MultiBandGPSolver
from spotgp.spectral import BlackbodyProvider, BandpassSet, SpectralContrastModel
np.random.seed(42)
# Geometric model and lightcurve
envelope = TrapezoidSymmetricEnvelope(lspot=8.0, tau_spot=3.0)
visibility = VisibilityFunction(peq=5.0, kappa=0.2, inc=np.pi / 3)
model = SpotEvolutionModel(envelope=envelope, visibility=visibility, sigma_k=0.01)
lc = LightcurveModel.from_spot_model(spot_model=model, nspot=25, tsim=60, tsamp=0.5)
delta_F = lc.flux - np.median(lc.flux)
t = lc.t
# Temperatures and bands
T_phot_true = 5800.0
T_spot_true = 4200.0
band_info = {
"SDSS_g": {"wavelength": 4640.0, "noise": 0.0015},
"Kepler": {"wavelength": 6400.0, "noise": 0.0008},
"TESS": {"wavelength": 7865.0, "noise": 0.0012},
}
colors = {"SDSS_g": "#4477AA", "Kepler": "#228833", "TESS": "#CC3311"}
bands = {}
for name, info in band_info.items():
c = float(contrast_factor(jnp.array(info["wavelength"]), T_spot_true, T_phot_true))
flux_band = 1.0 + c * delta_F + np.random.normal(0, info["noise"], len(t))
bands[name] = {
"x": t, "y": flux_band,
"yerr": np.full(len(t), info["noise"]),
"wavelength": info["wavelength"],
}
data = MultiBandData(bands)
hparam = dict(peq=5.0, kappa=0.2, inc=np.pi / 3, lspot=8.0, tau_spot=3.0, sigma_k=0.01)
bounds = {
"peq": (2.0, 15.0), "kappa": (0.001, 0.8), "inc": (0.1, np.pi - 0.1),
"lspot": (1.0, 20.0), "tau_spot": (0.1, 8.0),
"log_sigma_k": (-4.0, -0.5), "T_spot": (2500.0, 5500.0),
}
print(f"Data: {data.N} observations across {data.n_bands} bands")
2. Build a SpectralContrastModel#
The spectral contrast model has three components:
Spectrum provider – generates \(F_\lambda(T_{\rm eff})\). Two built-in options:
BlackbodyProvider: Planck function (no dependencies, fast)KorgProvider: 1D LTE spectral synthesis via Korg.jl (realistic, requires Julia)
Bandpass set – defines the photometric filters. Accepts either:
Custom arrays:
{"wavelength": ..., "throughput": ...}pyphot filter names:
"TESS","KEPLER_Kp","SDSS_g", etc.
Contrast table – precomputed \(f_{\rm spot}(T_{\rm spot},\, {\rm band})\) on a temperature grid, interpolated at inference time with JAX-traceable
jnp.interp.
Let’s start with BlackbodyProvider and custom Gaussian bandpasses.
# 1. Spectrum provider
provider = BlackbodyProvider(wl_range=(3000, 11000), n_points=3000)
# 2. Define bandpasses (Gaussian approximations to real filters)
wl_g = np.linspace(3800, 5500, 200)
tp_g = np.exp(-0.5 * ((wl_g - 4640) / 350) ** 2)
wl_kep = np.linspace(4200, 9000, 300)
tp_kep = np.exp(-0.5 * ((wl_kep - 6400) / 800) ** 2)
wl_tess = np.linspace(6000, 10500, 300)
tp_tess = np.exp(-0.5 * ((wl_tess - 7865) / 900) ** 2)
bps = BandpassSet({
"SDSS_g": {"wavelength": wl_g, "throughput": tp_g},
"Kepler": {"wavelength": wl_kep, "throughput": tp_kep},
"TESS": {"wavelength": wl_tess, "throughput": tp_tess},
})
print(f"Bands: {bps.band_names}")
print(f"Effective wavelengths: {bps.effective_wavelengths} Angstroms")
# 3. Precompute the contrast table
scm = SpectralContrastModel(
provider, bps,
T_phot=T_phot_true,
T_spot_grid=np.arange(2500, 5800, 25),
)
scm.summary()
3. Visualize the contrast table#
The precomputed table stores \(f_{\rm spot}(T_{\rm spot},\, {\rm band})\) – the bandpass-integrated flux ratio. The contrast factor is \(c = 1 - f_{\rm spot}\).
With BlackbodyProvider, the bandpass-integrated values are close to
the point-evaluation Planck ratio but not identical, because the finite
filter width introduces a weighted average across the bandpass.
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# Left: contrast factor from spectral model
for j, name in enumerate(bps.band_names):
axes[0].plot(scm.T_grid, 1 - scm.table[:, j],
label=name, color=list(colors.values())[j], lw=2)
axes[0].set_xlabel(r"$T_{\rm spot}$ [K]")
axes[0].set_ylabel(r"$c(\lambda) = 1 - f_{\rm spot}$")
axes[0].set_title("Bandpass-integrated contrast (SpectralContrastModel)")
axes[0].legend()
# Right: compare to point-evaluation Planck
T_grid = scm.T_grid
for j, name in enumerate(bps.band_names):
wl = band_info[name]["wavelength"]
c_spectral = 1 - scm.table[:, j]
c_planck = np.array([float(contrast_factor(jnp.array(wl), T, T_phot_true))
for T in T_grid])
residual = (c_spectral - c_planck) / np.where(c_planck > 0, c_planck, 1.0) * 100
axes[1].plot(T_grid, residual, label=name,
color=list(colors.values())[j], lw=1.5)
axes[1].axhline(0, color="gray", ls="--", lw=0.8)
axes[1].set_xlabel(r"$T_{\rm spot}$ [K]")
axes[1].set_ylabel("Residual [%]")
axes[1].set_title("Bandpass-integrated vs. point-evaluation Planck")
axes[1].legend()
fig.tight_layout()
plt.show()
4. Plug into MultiBandGPSolver#
Pass the SpectralContrastModel via the contrast_model argument.
The solver calls make_contrast_fn() internally to build a JAX-traceable
closure that replaces the default contrast_factor.
# Default solver (Planck point-evaluation)
gp_planck = MultiBandGPSolver(
data, hparam, T_phot=T_phot_true, T_spot_init=4000.0,
bounds=bounds, matrix_solver="cholesky_banded",
)
gp_planck.build_jax()
# Spectral contrast solver
gp_spectral = MultiBandGPSolver(
data, hparam, T_phot=T_phot_true, T_spot_init=4000.0,
bounds=bounds, contrast_model=scm,
matrix_solver="cholesky_banded",
)
gp_spectral.build_jax()
ll_planck = gp_planck.log_likelihood_at(gp_planck.theta0)
ll_spectral = gp_spectral.log_likelihood_at(gp_spectral.theta0)
print(f"log L (Planck point-eval): {ll_planck:.2f}")
print(f"log L (bandpass-integrated): {ll_spectral:.2f}")
print(f"Difference: {ll_spectral - ll_planck:.2f}")
When does this matter?
With BlackbodyProvider, the bandpass-integrated and point-evaluation
Planck results are very similar (within a few percent). The real advantage
comes from using KorgProvider with model-atmosphere spectra that include
molecular absorption features – especially for cool spots (\(T_{\rm spot}
\lesssim 4000\) K) where TiO bands dominate the optical.
5. Caching the contrast table#
For expensive spectrum providers like KorgProvider, you can cache the
precomputed table to disk with the cache_path argument. On subsequent
runs, the table is loaded from the .npz file if the configuration
(temperature grid, photosphere temperature, bandpass set) matches.
# Save to disk
scm_cached = SpectralContrastModel(
provider, bps, T_phot=T_phot_true,
T_spot_grid=np.arange(2500, 5800, 25),
cache_path="/tmp/spotgp_contrast_demo.npz",
)
# Load from disk (fast -- no recomputation)
scm_loaded = SpectralContrastModel(
provider, bps, T_phot=T_phot_true,
T_spot_grid=np.arange(2500, 5800, 25),
cache_path="/tmp/spotgp_contrast_demo.npz",
)
np.testing.assert_array_equal(scm_cached.table, scm_loaded.table)
print("Cache round-trip: OK")
6. Using Korg for model-atmosphere spectra#
KorgProvider uses Korg.jl for
1D LTE spectral synthesis with MARCS model atmospheres. This captures
molecular features (TiO, VO, H\(_2\)O) that dominate cool-star spectra
and significantly affect bandpass-integrated contrast below ~4000 K.
Installation
pip install juliacall juliapkg
python -c "import juliapkg; juliapkg.add('Korg', 'acafc109-a718-429c-b0e5-afd7f8c7ae46'); juliapkg.resolve()"
```python
from spotgp.spectral import KorgProvider
provider_korg = KorgProvider(
logg=4.44, # log g (cgs)
M_H=0.0, # [M/H]
wl_range=(3000, 11000),
linelist="vald_solar",
)
# First call triggers Julia JIT (~30-60 s), subsequent calls are fast
wl, flux = provider_korg.spectrum(5800.0)
scm_korg = SpectralContrastModel(
provider_korg, bps, T_phot=5800.0,
T_spot_grid=np.arange(2500, 5800, 50),
cache_path="contrast_korg.npz", # cache for reuse
)
scm_korg.summary()
7. Using pyphot for real filter curves#
If pyphot is installed, you can load real instrument filter profiles by name from its built-in library. This replaces the Gaussian approximations with the actual response curves.
Installation
pip install pyphot
```python
bps_real = BandpassSet({
"sdss_g": "SDSS_g",
"kepler": "KEPLER_Kp",
"tess": "TESS",
"2mass_j": "2MASS_J",
})
print(f"Effective wavelengths: {bps_real.effective_wavelengths}")
# Combine with any provider
scm_real = SpectralContrastModel(
provider_korg, bps_real, T_phot=5800.0,
T_spot_grid=np.arange(2500, 5800, 50),
cache_path="contrast_korg_real_filters.npz",
)
scm_real.summary()
You can also mix pyphot filters with custom curves:
bps_mixed = BandpassSet({
"tess": "TESS",
"my_filter": {"wavelength": wl_custom, "throughput": tp_custom},
})
Summary#
Concept |
Key class/function |
|---|---|
Planck spectrum (no dependencies) |
|
Model-atmosphere spectrum (Korg) |
|
Bandpass set (custom or pyphot) |
|
Precomputed contrast table |
|
Plug into solver |
|
Disk caching |
|
Next steps#
Compare Korg vs. blackbody contrast tables for M-dwarf temperatures
Run full MCMC with the spectral contrast model using BlackJAX sampling or Dynesty sampling
Extend to custom model grids (e.g., PHOENIX, BT-Settl) by writing a provider with a
.spectrum(Teff)method