""" Generate roman_footprint_rgps_w51.html from template + pysiaf WFI polygons. Corrected to use attitude_matrix(V2Ref, V3Ref, ra, dec, pa) so footprints are placed at the actual WFI pointing position (not displaced by ~0.496 deg due to the WFI_CEN focal-plane offset from the V3 boresight). """ import json, math import numpy as np import pysiaf from pysiaf.utils.rotations import attitude_matrix from astropy.table import Table from astropy.coordinates import SkyCoord import astropy.units as u ECSV = ( "/Users/adam/repos/roman_notebooks/notebooks/" "footprint_visualization/aux_data/roman_gps.sim.ecsv" ) print("Loading Roman SIAF …") rsiaf = pysiaf.Siaf("Roman") wfi_cen = rsiaf['WFI_CEN'] V2Ref = wfi_cen.V2Ref # +1546.38" (WFI focal-plane offset from V3 boresight) V3Ref = wfi_cen.V3Ref # -892.79" sensors = [rsiaf[f"WFI{j:02d}_FULL"] for j in range(1, 19)] print("Reading roman_gps.sim.ecsv …") t = Table.read(ECSV) # Roman WFI bandpass central wavelengths (µm) — for ordering _BP_WAVE = {'F062': 0.62, 'F087': 0.87, 'F106': 1.06, 'F129': 1.29, 'F146': 1.46, 'F158': 1.58, 'F184': 1.84, 'F213': 2.13, 'GRISM': 1.50, 'PRISM': 1.50} def unique_pointings(target_name): """Return sorted (ra, dec, pa) tuples, deduped on longest-wavelength bandpass.""" sub = t[t["TARGET_NAME"] == target_name] # keep only the longest-wavelength imaging bandpass available for this target bps = sorted(set(str(b) for b in sub["BANDPASS"]), key=lambda x: _BP_WAVE.get(x, 0.0)) sub = sub[sub["BANDPASS"] == bps[-1]] pts = {} for r in sub: key = (round(float(r["RA"]), 6), round(float(r["DEC"]), 6), round(float(r["PA"]), 3)) pts[key] = True return sorted(pts.keys()) def deduplicate_close_pointings(pointings, min_sep_arcmin=7.5): """ Greedy spatial deduplication: discard any pointing whose centre lies within min_sep_arcmin of a pointing already in the kept set. First occurrence wins. """ if len(pointings) <= 1: return list(pointings) ras = np.radians([p[0] for p in pointings]) decs = np.radians([p[1] for p in pointings]) sin_d = np.sin(decs) cos_d = np.cos(decs) thresh = np.radians(min_sep_arcmin / 60.0) kept = [0] for i in range(1, len(pointings)): cos_sep = (sin_d[i] * sin_d[kept] + cos_d[i] * cos_d[kept] * np.cos(ras[i] - ras[kept])) cos_sep = np.clip(cos_sep, -1.0, 1.0) if np.all(np.arccos(cos_sep) >= thresh): kept.append(i) return [pointings[i] for i in kept] def sca_polygons(ra_cen, dec_cen, pa): """Return list of 18 SCA polygons using corrected attitude_matrix.""" attmat = attitude_matrix(V2Ref, V3Ref, ra_cen, dec_cen, pa) polys = [] for sensor in sensors: sensor.set_attitude_matrix(attmat) c = sensor.corners("sky") # (ra_4pts, dec_4pts) polys.append([[round(float(c[0][i]), 5), round(float(c[1][i]), 5)] for i in range(4)]) return polys print("Computing Deep+Spec_W51 polygons …") _w51_all = unique_pointings("Deep+Spec_W51") w51_pts = deduplicate_close_pointings(_w51_all, min_sep_arcmin=7.5) print(f" {len(_w51_all)} → {len(w51_pts)} pointings after dedup (7.5\u2019 threshold) → {len(w51_pts)*18} polygons") print("Computing Disk7 polygons (may take ~30 s) …") _d7_all = unique_pointings("Disk7") d7_pts = deduplicate_close_pointings(_d7_all, min_sep_arcmin=7.5) print(f" {len(_d7_all)} → {len(d7_pts)} pointings after dedup (7.5\u2019 threshold) → {len(d7_pts)*18} polygons") PA_ROT = 45.0 # extra roll angle to show as a second overlay tiles = {} for i, (ra, dec, pa) in enumerate(w51_pts): coord = SkyCoord(ra=ra*u.deg, dec=dec*u.deg) gal = coord.galactic tiles[f"Deep+Spec_W51_{i+1:03d}"] = { "ra": ra, "dec": dec, "pa": pa, "l": round(float(gal.l.deg), 4), "b": round(float(gal.b.deg), 4), "polygons": sca_polygons(ra, dec, pa), } for i, (ra, dec, pa) in enumerate(d7_pts): if (i + 1) % 100 == 0: print(f" Disk7: {i+1}/{len(d7_pts)} …") coord = SkyCoord(ra=ra*u.deg, dec=dec*u.deg) gal = coord.galactic tiles[f"Disk7_{i+1:04d}"] = { "ra": ra, "dec": dec, "pa": pa, "l": round(float(gal.l.deg), 4), "b": round(float(gal.b.deg), 4), "polygons": sca_polygons(ra, dec, pa), } n_sca = sum(len(v['polygons']) for v in tiles.values()) print(f" {len(tiles)} tiles · {n_sca} SCA polygons") # PA-rotated footprints (nominal PA + PA_ROT) print(f"Computing PA+{PA_ROT:.0f}° rotated footprints …") tiles_rot = {} for key, tile in tiles.items(): tiles_rot[key] = { "ra": tile["ra"], "dec": tile["dec"], "l": tile["l"], "b": tile["b"], "polygons": sca_polygons(tile["ra"], tile["dec"], tile["pa"] + PA_ROT), } print(f" Done.") # Compact JSON for embedding rgps_js = json.dumps(tiles, separators=(',', ':')) rgps_rot_js = json.dumps(tiles_rot, separators=(',', ':')) html = f""" Roman RGPS Footprint — W51 — Aladin Lite

Roman RGPS — W51 Footprint

Background survey
Image overlay (click to toggle)
WFI footprint layers
Legend
Deep+Spec W51 (PA=0° cel.)
Disk7 (PA=0° cel.)
PA=+45° (celestial) overlay
Image overlay (active)
{len(w51_pts)} W51 · {len(d7_pts)} Disk7 pointings · 18 SCAs each
PA = V3-axis angle E of celestial North
At W51: galactic N ≈ 298° E of cel. N, so PA=0° appears
tilted ∼62° from galactic N in a galactic frame.
Source: roman_gps.sim.ecsv (APT)
Roman RGPS · pysiaf WFI geometry · {len(w51_pts)} W51 + {len(d7_pts)} Disk7 pointings · 18 SCAs each · F213 · roman_gps.sim.ecsv
""" with open("/Users/adam/Downloads/roman_footprint_rgps_w51.html", "w") as f: f.write(html) print(f"Written roman_footprint_rgps_w51.html") print(f" {len(tiles)} tiles, {sum(len(t['polygons']) for t in tiles.values())} total polygons")