Instrument centered coordinates for Moon/Sun sidelobe cuts

[chervias]

This change was implemented with Claude Code (Opus 5.5). Edited since the approach changed a little bit.

It requires the companion so3g PR, which adds projection engines that multiply the detector quaternion on the left ("_DetLeft" engines, selected with the Projectionist constructor argument multiply_det_left).

Summary

Adds a way to build a pointing matrix in instrument-centered coordinates of a source (Moon, Sun, or any name accepted by planets.get_source_pos), and uses it to flag samples from per-wafer sidelobe masks delivered as observation metadata.

Instrument-centered coordinates give the position of the source relative to each detector, on axes fixed to the telescope structure (the focal plane xi/eta axes at zero roll). Unlike the previous object-centered frame (detector pointing relative to the source), a sidelobe that is fixed to the instrument stays fixed in this frame.

Quaternion chain

For detector d and sample t:

q[d,t] = A[d] * X[t]
X[t] = Rz(roll[t]) * q_bore[t]^-1 * q_src[t]            (sight.Q)
A[d] = Ry(pi/2) * Rz(roll0) * rotation_xieta(xi, eta, 0)^-1
       * Rz(-roll0)                                      (fp.quats)

The per-detector factor is on the left, which is why the so3g "_DetLeft" engines are needed.

• q_bore is the standard CelestialSightLine.az_el boresight (qpoint * Rz(pi + roll)); q_src = rotation_lonlat(ra, dec) of the source. Only q_bore^-1 * q_src enters, so the horizon->celestial rotation cancels and the full-precision celestial pointing is used.
• Rz(roll[t]) removes the focal-plane roll so the axes are fixed to the telescope structure; conjugating q_det by Rz(roll0) (circular mean roll) expresses the detector offsets in the same axes. Exact for constant roll; a warning is logged if roll varies by more than 0.1 deg from its mean.
• Detector quaternions are rebuilt with gamma = 0, since with q_det^-1 on the left gamma would rotate each detector's map about the detector.
• Ry(pi/2) moves the detector from the pole to (lon, lat) = (0, 0).

Map conventions: each detector is at (0, 0); for small offsets lat ~ eta_src - eta_det and lon ~ -(xi_src - xi_det), so in a standard CAR geometry right = +az and up = +el (commanded). For LAT observations past the zenith (el > 90) the axes follow the structure and are therefore flipped relative to the sky.

Changes

• coords/pmat.py:
  • P has a multiply_det_left attribute/constructor argument. When set, _get_proj passes multiply_det_left=True to Projectionist.for_geom / for_tiled / ProjectionistHealpix.for_healpix, so the so3g "_DetLeft" engines are used (CAR and HEALPix, comps='T' only). It is only passed when set, so P keeps working with older so3g in the default mode.
  • P.for_tod(..., instrument_centered=<source name>) builds the pointing above and sets multiply_det_left. Requires comps='T', geom given, rot=None and the boresight in az/el/roll (boresight_equ is not supported); raises ValueError otherwise.
  • New get_instrument_centered(sight, fp, roll, source, timestamps, site) returning the (sight, fp) pair for multiply_det_left projection.
• coords/planets.py: new get_source_quat(source, timestamps, site, dt=600): source RA/Dec from get_source_pos on a grid of at most dt seconds, RA unwrapped and linearly interpolated to each sample, returned as rotation_lonlat quaternions.
• coords/sidelobes.py: rewritten. Masks are no longer fixed files per object; they come from the observation metadata as aman.sidelobes (Enmap metadata loader): .map is an enmap [nfields, ny, nx] of 0/1 in the instrument-centered coordinates of each object, labelled by the "fields" LabelAxis (e.g. ['moon', 'sun']).
  • get_cuts(aman, sidelobes="sidelobes", object_list=None, rise_tol=1 deg) returns one RangesMatrix [ndet, nsamp] per object. Clear errors if the metadata is missing or the map has no WCS.
  • make_cuts: skips objects below rise_tol at the first, middle and last sample; projects the mask with P.for_tod(instrument_centered=...); uses unit responsivities so the mask value is independent of the detector T response (the old astype(int) could round 1*T < 1 down to 0); cuts samples where the projected mask is non-zero.
• site_pipeline/make_ml_map.py: call sidelobes.get_cuts(obs) with the new API. The --sun-mask/--moon-mask options are no longer used.
• tests/test_pmat.py: test_pmat_instrument_centered checks gamma independence, that the angular distance from the origin equals the detector-source distance on the sky, the orientation for the boresight detector against the Moon's az/el, that from_map/to_map run, and that comps='TQU' is rejected.

Validation on real data

• LAT observation at el = 120 deg (past the zenith), roll = 60 deg: the spread across detectors at fixed time is ~0.6 x 0.5 deg (one wafer), and the boresight orientation check (instrument-centered position of a detector at the boresight vs. skyfield's Moon az/el) agrees to ~1e-4 at all tested samples, consistent with refraction. This confirms the chain, the axis convention and the roll handling of X[t].
• LAT observation scanning across the Moon (roll = 0): using the samples where each detector sees the Moon in its main beam, the Moon centroids of 681 detectors coincide to ~5 arcmin, and a scan over the rotation applied to the detector offsets is minimized within ~2 deg of zero, confirming the xi/eta orientation of the detector offsets at zero roll. (The common offset of ~(+3, +17) arcmin is a pointing calibration offset of that observation, not related to this change.)
• Sidelobe cuts on the same observation: the Moon is within 0-27.5 deg of every detector throughout, inside the mask along the scan direction, so all samples are cut by design; the Sun is below the horizon and nothing is cut.

WARNING

This implementation can ONLY be used to go from a mask to flagged samples, "T only". It CANNOT be used to make maps in instrument-centered coordinates: polarization is not supported (gamma is set to 0, and so3g only provides T engines for this mode), and mapmaking in these coordinates has not been designed or validated. That will need a future implementation.

Known open point: the direction in which detector offsets are rotated by the roll (+roll vs -roll) has not been verified independently; the roll = 0 observation above cannot distinguish the two. The difference between the two hypotheses is ~2 r |sin(roll)|, where r is the detector's distance from the boresight (~2 deg for the wafer tested), so a Moon crossing with a few degrees of roll is enough to settle it.

🤖 Generated with Claude Code

[chervias]

this will fail tests without the so3g PR, so we need to merge that first