SAJAX vs. Fleck vs. Spotter: Light Curve Comparison#
This notebook benchmarks three stellar-activity codes across a range of configurations.
All three codes are tested on identical parameters. Residuals are plotted relative to SAJAX.
Tests:
Rotational modulation — equatorial spot, equator-on view
Rotational modulation — effect of stellar inclination
Transit light curve, no spots (LDC geometry test)
Spot-crossing anomaly (dark spot on transit chord)
Facula-crossing anomaly (bright region on transit chord)
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import warnings
warnings.filterwarnings('ignore')
import jax
import jax.numpy as jnp
# sajax (must be installed/on PYTHONPATH)
from sajax import compute_light_curve, compute_combined_light_curve
# ── Optional dependencies ────────────────────────────────────────────────────
FLECK_AVAILABLE = False
SPOTTER_AVAILABLE = False
BATMAN_AVAILABLE = False
try:
from fleck import Star as FleckStar
import astropy.units as u
FLECK_AVAILABLE = True
print("fleck : available")
except ImportError:
print("fleck : NOT available → pip install fleck")
try:
import batman
BATMAN_AVAILABLE = True
print("batman : available")
except ImportError:
print("batman : NOT available → pip install batman-package")
try:
from spotter import Star as SpotterStar
from spotter import core
from spotter.light_curves import light_curve
SPOTTER_AVAILABLE = True
print("spotter : available")
except ImportError:
print("spotter : NOT available → pip install spotter")
print(f"\nJAX version: {jax.__version__} | devices: {jax.devices()}")
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 11
7 import jax
8 import jax.numpy as jnp
9
10 # sajax (must be installed/on PYTHONPATH)
---> 11 from sajax import compute_light_curve, compute_combined_light_curve
12
13 # ── Optional dependencies ────────────────────────────────────────────────────
14 FLECK_AVAILABLE = False
ModuleNotFoundError: No module named 'sajax'
# ── Stellar ──────────────────────────────────────────────────────────────────
WAVELENGTH = np.array([550.0]) # [nm] single broadband bin
FLUX_QUIET = np.array([1.0]) # normalised quiet-star flux
U1, U2 = 0.1, 0.5 # quadratic LDC
P_ROT = 10.0 # [days] stellar rotation period
INC_STAR = 90.0 # [deg] stellar inclination (90 = equator-on)
VE = 2.0 # [km/s] equatorial velocity
GRID_SIZE = 250 # SAJAX stellar radius in pixels
# ── Active region ────────────────────────────────────────────────────────────
SPOT_LAT = 0.0 # [deg] equatorial spot
SPOT_LON = 0.0 # [deg] faces observer at t=0 / phase=0
THETA_DEG = 10.0 # [deg] angular radius
SPOT_CONTRAST = 0.7 # spot / quiet flux ratio
FACULA_CONTRAST = 1.15 # facula / quiet flux ratio
# fleck uses a *linear* radius fraction r_lin = sin(θ_angular)
R_LIN = float(np.sin(np.deg2rad(THETA_DEG)))
FLUX_SPOT = np.array([SPOT_CONTRAST])
FLUX_FAC = np.array([FACULA_CONTRAST])
# ── Planetary transit ────────────────────────────────────────────────────────
K = 0.1 # Rp / R*
P_ORB = 1.0 # [days]
T0 = 0.0 # [days] mid-transit epoch
A_OVER_R = 4.2 # a / R*
ORB_INC = 90.0 # [deg] edge-on orbit
# Approximate full transit duration T14 ≈ (P/π) arcsin((1+k)/a)
T14 = P_ORB / np.pi * np.arcsin((1.0 + K) / A_OVER_R)
print(f"Approximate T14 = {T14*24:.2f} h")
# ── Time / phase arrays ──────────────────────────────────────────────────────
N_ROT = 300
TIMES_ROT = np.linspace(0.0, P_ROT, N_ROT, endpoint=False) # [days]
PHASES_ROT = (TIMES_ROT / P_ROT * 360.0) % 360.0 # [deg], for SAJAX
N_TRANSIT = 400
TIMES_TRANSIT = np.linspace(-2.2 * T14, 2.2 * T14, N_TRANSIT) # [days]
# ── SAJAX parameter dicts ────────────────────────────────────────────────────
SAJAX_PARAMS = {
"inc_star" : INC_STAR,
"ldc_coeffs" : [U1, U2],
}
SAJAX_TRANSIT_PARAMS = {
"t0" : T0,
"period" : P_ORB,
"a_over_rstar": A_OVER_R,
"inclination" : np.deg2rad(ORB_INC),
"k" : K,
"ecc" : 0.0,
"omega_peri" : 0.0,
}
print(f"Spot angular radius : {THETA_DEG:.1f}° → fleck linear radius : {R_LIN:.4f}")
print(f"Transit times : {N_TRANSIT} pts over [{TIMES_TRANSIT[0]*24:.1f}, "
f"{TIMES_TRANSIT[-1]*24:.1f}] h")
Approximate T14 = 2.02 h
Spot angular radius : 10.0° → fleck linear radius : 0.1736
Transit times : 400 pts over [-4.5, 4.5] h
1. Rotational Modulation — Equatorial Spot, Equator-on View#
A single dark circular spot at latitude 0°, centred on the sub-observer meridian at t = 0. All three codes use the same limb-darkening coefficients and spot geometry.
Parameter |
Value |
|---|---|
Spot latitude |
0° (equator) |
Spot longitude |
0° (sub-observer at t=0) |
Angular radius |
10° |
Contrast (spot/quiet) |
0.7 |
Stellar inclination |
90° (equator-on) |
Rotation period |
10 days |
# ── SAJAX ────────────────────────────────────────────────────────────────────
result_sajax_rot = compute_light_curve(
wavelength = WAVELENGTH,
flux_quiet = FLUX_QUIET,
flux_active = FLUX_SPOT,
params = SAJAX_PARAMS,
ar_lat = [SPOT_LAT],
ar_long = [SPOT_LON],
ar_size = [THETA_DEG],
phases_rot = PHASES_ROT,
stellar_grid_size = GRID_SIZE,
ve = VE,
ldc_mode = "quadratic",
)
lc_sajax_rot = np.array(result_sajax_rot["lc"])
amp_sajax = (lc_sajax_rot.max() - lc_sajax_rot.min()) * 1e6
print(f"SAJAX amplitude: {amp_sajax:.1f} ppm")
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
SAJAX amplitude: 10243.5 ppm
# ── Fleck ────────────────────────────────────────────────────────────────────
if FLECK_AVAILABLE:
fleck_star = FleckStar(
spot_contrast = SPOT_CONTRAST,
n_phases = N_ROT,
u_ld = [U1, U2],
rotation_period = P_ROT,
)
# Arrays must be (n_spots, n_stars)
lons_f = np.array([[SPOT_LON]]) * u.deg
lats_f = np.array([[SPOT_LAT]]) * u.deg
radii_f = np.array([[R_LIN]]) # dimensionless linear radius
inc_f = np.array([INC_STAR]) * u.deg
lc_fleck_rot = np.squeeze(fleck_star.light_curve(lons_f, lats_f, radii_f, inc_f))
times_fleck_rot = fleck_star.phases / (2 * np.pi) * P_ROT # [days]
amp_fleck = (lc_fleck_rot.max() - lc_fleck_rot.min()) * 1e6
print(f"fleck amplitude: {amp_fleck:.1f} ppm")
else:
lc_fleck_rot = times_fleck_rot = None
print("fleck not available — skipping")
fleck amplitude: 10240.9 ppm
# ── Spotter ──────────────────────────────────────────────────────────────────
if SPOTTER_AVAILABLE:
nside = 200
sp_star = SpotterStar.from_sides(
nside,
inc = np.deg2rad(INC_STAR),
u = (U1, U2),
period = P_ROT,
)
_spot_mask = sp_star.spot(lat=np.deg2rad(SPOT_LAT), lon=np.deg2rad(SPOT_LON),
radius=np.deg2rad(THETA_DEG), sharpness=1000)
sp_star = sp_star.set(y=sp_star.y * (1.0 - (1.0 - SPOT_CONTRAST) * _spot_mask))
lc_spotter_rot = np.array(
light_curve(sp_star, jnp.asarray(TIMES_ROT, dtype=jnp.float32))
).ravel()
lc_spotter_rot = lc_spotter_rot / lc_spotter_rot.max()
amp_spotter = (lc_spotter_rot.max() - lc_spotter_rot.min()) * 1e6
print(f"spotter amplitude: {amp_spotter:.1f} ppm")
else:
lc_spotter_rot = None
print("spotter not available — skipping")
W0501 01:13:53.102563 16068316 cpp_gen_intrinsics.cc:74] Empty bitcode string provided for eigen. Optimizations relying on this IR will be disabled.
spotter amplitude: 10231.3 ppm
fig, axes = plt.subplots(2, 1, figsize=(12, 7), sharex=False,
gridspec_kw={"height_ratios": [3, 1.5]})
fig.suptitle(
f"Rotational Modulation — equatorial spot (θ={THETA_DEG}°, contrast={SPOT_CONTRAST})",
fontsize=13, fontweight="bold",
)
ax = axes[0]
ax.plot(TIMES_ROT, lc_sajax_rot, lw=2.5, color="steelblue", label="SAJAX", zorder=4)
if lc_fleck_rot is not None:
ax.plot(times_fleck_rot, lc_fleck_rot, lw=2, ls="--", color="darkorange", label="fleck")
if lc_spotter_rot is not None:
ax.plot(TIMES_ROT, lc_spotter_rot, lw=2, ls=":", color="forestgreen", label="spotter")
ax.set_ylabel("Normalised Flux")
ax.legend()
ax.grid(alpha=0.3)
ax2 = axes[1]
if lc_fleck_rot is not None:
ax2.plot(TIMES_ROT, (lc_sajax_rot - lc_fleck_rot) * 1e6,
lw=1.5, color="darkorange", label="SAJAX - fleck")
if lc_spotter_rot is not None:
ax2.plot(TIMES_ROT, (lc_sajax_rot - lc_spotter_rot) * 1e6,
lw=1.5, ls=":", color="forestgreen", label="SAJAX - spotter")
ax2.axhline(0, color="k", lw=1, ls="--")
ax2.set_xlabel("Time [days]")
ax2.set_ylabel("Residual [ppm]")
ax2.set_title("Residuals relative to SAJAX")
ax2.legend(fontsize=9)
ax2.grid(alpha=0.3)
plt.tight_layout()
plt.show()
print("Amplitudes:")
print(f" SAJAX : {amp_sajax:.1f} ppm")
if lc_fleck_rot is not None:
print(f" fleck : {amp_fleck:.1f} ppm ({abs(amp_fleck-amp_sajax)/amp_sajax*100:.2f}% diff)")
if lc_spotter_rot is not None:
print(f" spotter: {amp_spotter:.1f} ppm ({abs(amp_spotter-amp_sajax)/amp_sajax*100:.2f}% diff)")
Amplitudes:
SAJAX : 10243.5 ppm
fleck : 10240.9 ppm (0.03% diff)
spotter: 10231.3 ppm (0.12% diff)
2. Rotational Modulation — Effect of Stellar Inclination#
Same equatorial spot, but now we vary the stellar inclination from 30° (nearly pole-on) to 90° (equator-on). As the inclination decreases, the spot spends more time near the limb and the amplitude decreases. At exactly 0° (pole-on), an equatorial spot would rotate at the stellar equator and contribute no net modulation.
inclinations = [10.0, 20.0,30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0]
colors_inc = ["#e15759", "#f28e2b", "steelblue", "#4e79a7", "#76b7b2", "#59a14f", "#edc949", "#af7aa1", "#ff9da7"]
fig, axes = plt.subplot_mosaic([["A", "B"],["A", "C"],["D", "E"],["D", "F"]], figsize=(14, 10))
fig.suptitle("Rotational Modulation vs. Stellar Inclination", fontsize=13, fontweight="bold")
amps_by_inc = {"SAJAX": [], "fleck": [], "spotter": []}
for inc, col in zip(inclinations, colors_inc):
params_i = {"inc_star": inc, "ldc_coeffs": [U1, U2]}
# SAJAX
res = compute_light_curve(
wavelength=WAVELENGTH, flux_quiet=FLUX_QUIET, flux_active=FLUX_SPOT,
params=params_i, ar_lat=[SPOT_LAT], ar_long=[SPOT_LON], ar_size=[THETA_DEG],
phases_rot=PHASES_ROT, stellar_grid_size=GRID_SIZE, ve=VE, ldc_mode="quadratic",
)
lc_s = np.array(res["lc"])
amps_by_inc["SAJAX"].append((lc_s.max() - lc_s.min()) * 1e6)
# Fleck
if FLECK_AVAILABLE:
fs_i = FleckStar(spot_contrast=SPOT_CONTRAST, n_phases=N_ROT,
u_ld=[U1, U2], rotation_period=P_ROT)
lc_f = np.squeeze(fs_i.light_curve(
np.array([[SPOT_LON]]) * u.deg,
np.array([[SPOT_LAT]]) * u.deg,
np.array([[R_LIN]]),
np.array([inc]) * u.deg,
))
t_f = fs_i.phases / (2 * np.pi) * P_ROT
amps_by_inc["fleck"].append((lc_f.max() - lc_f.min()) * 1e6)
axes['A'].plot(t_f, (lc_s - lc_f)*1e6, lw=1.5, ls="--", color=col, label=f"fleck i={inc:.0f}°")
# Spotter
if SPOTTER_AVAILABLE:
sp_i = SpotterStar.from_sides(nside, inc=np.deg2rad(inc), u=(U1, U2), period=P_ROT)
_spot_mask = sp_i.spot(lat=np.deg2rad(SPOT_LAT), lon=np.deg2rad(SPOT_LON),
radius=np.deg2rad(THETA_DEG), sharpness=1000)
sp_i = sp_i.set(y=sp_i.y * (1.0 - (1.0 - SPOT_CONTRAST) * _spot_mask))
lc_sp = np.array(
light_curve(sp_i, jnp.asarray(TIMES_ROT, dtype=jnp.float32))
).ravel()
lc_sp = lc_sp / lc_sp.max()
amps_by_inc["spotter"].append((lc_sp.max() - lc_sp.min()) * 1e6)
axes['D'].plot(TIMES_ROT, (lc_s - lc_sp)*1e6, lw=1.5, ls=":", color=col, label=f"spotter i={inc:.0f}°")
axes['A'].set_xlabel("Time [days]")
axes['A'].set_ylabel("Residuals [ppm]")
axes['A'].set_title("Light Curves")
axes['A'].legend(fontsize=8, ncol=2)
axes['A'].grid(alpha=0.3)
axes['D'].set_xlabel("Time [days]")
axes['D'].set_ylabel("Residuals [ppm]")
axes['D'].set_title("Light Curves")
axes['D'].legend(fontsize=8, ncol=2)
axes['D'].grid(alpha=0.3)
# Amplitude vs inclination
axes['B'].plot(inclinations, amps_by_inc["SAJAX"], "o-", lw=2, ms=8, color="steelblue", label="SAJAX")
axes['E'].plot(inclinations, amps_by_inc["SAJAX"], "o-", lw=2, ms=8, color="steelblue", label="SAJAX")
if amps_by_inc["fleck"]:
axes['B'].plot(inclinations, amps_by_inc["fleck"], "s--", lw=2, ms=8, color="darkorange", label="fleck")
axes['C'].plot(inclinations, np.array(amps_by_inc["fleck"]) - np.array(amps_by_inc["SAJAX"]), "s--", lw=2, ms=8, color="darkorange", label="fleck - SAJAX")
if amps_by_inc["spotter"]:
axes['E'].plot(inclinations, amps_by_inc["spotter"], "^:", lw=2, ms=8, color="forestgreen", label="spotter")
axes['F'].plot(inclinations, np.array(amps_by_inc["spotter"]) - np.array(amps_by_inc["SAJAX"]), "^:", lw=2, ms=8, color="forestgreen", label="spotter - SAJAX")
axes['B'].set_ylabel("LC amplitude [ppm]")
axes['C'].set_ylabel("Residuals [ppm]")
axes['B'].sharex(axes['C'])
axes['B'].set_title("Amplitude vs. Stellar Inclination")
axes['C'].set_xlabel("Stellar inclination [deg]")
axes['B'].legend()
axes['C'].legend()
axes['B'].grid(alpha=0.3)
axes['C'].grid(alpha=0.3)
axes['E'].set_ylabel("LC amplitude [ppm]")
axes['F'].set_ylabel("Residuals [ppm]")
axes['E'].sharex(axes['F'])
axes['E'].set_title("Amplitude vs. Stellar Inclination")
axes['F'].set_xlabel("Stellar inclination [deg]")
axes['E'].legend()
axes['F'].legend()
axes['E'].grid(alpha=0.3)
axes['F'].grid(alpha=0.3)
plt.tight_layout()
plt.show()
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
3. Rotational Modulation — Effect of Spot size#
Same stellar inclination, but now we vary the spot size from small (5°) to very large (45°).
spot_sizes = [5.0, 10.0, 20.0, 30.0, 40.0, 45.0]
colors_inc = ["#4e79a7", "#76b7b2", "#59a14f", "#edc949", "#af7aa1", "#ff9da7"]
fig, axes = plt.subplot_mosaic([["A", "B"],["A", "C"],["D", "E"],["D", "F"]], figsize=(14, 10))
fig.suptitle("Rotational Modulation vs. Spot Size", fontsize=13, fontweight="bold")
amps_by_spot_size = {"SAJAX": [], "fleck": [], "spotter": []}
for size, col in zip(spot_sizes, colors_inc):
# SAJAX
res = compute_light_curve(
wavelength=WAVELENGTH, flux_quiet=FLUX_QUIET, flux_active=FLUX_SPOT,
params=SAJAX_PARAMS, ar_lat=[SPOT_LAT], ar_long=[SPOT_LON], ar_size=[size],
phases_rot=PHASES_ROT, stellar_grid_size=GRID_SIZE, ve=VE, ldc_mode="quadratic",
)
lc_s = np.array(res["lc"])
amps_by_spot_size["SAJAX"].append((lc_s.max() - lc_s.min()) * 1e6)
# Fleck
if FLECK_AVAILABLE:
fs_i = FleckStar(spot_contrast=SPOT_CONTRAST, n_phases=N_ROT,
u_ld=[U1, U2], rotation_period=P_ROT)
r_lin = float(np.sin(np.deg2rad(size)))
lc_f = np.squeeze(fs_i.light_curve(
np.array([[SPOT_LON]]) * u.deg,
np.array([[SPOT_LAT]]) * u.deg,
np.array([[r_lin]]),
np.array([INC_STAR]) * u.deg,
))
t_f = fs_i.phases / (2 * np.pi) * P_ROT
amps_by_spot_size["fleck"].append((lc_f.max() - lc_f.min()) * 1e6)
axes['A'].plot(t_f, (lc_s - lc_f)*1e6, lw=1.5, ls="--", color=col, label=f"spot size={size:.0f}°")
# Spotter
if SPOTTER_AVAILABLE:
spotsize_star = SpotterStar.from_sides(nside, inc=np.deg2rad(INC_STAR), u=(U1, U2), period=P_ROT)
_spot_mask = spotsize_star.spot(lat=np.deg2rad(SPOT_LAT), lon=np.deg2rad(SPOT_LON),
radius=np.deg2rad(size), sharpness=1000)
spotsize_star = spotsize_star.set(y=spotsize_star.y * (1.0 - (1.0 - SPOT_CONTRAST) * _spot_mask))
lc_sp = np.array(
light_curve(spotsize_star, jnp.asarray(TIMES_ROT, dtype=jnp.float32))
).ravel()
lc_sp = lc_sp / lc_sp.max()
amps_by_spot_size["spotter"].append((lc_sp.max() - lc_sp.min()) * 1e6)
axes['D'].plot(TIMES_ROT, (lc_s - lc_sp)*1e6, lw=1.5, ls=":", color=col, label=f"spot size={size:.0f}°")
axes['A'].set_xlabel("Time [days]")
axes['A'].set_ylabel("Residuals [ppm]")
axes['A'].set_title("Light Curves")
axes['A'].legend(fontsize=8, ncol=2)
axes['A'].grid(alpha=0.3)
axes['D'].set_xlabel("Time [days]")
axes['D'].set_ylabel("Residuals [ppm]")
axes['D'].set_title("Light Curves")
axes['D'].legend(fontsize=8, ncol=2)
axes['D'].grid(alpha=0.3)
# Amplitude vs spot size
axes['B'].plot(spot_sizes, amps_by_spot_size["SAJAX"], "o-", lw=2, ms=8, color="steelblue", label="SAJAX")
axes['E'].plot(spot_sizes, amps_by_spot_size["SAJAX"], "o-", lw=2, ms=8, color="steelblue", label="SAJAX")
if amps_by_spot_size["fleck"]:
axes['B'].plot(spot_sizes, amps_by_spot_size["fleck"], "s--", lw=2, ms=8, color="darkorange", label="fleck")
axes['C'].plot(spot_sizes, np.array(amps_by_spot_size["fleck"]) - np.array(amps_by_spot_size["SAJAX"]), "s--", lw=2, ms=8, color="darkorange", label="fleck - SAJAX")
if amps_by_spot_size["spotter"]:
axes['E'].plot(spot_sizes, amps_by_spot_size["spotter"], "^:", lw=2, ms=8, color="forestgreen", label="spotter")
axes['F'].plot(spot_sizes, np.array(amps_by_spot_size["spotter"]) - np.array(amps_by_spot_size["SAJAX"]), "^:", lw=2, ms=8, color="forestgreen", label="spotter - SAJAX")
axes['B'].set_ylabel("LC amplitude [ppm]")
axes['C'].set_ylabel("Residuals [ppm]")
axes['B'].sharex(axes['C'])
axes['B'].set_title("Amplitude vs. Spot size")
axes['C'].set_xlabel("Spot size [deg]")
axes['B'].legend()
axes['B'].grid(alpha=0.3)
axes['C'].grid(alpha=0.3)
axes['E'].set_ylabel("LC amplitude [ppm]")
axes['F'].set_ylabel("Residuals [ppm]")
axes['E'].sharex(axes['F'])
axes['E'].set_title("Amplitude vs. Spot size")
axes['F'].set_xlabel("Spot size [deg]")
axes['E'].legend()
axes['E'].grid(alpha=0.3)
axes['F'].grid(alpha=0.3)
plt.tight_layout()
plt.show()
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
3. Transit Light Curve — No Spots (LDC Geometry Test)#
A pure transit with no active regions. This checks that the limb-darkening profile and geometric transit depth agree across codes.
Parameter |
Value |
|---|---|
Radius ratio k |
0.1 (depth ≈ 1%) |
a / R★ |
4.2 |
Orbital inclination |
90° (edge-on) |
LDC (quadratic) |
u₁=0.5079, u₂=0.2239 |
For SAJAX, a tiny placeholder active region (size → 0) is needed to satisfy the API; it contributes no flux perturbation.
# ── SAJAX — transit, no spots ─────────────────────────────────────────────────
res_sajax_tr = compute_combined_light_curve(
wavelength = WAVELENGTH,
flux_quiet = FLUX_QUIET,
flux_active = FLUX_QUIET, # same as quiet → no perturbation
params = SAJAX_PARAMS,
ar_lat = [0.0],
ar_long = [0.0],
ar_size = [1e-4], # effectively zero
times = TIMES_TRANSIT,
P_rot = P_ROT,
transit_params = SAJAX_TRANSIT_PARAMS,
stellar_grid_size = GRID_SIZE,
ve = VE,
ldc_mode = "quadratic",
)
lc_sajax_tr = np.array(res_sajax_tr["lc"])
depth_sajax = 1.0 - lc_sajax_tr.min()
print(f"SAJAX transit depth: {depth_sajax*1e6:.0f} ppm ({depth_sajax*100:.4f}%)")
# ── Batman (reference analytical transit) ────────────────────────────────────
if BATMAN_AVAILABLE:
bp = batman.TransitParams()
bp.t0 = T0
bp.per = P_ORB
bp.rp = K
bp.a = A_OVER_R
bp.inc = ORB_INC
bp.ecc = 0.0
bp.w = 90.0
bp.limb_dark = "quadratic"
bp.u = [U1, U2]
bm = batman.TransitModel(bp, TIMES_TRANSIT)
lc_batman = bm.light_curve(bp)
depth_batman = 1.0 - lc_batman.min()
print(f"batman transit depth: {depth_batman*1e6:.0f} ppm ({depth_batman*100:.4f}%)")
else:
lc_batman = None
print("batman not available — skipping")
# ── Fleck — transit, no spots ─────────────────────────────────────────────────
# if FLECK_AVAILABLE and BATMAN_AVAILABLE:
# fleck_star_notspot = FleckStar(
# spot_contrast = 0.0,
# n_phases = N_ROT,
# u_ld = [U1, U2],
# rotation_period = P_ROT,
# )
# lc_fleck_tr = np.squeeze(
# fleck_star_notspot.light_curve(
# np.array([[1e-9]]) * u.deg,
# np.array([[1e-9]]) * u.deg,
# np.array([[1e-9]]), # 1e-9 not 0.0: zero radius creates a degenerate
# # shapely ellipse that breaks affinity.rotate
# INC_STAR * u.deg, # scalar Quantity — fleck requires this when planet= is set
# times = TIMES_TRANSIT,
# planet = bp,
# )
# )
# depth_fleck = 1.0 - lc_fleck_tr.min()
# print(f"fleck transit depth: {depth_fleck*1e6:.0f} ppm ({depth_fleck*100:.4f}%)")
# else:
lc_fleck_tr = None
# spotter has no pixel-level transit integration API
lc_spotter_tr = None
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
SAJAX transit depth: 11312 ppm (1.1312%)
batman transit depth: 11318 ppm (1.1318%)
fig, axes = plt.subplots(2, 1, figsize=(12, 7), sharex=True,
gridspec_kw={"height_ratios": [3, 1.5]})
fig.suptitle(f"Transit (no spots) — k={K}, a/R*={A_OVER_R}, u=({U1},{U2})",
fontsize=13, fontweight="bold")
ax = axes[0]
ax.plot(TIMES_TRANSIT * 24, lc_sajax_tr, lw=2.5, color="steelblue", label="SAJAX", zorder=5)
if lc_batman is not None:
ax.plot(TIMES_TRANSIT * 24, lc_batman, lw=1.5, ls="--", color="k", label="batman (ref.)")
if lc_fleck_tr is not None:
ax.plot(TIMES_TRANSIT * 24, lc_fleck_tr, lw=2, ls="-.", color="darkorange", label="fleck")
if lc_spotter_tr is not None:
ax.plot(TIMES_TRANSIT * 24, lc_spotter_tr, lw=2, ls=":", color="forestgreen", label="spotter")
ax.set_ylabel("Normalised Flux")
ax.legend()
ax.grid(alpha=0.3)
ax2 = axes[1]
if lc_batman is not None:
ax2.plot(TIMES_TRANSIT * 24, (lc_sajax_tr - lc_batman) * 1e6,
lw=1.5, color="steelblue", label="SAJAX − batman")
if lc_fleck_tr is not None:
ax2.plot(TIMES_TRANSIT * 24, (lc_fleck_tr - lc_batman) * 1e6,
lw=1.5, ls="-.", color="darkorange", label="fleck − batman")
if lc_spotter_tr is not None and lc_batman is not None:
ax2.plot(TIMES_TRANSIT * 24, (lc_spotter_tr - lc_batman) * 1e6,
lw=1.5, ls=":", color="forestgreen", label="spotter − batman")
ax2.axhline(0, color="k", lw=1, ls="--")
ax2.set_xlabel("Time from mid-transit [hours]")
ax2.set_ylabel("Residual [ppm]")
ax2.set_title("Residuals relative to batman (analytical)")
ax2.legend(fontsize=9)
ax2.grid(alpha=0.3)
plt.tight_layout()
plt.show()
4. Spot-Crossing Anomaly#
A dark spot is placed on the transit chord (lat=0°, lon=0°, facing the observer at t₀=0). As the planet crosses the spot, it occults darker-than-average stellar surface, so the observed flux rises above the unperturbed transit — the classic spot-crossing bump.
In SAJAX, the planet mask is applied at the pixel level inside the flux integrator, so the crossing is computed exactly without any post-multiplication approximation. Fleck uses a similar pixel-level approach via batman.
Parameter |
Value |
|---|---|
Spot lat / lon |
0° / 0° |
Spot angular radius |
10° |
Spot contrast |
0.7 |
Spot position relative to transit |
on transit chord (direct crossing) |
# ── SAJAX — spot crossing ────────────────────────────────────────────────────
res_sajax_sc = compute_combined_light_curve(
wavelength = WAVELENGTH,
flux_quiet = FLUX_QUIET,
flux_active = FLUX_SPOT,
params = SAJAX_PARAMS,
ar_lat = [SPOT_LAT],
ar_long = [SPOT_LON],
ar_size = [THETA_DEG],
times = TIMES_TRANSIT,
P_rot = P_ROT,
transit_params = SAJAX_TRANSIT_PARAMS,
stellar_grid_size = GRID_SIZE,
ve = VE,
ldc_mode = "quadratic",
)
lc_sajax_sc = np.array(res_sajax_sc["lc"])
bump_sajax = (lc_sajax_sc - lc_sajax_tr).max() * 1e6
print(f"SAJAX spot-crossing bump height: {bump_sajax:.1f} ppm")
# ── Fleck — spot crossing ─────────────────────────────────────────────────────
# if FLECK_AVAILABLE and BATMAN_AVAILABLE:
# fleck_star_sc = FleckStar(
# spot_contrast = SPOT_CONTRAST,
# n_phases = N_ROT,
# u_ld = [U1, U2],
# rotation_period = P_ROT,
# )
# lc_fleck_sc = np.squeeze(
# fleck_star_sc.light_curve(
# np.array([[SPOT_LON]]) * u.deg,
# np.array([[SPOT_LAT]]) * u.deg,
# np.array([[R_LIN]]),
# INC_STAR * u.deg, # scalar — required when planet= is set
# times = TIMES_TRANSIT,
# planet = bp,
# )
# )
# bump_fleck = (lc_fleck_sc - lc_fleck_tr).max() * 1e6
# print(f"fleck spot-crossing bump height: {bump_fleck:.1f} ppm")
# else:
lc_fleck_sc = None
# spotter has no pixel-level transit integration API
lc_spotter_sc = None
build_model: scalar LDCs provided for 'quadratic' law ([0.1000, 0.5000]) — broadcasting across all 1 wavelength bins.
build_model: active region overlap mode: 'hottest_wins' (overlaps take flux from hottest AR)
SAJAX spot-crossing bump height: -6839.8 ppm
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
fig.suptitle(f"Spot-Crossing Anomaly (spot contrast={SPOT_CONTRAST}, θ={THETA_DEG}°)",
fontsize=13, fontweight="bold")
t_h = TIMES_TRANSIT * 24 # hours
# ── Panel A: full transit with crossing ──────────────────────────────────────
ax = axes[0]
ax.plot(t_h, lc_sajax_sc, lw=2.5, color="steelblue", label="SAJAX", zorder=5)
if lc_fleck_sc is not None:
ax.plot(t_h, lc_fleck_sc, lw=2, ls="--", color="darkorange", label="fleck")
if lc_spotter_sc is not None:
ax.plot(t_h, lc_spotter_sc, lw=2, ls=":", color="forestgreen", label="spotter")
if lc_batman is not None:
ax.plot(t_h, lc_batman, lw=1, ls=":", color="grey", alpha=0.6, label="no-spot transit")
ax.set_xlabel("Time [hours]")
ax.set_ylabel("Normalised Flux")
ax.set_title("Transit + Spot Crossing")
ax.legend(fontsize=9)
ax.grid(alpha=0.3)
# ── Panel B: crossing anomaly = LC(with spot) − LC(no spot) ──────────────────
ax = axes[1]
anom_sajax = (lc_sajax_sc - lc_sajax_tr) * 1e6
ax.plot(t_h, anom_sajax, lw=2.5, color="steelblue", label=f"SAJAX peak={anom_sajax.max():.0f} ppm")
if lc_fleck_sc is not None:
anom_fleck = (lc_fleck_sc - lc_fleck_tr) * 1e6
ax.plot(t_h, anom_fleck, lw=2, ls="--", color="darkorange",
label=f"fleck peak={anom_fleck.max():.0f} ppm")
if lc_spotter_sc is not None and lc_spotter_tr is not None:
anom_spotter = (lc_spotter_sc - lc_spotter_tr) * 1e6
ax.plot(t_h, anom_spotter, lw=2, ls=":", color="forestgreen",
label=f"spotter peak={anom_spotter.max():.0f} ppm")
ax.axhline(0, color="k", lw=1, ls="--")
ax.set_xlabel("Time [hours]")
ax.set_ylabel("Anomaly [ppm]")
ax.set_title("Spot-Crossing Anomaly")
ax.legend(fontsize=9)
ax.grid(alpha=0.3)
# ── Panel C: code-to-code residuals ──────────────────────────────────────────
ax = axes[2]
if lc_fleck_sc is not None:
ax.plot(t_h, (lc_sajax_sc - lc_fleck_sc) * 1e6, lw=2, color="darkorange", label="SAJAX − fleck")
if lc_spotter_sc is not None:
ax.plot(t_h, (lc_sajax_sc - lc_spotter_sc) * 1e6, lw=2, ls=":", color="forestgreen",
label="SAJAX − spotter")
ax.axhline(0, color="k", lw=1, ls="--")
ax.set_xlabel("Time [hours]")
ax.set_ylabel("Residual [ppm]")
ax.set_title("Code-to-Code Residuals")
ax.legend(fontsize=9)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
5. Facula-Crossing Anomaly#
Same geometry as the spot-crossing test, but the active region is now a facula —
a bright region with contrast = 1.15 (15% brighter than the quiet star).
When the planet occults the facula, the observed flux dips below the unperturbed
transit level — the signal is the mirror image of the spot-crossing bump.
This tests whether each code correctly handles both the sign and magnitude of active-region contrast.
# ── SAJAX — facula crossing ───────────────────────────────────────────────────
res_sajax_fc = compute_combined_light_curve(
wavelength = WAVELENGTH,
flux_quiet = FLUX_QUIET,
flux_active = FLUX_FAC,
params = SAJAX_PARAMS,
ar_lat = [SPOT_LAT],
ar_long = [SPOT_LON],
ar_size = [THETA_DEG],
times = TIMES_TRANSIT,
P_rot = P_ROT,
transit_params = SAJAX_TRANSIT_PARAMS,
stellar_grid_size = GRID_SIZE,
ve = VE,
ldc_mode = "quadratic",
)
lc_sajax_fc = np.array(res_sajax_fc["lc"])
dip_sajax = (lc_sajax_fc - lc_sajax_tr).min() * 1e6
print(f"SAJAX facula-crossing dip: {dip_sajax:.1f} ppm")
# ── Fleck — facula crossing ───────────────────────────────────────────────────
if FLECK_AVAILABLE and BATMAN_AVAILABLE:
fleck_star_fc = FleckStar(
spot_contrast = FACULA_CONTRAST,
n_phases = N_ROT,
u_ld = [U1, U2],
rotation_period = P_ROT,
)
lc_fleck_fc = np.squeeze(
fleck_star_fc.light_curve(
np.array([[SPOT_LON]]) * u.deg,
np.array([[SPOT_LAT]]) * u.deg,
np.array([[R_LIN]]),
INC_STAR * u.deg, # scalar — required when planet= is set
times = TIMES_TRANSIT,
planet = bp,
)
)
dip_fleck = (lc_fleck_fc - lc_fleck_tr).min() * 1e6
print(f"fleck facula-crossing dip: {dip_fleck:.1f} ppm")
else:
lc_fleck_fc = None
# spotter has no pixel-level transit integration API
lc_spotter_fc = None
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
fig.suptitle(f"Facula-Crossing Anomaly (facula contrast={FACULA_CONTRAST}, θ={THETA_DEG}°)",
fontsize=13, fontweight="bold")
t_h = TIMES_TRANSIT * 24
# ── Panel A: full transit ─────────────────────────────────────────────────────
ax = axes[0]
ax.plot(t_h, lc_sajax_fc, lw=2.5, color="steelblue", label="SAJAX", zorder=5)
if lc_fleck_fc is not None:
ax.plot(t_h, lc_fleck_fc, lw=2, ls="--", color="darkorange", label="fleck")
if lc_spotter_fc is not None:
ax.plot(t_h, lc_spotter_fc, lw=2, ls=":", color="forestgreen", label="spotter")
if lc_batman is not None:
ax.plot(t_h, lc_batman, lw=1, ls=":", color="grey", alpha=0.6, label="no-fac transit")
ax.set_xlabel("Time [hours]")
ax.set_ylabel("Normalised Flux")
ax.set_title("Transit + Facula Crossing")
ax.legend(fontsize=9)
ax.grid(alpha=0.3)
# ── Panel B: anomaly = LC(with facula) − LC(no active region) ────────────────
ax = axes[1]
anom_sajax_fc = (lc_sajax_fc - lc_sajax_tr) * 1e6
ax.plot(t_h, anom_sajax_fc, lw=2.5, color="steelblue",
label=f"SAJAX dip={anom_sajax_fc.min():.0f} ppm")
if lc_fleck_fc is not None:
anom_fleck_fc = (lc_fleck_fc - lc_fleck_tr) * 1e6
ax.plot(t_h, anom_fleck_fc, lw=2, ls="--", color="darkorange",
label=f"fleck dip={anom_fleck_fc.min():.0f} ppm")
if lc_spotter_fc is not None and lc_spotter_tr is not None:
anom_spotter_fc = (lc_spotter_fc - lc_spotter_tr) * 1e6
ax.plot(t_h, anom_spotter_fc, lw=2, ls=":", color="forestgreen",
label=f"spotter dip={anom_spotter_fc.min():.0f} ppm")
ax.axhline(0, color="k", lw=1, ls="--")
ax.set_xlabel("Time [hours]")
ax.set_ylabel("Anomaly [ppm]")
ax.set_title("Facula-Crossing Anomaly")
ax.legend(fontsize=9)
ax.grid(alpha=0.3)
# ── Panel C: code-to-code residuals ──────────────────────────────────────────
ax = axes[2]
if lc_fleck_fc is not None:
ax.plot(t_h, (lc_sajax_fc - lc_fleck_fc) * 1e6, lw=2, color="darkorange",
label="SAJAX − fleck")
if lc_spotter_fc is not None:
ax.plot(t_h, (lc_sajax_fc - lc_spotter_fc) * 1e6, lw=2, ls=":", color="forestgreen",
label="SAJAX − spotter")
ax.axhline(0, color="k", lw=1, ls="--")
ax.set_xlabel("Time [hours]")
ax.set_ylabel("Residual [ppm]")
ax.set_title("Code-to-Code Residuals")
ax.legend(fontsize=9)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
6. Off-Chord Spot (Impact Parameter Mismatch)#
The planet now has an impact parameter b = 0.4, but the spot remains at the equator (lat=0°). The transit chord misses the spot centre, so any anomaly is strongly suppressed. This tests limb geometry and partial occultation handling.
B_IMPACT = 0.4
ORB_INC_B = float(np.rad2deg(np.arccos(B_IMPACT / A_OVER_R))) # inclination for b=0.4
transit_params_b = {
"t0" : T0,
"period" : P_ORB,
"a_over_rstar": A_OVER_R,
"inclination" : np.deg2rad(ORB_INC_B),
"k" : K,
"ecc" : 0.0,
"omega_peri" : 0.0,
}
print(f"Impact parameter b = {B_IMPACT} → inclination = {ORB_INC_B:.3f}°")
# ── SAJAX — b=0.4, with spot ──────────────────────────────────────────────────
res_sajax_b = compute_combined_light_curve(
wavelength = WAVELENGTH,
flux_quiet = FLUX_QUIET,
flux_active = FLUX_SPOT,
params = SAJAX_PARAMS,
ar_lat = [SPOT_LAT],
ar_long = [SPOT_LON],
ar_size = [THETA_DEG],
times = TIMES_TRANSIT,
P_rot = P_ROT,
transit_params = transit_params_b,
stellar_grid_size = GRID_SIZE,
ve = VE,
ldc_mode = "quadratic",
)
lc_sajax_b = np.array(res_sajax_b["lc"])
# ── SAJAX — b=0.4, no spot ────────────────────────────────────────────────────
res_sajax_b0 = compute_combined_light_curve(
wavelength=WAVELENGTH, flux_quiet=FLUX_QUIET, flux_active=FLUX_QUIET,
params=SAJAX_PARAMS, ar_lat=[0.0], ar_long=[0.0], ar_size=[1e-4],
times=TIMES_TRANSIT, P_rot=P_ROT, transit_params=transit_params_b,
stellar_grid_size=GRID_SIZE, ve=VE, ldc_mode="quadratic",
)
lc_sajax_b0 = np.array(res_sajax_b0["lc"])
# ── Fleck — b=0.4 ─────────────────────────────────────────────────────────────
if FLECK_AVAILABLE and BATMAN_AVAILABLE:
bp_b = batman.TransitParams()
bp_b.t0 = T0; bp_b.per = P_ORB; bp_b.rp = K; bp_b.a = A_OVER_R
bp_b.inc = ORB_INC_B; bp_b.ecc = 0.0; bp_b.w = 90.0
bp_b.limb_dark = "quadratic"; bp_b.u = [U1, U2]
lc_fleck_b = np.squeeze(
fleck_star_sc.light_curve(
np.array([[SPOT_LON]]) * u.deg,
np.array([[SPOT_LAT]]) * u.deg,
np.array([[R_LIN]]),
INC_STAR * u.deg, # scalar — required when planet= is set
times = TIMES_TRANSIT,
planet = bp_b,
)
)
else:
lc_fleck_b = None
# ── Plot ──────────────────────────────────────────────────────────────────────
fig, axes = plt.subplots(1, 2, figsize=(13, 5))
fig.suptitle(f"Off-Chord Spot (b={B_IMPACT}, spot at equator)",
fontsize=13, fontweight="bold")
t_h = TIMES_TRANSIT * 24
ax = axes[0]
ax.plot(t_h, lc_sajax_b, lw=2.5, color="steelblue", label="SAJAX (spot)")
ax.plot(t_h, lc_sajax_b0, lw=1.5, ls="--", color="steelblue", alpha=0.5, label="SAJAX (no spot)")
if lc_fleck_b is not None:
ax.plot(t_h, lc_fleck_b, lw=2, ls="-.", color="darkorange", label="fleck (spot)")
ax.set_xlabel("Time [hours]"); ax.set_ylabel("Normalised Flux")
ax.set_title("Transit Light Curves")
ax.legend(fontsize=9); ax.grid(alpha=0.3)
ax = axes[1]
anom_b = (lc_sajax_b - lc_sajax_b0) * 1e6
ax.plot(t_h, anom_b, lw=2.5, color="steelblue",
label=f"SAJAX anomaly peak={anom_b.max():.0f} ppm")
ax.plot(t_h, anom_sajax, lw=1.5, ls=":", color="steelblue", alpha=0.5,
label="SAJAX anomaly b=0 (ref)")
if lc_fleck_b is not None:
bm_b = batman.TransitModel(bp_b, TIMES_TRANSIT)
lc_batman_b = bm_b.light_curve(bp_b)
anom_fleck_b = (lc_fleck_b - lc_batman_b) * 1e6
ax.plot(t_h, anom_fleck_b, lw=2, ls="-.", color="darkorange",
label=f"fleck anomaly peak={anom_fleck_b.max():.0f} ppm")
ax.axhline(0, color="k", lw=1, ls="--")
ax.set_xlabel("Time [hours]"); ax.set_ylabel("Anomaly [ppm]")
ax.set_title("Spot-Crossing Anomaly (suppressed for b=0.4)")
ax.legend(fontsize=9); ax.grid(alpha=0.3)
plt.tight_layout()
plt.show()
Summary#
The table below collects the key metrics from each test.
rows = []
# Test 1: rotational amplitude
row = {"Test": "1. Rotational amplitude [ppm]",
"SAJAX": f"{amp_sajax:.1f}"}
if lc_fleck_rot is not None:
row["fleck"] = f"{amp_fleck:.1f} ({(amp_fleck-amp_sajax)/amp_sajax*100:+.2f}%)"
if lc_spotter_rot is not None:
row["spotter"] = f"{amp_spotter:.1f} ({(amp_spotter-amp_sajax)/amp_sajax*100:+.2f}%)"
rows.append(row)
# Test 3: transit depth
row = {"Test": "3. Transit depth (no spots) [ppm]",
"SAJAX": f"{depth_sajax*1e6:.0f}"}
if lc_batman is not None:
row["batman"] = f"{depth_batman*1e6:.0f} (ref)"
if lc_fleck_tr is not None:
row["fleck"] = f"{depth_fleck*1e6:.0f} ({(depth_fleck-depth_sajax)/depth_sajax*100:+.3f}%)"
rows.append(row)
# Test 4: spot-crossing bump
row = {"Test": "4. Spot-crossing bump [ppm]",
"SAJAX": f"{bump_sajax:.1f}"}
if lc_fleck_sc is not None:
row["fleck"] = f"{bump_fleck:.1f} ({(bump_fleck-bump_sajax)/bump_sajax*100:+.2f}%)"
rows.append(row)
# Test 5: facula-crossing dip
row = {"Test": "5. Facula-crossing dip [ppm]",
"SAJAX": f"{dip_sajax:.1f}"}
if lc_fleck_fc is not None:
row["fleck"] = f"{dip_fleck:.1f} ({(dip_fleck-dip_sajax)/dip_sajax*100:+.2f}%)"
rows.append(row)
# Print as table
col_w = 42
header_keys = list({k for r in rows for k in r if k != "Test"})
header_keys.sort()
print(f"{'Test':<38}" + "".join(f"{h:<{col_w}}" for h in header_keys))
print("-" * (38 + col_w * len(header_keys)))
for r in rows:
print(f"{r['Test']:<38}" + "".join(f"{r.get(h,'—'):<{col_w}}" for h in header_keys))
print("""
Notes
-----
* Rotational modulation differences arise from grid resolution:
SAJAX uses a pixel grid (radius=100px), fleck uses an analytic cap model,
spotter uses HEALPix nside=32 (~1.8° pixel scale).
* Transit depth differences reflect the limb-darkening integration method.
Batman (analytic) is the ground truth; pixel-grid codes introduce ≲1 ppm error.
* Spot/facula-crossing anomaly differences include both the rotational-modulation
offset at transit time and the occultation geometry.
""")