"""Per-shutter NIRSpec dispersion model for MSA shutters.
All supported (disperser, filter) combinations use a *per-shutter*
lookup table derived from spacetelescope/msaviz's numerical
integration of the pipeline dispersion polynomials. The table lives
in `data/dispersion_cutoffs.npz` and is regenerated by
`scripts/precompute_dispersion_cutoffs.py` (re-run when the
underlying msaviz reference files change).
For each (disperser, filter) the table stores four (4, 171, 365)
float32 arrays — one slice per quadrant — under keys
{DISPERSER}_{FILTER}_blue_edge
{DISPERSER}_{FILTER}_gap_lo
{DISPERSER}_{FILTER}_gap_hi
{DISPERSER}_{FILTER}_red_edge
NaN means the shutter's spectrum doesn't reach the corresponding
detector. A linear V2-shift fallback is retained for the case where
the table is missing on disk.
"""
from pathlib import Path
import numpy as np
from .coords import MSA_V2_REF
# Per-grating fiducial wavelength ranges (microns), sourced from
# spacetelescope/msaviz's `ranges.json`, which is derived from the
# pipeline reference files. These are the wavelength bounds the
# detector ACTUALLY sees through each (disperser, filter) — wider
# than the JDox "useful range" in a few places (e.g. G140M/F100LP
# JDox quotes 0.97-1.84 but the pipeline reference is 0.97-1.89).
# We use the pipeline values so the per-shutter table lookups and
# this fallback agree.
GRATING_RANGES: dict[str, dict[str, tuple[float, float]]] = {
"PRISM": {"CLEAR": (0.60, 5.30)},
"G140M": {"F070LP": (0.70, 1.27), "F100LP": (0.97, 1.89)},
"G235M": {"F170LP": (1.66, 3.17)},
"G395M": {"F290LP": (2.87, 5.27)},
"G140H": {"F070LP": (0.70, 1.27), "F100LP": (0.97, 1.89)},
"G235H": {"F170LP": (1.66, 3.17)},
"G395H": {"F290LP": (2.87, 5.27)},
}
# Filter blue cutoffs (microns) — hard lower bounds.
FILTER_BLUE_CUTOFF: dict[str, float] = {
"CLEAR": 0.60,
"F070LP": 0.70,
"F100LP": 0.97,
"F170LP": 1.66,
"F290LP": 2.87,
}
# Projected MSA-to-detector dispersion span along V2 (arcsec). Placeholder per PLAN.md.
V2_DISP_EXTENT: float = 180.0
LAMBDA_PER_ARCSEC: float = (5.30 - 0.60) / V2_DISP_EXTENT # ~0.026 μm/arcsec
# Half-width of the on-detector spectrum, in V2 arcsec, per disperser.
#
# Source: eMPT (Bonaventura et al. 2023, A&A 672, A40),
# `reference_files/prism_sep.dat`, which tabulates per-shutter +/- column
# separations for the PRISM spectrum. At a central shutter the values are
# sep_p ≈ +176, sep_m ≈ -177 columns. Each MSA column is ~0.20" along V2,
# so the projected V2 half-extent of PRISM's spectrum is ~0.20 × 176 ≈ 35".
# Within ~10 % across the MSA.
#
# For the grating modes, eMPT applies no column cutoff at all — any pair
# of shutters at the same detector-y row collide (the spectrum spans the
# entire detector). For an H-grating the spectrum reaches beyond the MSA
# in V2, so cross-quadrant pairs collide as well. We approximate both
# behaviours with a single (large) V2 half-extent below.
#
# IMPORTANT — what the visualisation actually shows. We compute the
# spectral-conflict shutters as "same s row AND share a primary
# detector AND within `v2_overlap_distance` of the open shutter in V2."
# So `v2_overlap_distance(disperser, filt)` must return the *full*
# V2 extent of the dispersed spectrum on the detector — when two
# same-detector spectra of length L are offset by ΔV2, they share at
# least one detector pixel iff |ΔV2| < L.
#
# Values measured directly from `slit_frame → detector` traces in
# stenv: detector x span of the spectrum × ≈ 0.077 ″/V2-px. Per-combo
# because the filter sets the usable wavelength range, which the
# spectrum's x-span follows.
#
# Cross-detector pairs (e.g. Q4 G395M on NRS1 vs Q2 G395M on NRS2 at
# same d) are filtered out by the primary-detector lookup BEFORE this
# distance check fires — so two spectra at ΔV2 ≈ 87″ that live on
# different detector halves correctly report no overlap, even though
# 87″ < the 103″ full extent.
SPECTRUM_V2_EXTENT: dict[tuple[str, str], float] = {
("PRISM", "CLEAR"): 32.0,
("G140M", "F070LP"): 98.0,
("G140M", "F100LP"): 109.0,
("G235M", "F170LP"): 110.0,
("G395M", "F290LP"): 103.0,
("G140H", "F070LP"): 185.0,
("G140H", "F100LP"): 307.0,
("G235H", "F170LP"): 300.0,
("G395H", "F290LP"): 281.0,
}
# Fallback by disperser only, for unrecognised filter combos.
_FALLBACK_V2_EXTENT: dict[str, float] = {
"PRISM": 32.0,
"G140M": 109.0,
"G235M": 110.0,
"G395M": 103.0,
"G140H": 307.0,
"G235H": 300.0,
"G395H": 281.0,
}
# Backwards-compatibility aliases for the old "half-extent" names —
# downstream code (optimizer, tests) that imports them just gets the
# (now full) values.
SPECTRUM_V2_HALFEXTENT = SPECTRUM_V2_EXTENT
_FALLBACK_V2_HALFEXTENT = _FALLBACK_V2_EXTENT
# ── Grating diagonal-step relaxation ──────────────────────────────────────
# A "no-buffer" conflict between two slitlets whose rows are EXACTLY adjacent
# (1 row apart, no real row overlap) is painted purple (Mask Conflict) like
# APT MPT. But on the GRATING side a *column-offset* diagonal step (e.g.
# following an elongated galaxy) isn't a same-column collision and we demote
# purple → orange (Masked).
#
# Why the threshold is 1 column (any offset), not a large margin:
# two slitlets at adjacent ROWS stay separated by a FIXED ~1-row gap in the
# cross-dispersion direction no matter how far apart they are in COLUMNS —
# both spectra share the same trace tilt, so a column step only slides them
# along dispersion (detector x), never in cross-dispersion (detector y). And
# the spectra are far longer (M ~510 col, H ~1260 col) than the MSA is wide
# (365 col), so the x-overlap never closes within reach either. So the column-
# offset MAGNITUDE doesn't change the (marginal, 1-row) overlap at all — the
# only physically meaningful split is same-column (Δd=0, spectra fully stacked
# → real conflict, purple) vs different-column (Δd≥1, a deliberate diagonal
# step → orange). H and M therefore use the same value. PRISM is never relaxed
# (matches APT/MPT). Raise these if you want to require a wider deliberate step.
GRATING_ADJ_MIN_COLSEP_H = 1
GRATING_ADJ_MIN_COLSEP_M = 1
[docs]
def grating_adjacency_min_colsep(disperser: str) -> int | None:
"""Minimum column separation that demotes a no-buffer ADJACENCY conflict
(rows exactly 1 apart, no real row overlap) from purple (Mask Conflict)
to orange (Masked) for this disperser. ``1`` for any M/H grating — any
nonzero column offset is a deliberate diagonal step (see module note);
only exact same-column stacking (Δd=0) stays purple. Returns ``None`` when
adjacency must always stay purple — PRISM and anything not an M/H grating.
"""
d = (disperser or "").upper()
if d.startswith("G") and d.endswith("H"):
return GRATING_ADJ_MIN_COLSEP_H
if d.startswith("G") and d.endswith("M"):
return GRATING_ADJ_MIN_COLSEP_M
return None
[docs]
def v2_overlap_distance(disperser: str, filt: str) -> float:
"""Full V2 extent (arcsec) of the spectrum on the detector. Two
same-row shutters whose spectra land on the SAME detector and
whose V2 separation is less than this value share at least one
detector pixel — i.e. their spectra collide on the pipeline.
The same-detector test is the responsibility of the caller (use
:func:`primary_detector`). This function only returns the distance.
Looked up per (disperser, filter) combo because the filter sets
the usable wavelength range, which controls the spectrum's
detector x-span. Falls back to the disperser-only value if the
combo isn't in the measured table, and to a safe upper bound
(300″) for unknown dispersers."""
d = disperser.upper()
f = filt.upper()
if (d, f) in SPECTRUM_V2_EXTENT:
return SPECTRUM_V2_EXTENT[(d, f)]
if d in _FALLBACK_V2_EXTENT:
return _FALLBACK_V2_EXTENT[d]
return 300.0
[docs]
def tilt_slope_map(
disperser: str, filt: str,
) -> tuple[np.ndarray, np.ndarray, np.ndarray] | None:
"""Return the precomputed tilt-slope map for a (disperser, filter)
combo, or ``None`` if the table doesn't ship it.
The map describes how much the spectrum's cross-dispersion row
drifts as a function of V2 distance from the open shutter::
s_expected(ΔV2_arcsec) = s_open + slope * ΔV2_arcsec
where `s_open` is the MSA s-row (1-based) of the dispersing
shutter and `ΔV2_arcsec` is the candidate shutter's V2 offset
from it. The slope is in [MSA rows per arcsec V2] and is computed
from `slit_frame → detector` traces on a 10×10 grid per MSA
quadrant (see ``scripts/precompute_trace_tilt.py``).
Returns
-------
tuple or None
``(grid_rows, grid_cols, slope)`` — its parts:
- ``grid_rows`` — (Ngrid_r,) int16, vMPT row centres (1..171)
- ``grid_cols`` — (Ngrid_c,) int16, vMPT col centres (1..365)
- ``slope`` — (4, Ngrid_r, Ngrid_c) float32, tilt slope in rows
per arcsec V2; NaN where the shutter projects onto neither
detector for this combo.
``None`` when the precomputed table is missing or ships no tilt
arrays for this combo (older builds) — callers should fall back
to slope = 0 (flat-row overlap check).
"""
tbl = _load_dispersion_table()
if tbl is None:
return None
d = disperser.upper()
f = filt.upper()
key_r = f"{d}_{f}_tilt_grid_rows"
key_c = f"{d}_{f}_tilt_grid_cols"
key_s = f"{d}_{f}_tilt_slope"
if key_r not in tbl or key_c not in tbl or key_s not in tbl:
return None
return tbl[key_r], tbl[key_c], tbl[key_s]
[docs]
def tilt_slope_grid(disperser: str, filt: str) -> np.ndarray | None:
"""Return the full (4, 171, 365) tilt slope grid for (disperser,
filter), bilinearly interpolated from the 10×10 quadrant sample.
Returns ``None`` if the table or this combo isn't available; the
caller is expected to fall back to 0-slope (flat) in that case.
The cached interpolated grid lets the runtime overlap check
estimate each candidate shutter's spectrum-y-drift (used to widen
the y tolerance for cross-quadrant pairs), without recomputing
the per-shutter bilinear from scratch."""
cache_key = (disperser.upper(), filt.upper())
if cache_key in _TILT_SLOPE_GRID_CACHE:
return _TILT_SLOPE_GRID_CACHE[cache_key]
m = tilt_slope_map(disperser, filt)
if m is None:
_TILT_SLOPE_GRID_CACHE[cache_key] = None
return None
grid_rows, grid_cols, slope = m
s_targets = np.arange(1, 172, dtype=float)
d_targets = np.arange(1, 366, dtype=float)
s_mesh, d_mesh = np.meshgrid(s_targets, d_targets, indexing="ij")
out = np.full((4, 171, 365), np.nan, dtype=np.float32)
for qi in range(4):
out[qi] = _bilinear_finite_aware(
np.asarray(grid_rows, dtype=float),
np.asarray(grid_cols, dtype=float),
slope[qi],
s_mesh, d_mesh,
)
# NaN cells (all-NaN neighbourhoods) → 0 so callers can use the
# value unconditionally; same convention as `tilt_slope_for_shutter`.
out = np.where(np.isnan(out), 0.0, out)
_TILT_SLOPE_GRID_CACHE[cache_key] = out
return out
_TILT_SLOPE_GRID_CACHE: dict[tuple[str, str], np.ndarray | None] = {}
[docs]
def tilt_slope_for_shutter(
disperser: str, filt: str, q: int, s: int, d: int,
) -> float:
"""Bilinear-interpolated tilt slope (rows/arcsec V2) for a single
open shutter at (quadrant, s-row, d-col). Returns 0.0 if the table
is unavailable, the combo isn't covered, or all four grid corners
around (s, d) are NaN. Callers can use the returned value
unconditionally — a zero slope reduces the tilt-aware check back
to the flat row check.
"""
m = tilt_slope_map(disperser, filt)
if m is None:
return 0.0
grid_rows, grid_cols, slope = m
qi = int(q) - 1
if not (0 <= qi < slope.shape[0]):
return 0.0
s_clamped = float(max(grid_rows[0], min(grid_rows[-1], s)))
d_clamped = float(max(grid_cols[0], min(grid_cols[-1], d)))
# Bracket indices for the input (s, d) along each grid axis.
ri = int(np.searchsorted(grid_rows, s_clamped, side="right") - 1)
ri = max(0, min(ri, len(grid_rows) - 2))
ci = int(np.searchsorted(grid_cols, d_clamped, side="right") - 1)
ci = max(0, min(ci, len(grid_cols) - 2))
fr = (s_clamped - grid_rows[ri]) / (grid_rows[ri + 1] - grid_rows[ri])
fc = (d_clamped - grid_cols[ci]) / (grid_cols[ci + 1] - grid_cols[ci])
# Bilinear over the 4 corners — finite-aware: any NaN corner is
# given a zero weight; result is normalised by the surviving
# weights. If all four are NaN we return 0 (flat-row fallback).
corners = [
((1 - fr) * (1 - fc), slope[qi, ri, ci ]),
((1 - fr) * fc , slope[qi, ri, ci + 1]),
( fr * (1 - fc), slope[qi, ri + 1, ci ]),
( fr * fc , slope[qi, ri + 1, ci + 1]),
]
num = 0.0
den = 0.0
for w, v in corners:
if np.isfinite(v):
num += w * float(v)
den += w
if den <= 0:
return 0.0
return float(num / den)
[docs]
def disperser_range(disperser: str, filt: str) -> tuple[float, float] | None:
"""Nominal ``(lam_min, lam_max)`` in μm for a (disperser, filter)
combination.
Returns the wavelength range the chosen mode actually delivers to
the detector (the values in :data:`GRATING_RANGES`) — the same
numbers shown in the JDox "useful range" docs, except where
msaviz / the pipeline reference disagree (we follow the pipeline).
Returns ``None`` when the combination isn't supported (e.g.
``G140H`` + ``F290LP``); the optimizer treats that as "no
constraint can pass" and drops the affected sources.
"""
if disperser is None or filt is None:
return None
d = disperser.upper()
f = filt.upper()
bands = GRATING_RANGES.get(d)
if bands is None:
return None
return bands.get(f)
[docs]
def disperser_min_lambda(disperser: str, filt: str) -> float | None:
"""Bluest wavelength the (disperser, filter) combo can deliver.
``None`` when the combination isn't recognised."""
r = disperser_range(disperser, filt)
return None if r is None else r[0]
[docs]
def disperser_max_lambda(disperser: str, filt: str) -> float | None:
"""Reddest wavelength the (disperser, filter) combo can deliver.
``None`` when the combination isn't recognised."""
r = disperser_range(disperser, filt)
return None if r is None else r[1]
[docs]
def interval_covered(
lo: float, hi: float,
blue: float, gap_lo: float, gap_hi: float, red: float,
) -> bool:
"""Does the spectrum ``[blue, red]`` (with the detector gap
``[gap_lo, gap_hi]`` excluded) fully cover the requested
``[lo, hi]`` range?
Used by the per-target "required wavelength range" constraint.
Parameters
----------
lo, hi : float
Requested interval, in μm. Must satisfy ``lo <= hi``.
blue, red : float
Bluest / reddest λ the centre shutter delivers to the
detector. Both NaN when the source doesn't reach the
detector at all (cutoffs returns NaN for off-grid shutters).
gap_lo, gap_hi : float
NRS1/NRS2 detector-gap wavelength bounds for this shutter.
Both NaN ⇒ no gap (e.g. M-grating modes within their nominal
range).
Returns
-------
bool
True iff every wavelength in ``[lo, hi]`` lands somewhere on
the detector (i.e. inside ``[blue, gap_lo] ∪ [gap_hi, red]``,
where the gap is skipped iff it has finite bounds).
"""
if not (np.isfinite(lo) and np.isfinite(hi)):
return False
if not (np.isfinite(blue) and np.isfinite(red)):
return False
if lo > hi:
lo, hi = hi, lo
if lo < blue or hi > red:
return False
has_gap = np.isfinite(gap_lo) and np.isfinite(gap_hi)
if not has_gap:
return True
# Gap is INSIDE [blue, red]; reject if [lo, hi] dips into it.
# Strict: any wavelength in [lo, hi] that falls in (gap_lo, gap_hi)
# is missing.
return not (lo < gap_hi and hi > gap_lo)
# Fallback fractional gap parameters used when neither the per-
# shutter PRISM table NOR a per-disperser fiducial value is
# available. Width tightened from the previous 10 % so the gratings
# show a more believably narrow detector gap.
GAP_CENTER_REL: float = 0.50
GAP_WIDTH_REL: float = 0.04
# Per-disperser fiducial gap fallback (used only when the per-shutter
# PRISM table can't be located — kept as a safety net so vMPT still
# loads if `data/prism_cutoffs.npz` goes missing).
DETECTOR_GAP_FIDUCIAL: dict[str, tuple[float, float]] = {
"PRISM": (1.87, 3.93),
}
# Per-shutter dispersion bounds for every supported (disperser,
# filter) combo, generated by scripts/precompute_dispersion_cutoffs.py
# and shipped at vmpt/data/dispersion_cutoffs.npz. Loaded lazily on
# first lookup so import stays cheap.
_DISPERSION_TABLE_PATH = (
Path(__file__).resolve().parent / "data" / "dispersion_cutoffs.npz"
)
_DISPERSION_TABLE: dict | None = None
_DISPERSION_TABLE_LOAD_ATTEMPTED: bool = False
def _load_dispersion_table() -> dict | None:
"""Lazy-load the per-shutter dispersion bounds table. Returns
None if the npz isn't present (fresh checkout, fall back to the
linear-shift model)."""
global _DISPERSION_TABLE, _DISPERSION_TABLE_LOAD_ATTEMPTED
if _DISPERSION_TABLE_LOAD_ATTEMPTED:
return _DISPERSION_TABLE
_DISPERSION_TABLE_LOAD_ATTEMPTED = True
if not _DISPERSION_TABLE_PATH.exists():
return None
try:
# Eager-copy the relevant arrays so the npz handle can close.
with np.load(_DISPERSION_TABLE_PATH) as data:
_DISPERSION_TABLE = {key: data[key] for key in data.files}
except (OSError, KeyError, ValueError):
_DISPERSION_TABLE = None
return _DISPERSION_TABLE
def _table_lookup(
disperser: str, filt: str, q: int, s: int, d: int,
) -> tuple[float, float, float, float] | None:
"""Read the precomputed wavelength bounds for shutter (q, s, d)
under the given (disperser, filter). Returns None if the table or
the specific combo isn't loadable. Per-value NaN means the
spectrum doesn't reach that detector."""
tbl = _load_dispersion_table()
if tbl is None:
return None
qi, si, di = q - 1, s - 1, d - 1
if not (0 <= qi < 4 and 0 <= si < 171 and 0 <= di < 365):
return None
key_blue = f"{disperser}_{filt}_blue_edge"
key_glo = f"{disperser}_{filt}_gap_lo"
key_ghi = f"{disperser}_{filt}_gap_hi"
key_red = f"{disperser}_{filt}_red_edge"
if key_blue not in tbl:
return None
return (
float(tbl[key_blue][qi, si, di]),
float(tbl[key_glo][qi, si, di]),
float(tbl[key_ghi][qi, si, di]),
float(tbl[key_red][qi, si, di]),
)
# Cache of per-(disperser, filter) interpolated (4, 171, 365) arrays
# for the spec-overlap detector-pixel check. Computed lazily on first
# lookup so import stays cheap.
_SHUTTER_XY_CACHE: dict[tuple[str, str], dict | None] = {}
[docs]
def shutter_xy_grids(disperser: str, filt: str) -> dict | None:
"""Return interpolated per-shutter on-detector x/y ranges for the
given (disperser, filter), or ``None`` if the precomputed sample
grid isn't present.
Result is a dict of six (4, 171, 365) float32 arrays:
``x_lo_nrs1, x_hi_nrs1, y_nrs1``
``x_lo_nrs2, x_hi_nrs2, y_nrs2``
NaN at a shutter index means the spectrum doesn't reach that
detector. The 10×10 quadrant grid samples from
``precompute_trace_tilt.py`` are bilinearly interpolated to every
(s, d) so the spec-overlap check can do detector-pixel-accurate
interval intersection without re-running ``AssignWcsStep``.
Cached per (disperser, filter); call cost is one quick bilinear
pass the first time and then a dict lookup."""
key = (disperser.upper(), filt.upper())
if key in _SHUTTER_XY_CACHE:
return _SHUTTER_XY_CACHE[key]
tbl = _load_dispersion_table()
if tbl is None:
_SHUTTER_XY_CACHE[key] = None
return None
d_up, f_up = key
needed = [f"{d_up}_{f_up}_x_lo_nrs1", f"{d_up}_{f_up}_x_hi_nrs1",
f"{d_up}_{f_up}_y_nrs1",
f"{d_up}_{f_up}_x_lo_nrs2", f"{d_up}_{f_up}_x_hi_nrs2",
f"{d_up}_{f_up}_y_nrs2",
f"{d_up}_{f_up}_tilt_grid_rows",
f"{d_up}_{f_up}_tilt_grid_cols"]
if any(k not in tbl for k in needed):
_SHUTTER_XY_CACHE[key] = None
return None
grid_rows = np.asarray(tbl[f"{d_up}_{f_up}_tilt_grid_rows"], dtype=float)
grid_cols = np.asarray(tbl[f"{d_up}_{f_up}_tilt_grid_cols"], dtype=float)
# All-shutter (s, d) target grid (1-based to match the table).
s_all = np.arange(1, 172, dtype=float)
d_all = np.arange(1, 366, dtype=float)
s_mesh, d_mesh = np.meshgrid(s_all, d_all, indexing="ij") # (171, 365)
out: dict[str, np.ndarray] = {}
for arr_name in ("x_lo_nrs1", "x_hi_nrs1", "y_nrs1",
"x_lo_nrs2", "x_hi_nrs2", "y_nrs2"):
src = np.asarray(tbl[f"{d_up}_{f_up}_{arr_name}"], dtype=float)
interp_full = np.full((4, 171, 365), np.nan, dtype=np.float32)
for qi in range(4):
interp_full[qi] = _bilinear_finite_aware(
grid_rows, grid_cols, src[qi], s_mesh, d_mesh,
)
out[arr_name] = interp_full
_SHUTTER_XY_CACHE[key] = out
return out
def _bilinear_finite_aware(
grid_rows: np.ndarray, grid_cols: np.ndarray,
src: np.ndarray, s_targets: np.ndarray, d_targets: np.ndarray,
) -> np.ndarray:
"""Bilinear interpolation that drops NaN corners by re-normalising
the surviving weights — so an interior NaN in the 10×10 sample
grid (e.g. a shutter where AssignWcsStep returned no on-detector
samples) doesn't poison the whole neighbourhood."""
# Pick the bracketing sample cell from the CLAMPED position (keeps the
# cell index valid), but take the interpolation fraction from the
# UNclamped target. Inside the sample span this is ordinary bilinear
# interpolation; beyond the outermost samples it LINEARLY EXTRAPOLATES
# the smooth detector-x/y trend (fr/fc run outside [0, 1]).
#
# The 10×10 grid samples rows [16…155] and cols [34…331] only
# (`np.linspace(1, 171, 12)[1:-1]`-style, endpoints dropped), so MSA
# rows 1–15 / 156–171 and cols 1–33 / 332–365 are off-grid. Clamping
# them (the old behaviour) pinned every off-grid edge row to ONE
# detector-y, so the cross-quadrant overlap check saw a whole block of
# edge rows as sharing a detector row — e.g. all of Q2 rows 1–16 read
# as colliding with Q4 opens, flooding lower Q2 with spurious "Masked"
# shutters. Extrapolation restores the ~5 px/row gradient so off-grid
# rows get distinct, physical detector-y values.
s_sel = np.clip(s_targets, grid_rows[0], grid_rows[-1])
d_sel = np.clip(d_targets, grid_cols[0], grid_cols[-1])
ri = np.clip(
np.searchsorted(grid_rows, s_sel, side="right") - 1,
0, len(grid_rows) - 2,
)
ci = np.clip(
np.searchsorted(grid_cols, d_sel, side="right") - 1,
0, len(grid_cols) - 2,
)
fr = (s_targets - grid_rows[ri]) / (grid_rows[ri + 1] - grid_rows[ri])
fc = (d_targets - grid_cols[ci]) / (grid_cols[ci + 1] - grid_cols[ci])
fr = fr[..., None] if fr.ndim < s_targets.ndim else fr
fc = fc[..., None] if fc.ndim < d_targets.ndim else fc
fr = np.broadcast_to(fr, s_targets.shape)
fc = np.broadcast_to(fc, d_targets.shape)
# Four corner values + weights.
v00 = src[ri, ci ]
v01 = src[ri, ci + 1]
v10 = src[ri + 1, ci ]
v11 = src[ri + 1, ci + 1]
w00 = (1 - fr) * (1 - fc)
w01 = (1 - fr) * fc
w10 = fr * (1 - fc)
w11 = fr * fc
num = np.zeros_like(s_targets, dtype=float)
den = np.zeros_like(s_targets, dtype=float)
for v, w in ((v00, w00), (v01, w01), (v10, w10), (v11, w11)):
m = np.isfinite(v)
num = num + np.where(m, w * np.where(m, v, 0.0), 0.0)
den = den + np.where(m, w, 0.0)
with np.errstate(invalid="ignore", divide="ignore"):
out = num / den
out = np.where(den > 0, out, np.nan)
return out.astype(np.float32)
[docs]
def primary_detector(
disperser: str, filt: str, q: int, s: int, d: int,
) -> int:
"""Per-shutter "primary detector" lookup. Returns::
0 → spectrum's main body is on NRS1
1 → spectrum's main body is on NRS2
-1 → table not loaded / not in shipped combo / off-detector
The main body is the detector that carries the larger wavelength
range of the dispersed spectrum. For shutters whose spectrum is
fully on one detector this is trivially that detector; for
gap-spanning shutters it's whichever side covers more wavelength.
Derived from the per-detector wavelength bounds in
``data/dispersion_cutoffs.npz`` (populated by
``scripts/precompute_dispersion_cutoffs.py``).
The spec-overlap check uses this lookup to skip cross-detector
pairs — two shutters whose spectra don't share a detector main
body cannot collide on the pipeline, regardless of how close they
are in V2."""
tbl = _load_dispersion_table()
if tbl is None:
return -1
qi, si, di = q - 1, s - 1, d - 1
if not (0 <= qi < 4 and 0 <= si < 171 and 0 <= di < 365):
return -1
d_up = disperser.upper()
f_up = filt.upper()
key1_lo = f"{d_up}_{f_up}_nrs1_lo"
key1_hi = f"{d_up}_{f_up}_nrs1_hi"
key2_lo = f"{d_up}_{f_up}_nrs2_lo"
key2_hi = f"{d_up}_{f_up}_nrs2_hi"
if key1_lo not in tbl or key2_lo not in tbl:
return -1
n1_lo = float(tbl[key1_lo][qi, si, di])
n1_hi = float(tbl[key1_hi][qi, si, di])
n2_lo = float(tbl[key2_lo][qi, si, di])
n2_hi = float(tbl[key2_hi][qi, si, di])
n1_span = (n1_hi - n1_lo) if (np.isfinite(n1_lo) and np.isfinite(n1_hi)) else -1.0
n2_span = (n2_hi - n2_lo) if (np.isfinite(n2_lo) and np.isfinite(n2_hi)) else -1.0
if n1_span < 0 and n2_span < 0:
return -1
return 0 if n1_span >= n2_span else 1
[docs]
def primary_detector_grid(disperser: str, filt: str) -> np.ndarray | None:
"""Return the full (4, 171, 365) int8 array of primary detectors
for the (disperser, filter) combo. 0=NRS1, 1=NRS2, -1=neither/unknown.
Returns ``None`` if the table or per-detector arrays aren't
available; the live overlay falls back to static quadrant-pairing
in that case."""
tbl = _load_dispersion_table()
if tbl is None:
return None
d_up = disperser.upper()
f_up = filt.upper()
key1_lo = f"{d_up}_{f_up}_nrs1_lo"
key1_hi = f"{d_up}_{f_up}_nrs1_hi"
key2_lo = f"{d_up}_{f_up}_nrs2_lo"
key2_hi = f"{d_up}_{f_up}_nrs2_hi"
if any(k not in tbl for k in (key1_lo, key1_hi, key2_lo, key2_hi)):
return None
n1_lo = tbl[key1_lo]
n1_hi = tbl[key1_hi]
n2_lo = tbl[key2_lo]
n2_hi = tbl[key2_hi]
n1_span = np.where(
np.isfinite(n1_lo) & np.isfinite(n1_hi),
n1_hi - n1_lo, -1.0,
)
n2_span = np.where(
np.isfinite(n2_lo) & np.isfinite(n2_hi),
n2_hi - n2_lo, -1.0,
)
out = np.where(n1_span >= n2_span, 0, 1).astype(np.int8)
out[(n1_span < 0) & (n2_span < 0)] = -1
return out
[docs]
def cutoffs(
v2_arcsec: float,
v3_arcsec: float,
disperser: str,
filt: str,
*,
q: int | None = None,
s: int | None = None,
d: int | None = None,
) -> dict:
"""Wavelength endpoints of the dispersed spectrum on the detector
for a shutter at (V2, V3).
If `q, s, d` are supplied AND the precomputed dispersion table is
available for the requested (disperser, filter), returns the
per-shutter values from `data/dispersion_cutoffs.npz` (derived
from msaviz's integration of the pipeline dispersion models).
This is the accurate path; the gap location and spectrum edges
vary substantially across the MSA for every disperser, but
especially for PRISM (non-linear dispersion).
Without shutter indices, or if the table is missing, the function
falls back to a linear V2-shift model. The fallback exists so
vMPT still loads on a fresh checkout that hasn't run the
precompute, and so the existing test suite (which calls
`cutoffs(V2, V3, ...)` without indices) keeps working.
"""
disperser = disperser.upper()
filt = filt.upper()
if disperser not in GRATING_RANGES or filt not in GRATING_RANGES[disperser]:
raise ValueError(f"Unsupported (disperser, filter) = ({disperser}, {filt})")
lam_min, lam_max = GRATING_RANGES[disperser][filt]
blue_cut = FILTER_BLUE_CUTOFF.get(filt, 0.0)
def _maybe(lam, is_blue_edge: bool = False, lam_red_ref: float | None = None):
"""Apply filter blue-cut and convert NaN/None → None."""
if lam is None:
return None
if isinstance(lam, float) and (lam != lam): # NaN check w/o numpy
return None
if lam < blue_cut:
if is_blue_edge and lam_red_ref is not None and lam_red_ref > blue_cut:
return float(blue_cut)
return None
return float(lam)
# Per-shutter lookup — preferred when indices + table are available
# for the given (disperser, filter).
if q is not None and s is not None and d is not None:
hit = _table_lookup(disperser, filt, int(q), int(s), int(d))
if hit is not None:
blue, glo, ghi, red = hit
# If both bounds are NaN the shutter has no spectrum at all
# on either detector for this combo.
if blue != blue and red != red:
return {"lam_blue": None, "lam_gap_lo": None,
"lam_gap_hi": None, "lam_red": None}
return {
"lam_blue": _maybe(blue, is_blue_edge=True, lam_red_ref=red),
"lam_gap_lo": _maybe(glo),
"lam_gap_hi": _maybe(ghi),
"lam_red": _maybe(red),
}
# Fall through to fiducial / linear path if the lookup missed.
# Fallback / grating path: linear V2 shift, fiducial or fractional gap.
span = lam_max - lam_min
dlam_dv2 = span / V2_DISP_EXTENT
shift = (v2_arcsec - MSA_V2_REF) * dlam_dv2
if disperser == "PRISM":
# PRISM dispersion isn't linear in V2; without the table we
# hold endpoints at the disperser's full range (per msaviz
# the spread across shutters is ~0.01 μm anyway).
lam_blue_raw = lam_min
lam_red_raw = lam_max
else:
lam_blue_raw = lam_min + shift
lam_red_raw = lam_max + shift
lam_blue = max(lam_min, min(lam_max, lam_blue_raw))
lam_red = max(lam_min, min(lam_max, lam_red_raw))
fixed_gap = DETECTOR_GAP_FIDUCIAL.get(disperser)
if fixed_gap is not None:
lam_gap_lo_raw, lam_gap_hi_raw = fixed_gap
else:
gap_half = 0.5 * GAP_WIDTH_REL * span
lam_gap_center_raw = lam_min + GAP_CENTER_REL * span + shift
lam_gap_lo_raw = lam_gap_center_raw - gap_half
lam_gap_hi_raw = lam_gap_center_raw + gap_half
if lam_gap_hi_raw < lam_blue or lam_gap_lo_raw > lam_red:
lam_gap_lo = None
lam_gap_hi = None
else:
lam_gap_lo = max(lam_blue, min(lam_red, lam_gap_lo_raw))
lam_gap_hi = max(lam_blue, min(lam_red, lam_gap_hi_raw))
return {
"lam_blue": _maybe(lam_blue, is_blue_edge=True, lam_red_ref=lam_red),
"lam_gap_lo": _maybe(lam_gap_lo),
"lam_gap_hi": _maybe(lam_gap_hi),
"lam_red": _maybe(lam_red),
}