laurence
c6299f41ab
Preparing to merge this repository into a combined astrophotography repo. session-scripts/ becomes pipeline/ because the scripts import layout.py from their own directory and must stay together, and because 'pipeline' says what it is rather than how it came about. observing/ stays at the top level: observing plans are not processing code.
228 lines
9.9 KiB
Python
228 lines
9.9 KiB
Python
"""Figure: the transient search and why its six finalists are not transients.
|
|
|
|
Top block, one column per finalist: the 60 min luminance master beside a
|
|
Digitized Sky Survey red plate of the same patch of sky from the 1990s. If a
|
|
source is on both, it is not a transient, and that is a judgement a reader can
|
|
make from the picture without trusting any of the catalogue machinery.
|
|
|
|
Bottom left: where the catalogues run out. Almost every source brighter than
|
|
G = 19 has a Gaia DR3 counterpart; below that the unmatched fraction climbs
|
|
steeply, not because transients appear but because Gaia's G < 20.5 point-source
|
|
catalogue stops being complete against a bright galaxy. This is what sets the
|
|
honest limit of the search.
|
|
|
|
Bottom middle: the point-source selection. Half-light radius against magnitude,
|
|
with the image's own PSF measured from Gaia stars. Note that Centaurus A's
|
|
globular clusters sit inside the stellar locus - at 3.8 Mpc they are unresolved
|
|
in 2.7 arcsec seeing - so morphology cannot separate a cluster from a
|
|
supernova here, and five of the six finalists are clusters.
|
|
|
|
Bottom right: where the finalists are, on the master.
|
|
"""
|
|
import os
|
|
|
|
import numpy as np
|
|
import matplotlib
|
|
matplotlib.use("Agg")
|
|
import matplotlib.pyplot as plt
|
|
from astropy import units as u
|
|
from astropy.coordinates import SkyCoord
|
|
from astropy.io import fits
|
|
from astropy.stats import sigma_clipped_stats
|
|
from astropy.wcs import WCS
|
|
|
|
import mo_common as C
|
|
|
|
OUTPNG = os.path.join(C.OUT, "NGC5128-mo-transient-search.png")
|
|
HALF = 26 # px in the master = 28 arcsec across
|
|
|
|
# Identifications established by cross-matching the six catalogue non-matches
|
|
# against published Centaurus A cluster catalogues, SIMBAD and the VizieR
|
|
# archive as a whole (see notes-transient-search.md).
|
|
IDENT = {
|
|
"201.672625": ("WHH-31", "Cen A globular cluster", "V=19.50 (Woodley+07)"),
|
|
"201.500583": ("pff_gc-081", "Cen A globular cluster",
|
|
"V=19.34 (Woodley+07)"),
|
|
"201.441625": ("WHH-24", "Cen A globular cluster", "V=19.63 (Woodley+07)"),
|
|
"201.407000": ("at CAVS 210",
|
|
"RESOLVED (1.34x PSF), archival 1996-2020",
|
|
"DENIS I=18.23; Gaia pt src G=21.03 at 1.32\" - see 4.3"),
|
|
"201.303542": ("WHH-10", "Cen A globular cluster", "V=19.57 (Woodley+07)"),
|
|
"201.256708": ("pff_gc-029", "Cen A globular cluster",
|
|
"V=19.37 (Woodley+07)"),
|
|
}
|
|
|
|
|
|
def ident_for(ra):
|
|
for k, v in IDENT.items():
|
|
if abs(float(k) - ra) < 1e-4:
|
|
return v
|
|
return ("unidentified", "", "")
|
|
|
|
|
|
def stretch(a, lo=25.0, hi=99.6):
|
|
a = np.asarray(a, dtype=float)
|
|
good = np.isfinite(a)
|
|
if not good.any():
|
|
return 0.0, 1.0
|
|
return np.percentile(a[good], lo), np.percentile(a[good], hi)
|
|
|
|
|
|
def main():
|
|
rows = list(np.load(os.path.join(C.OUT, "_mo_vetted.npy"),
|
|
allow_pickle=True))
|
|
rows.sort(key=lambda r: r["dist_cenA_arcmin"])
|
|
z = np.load(os.path.join(C.OUT, "_mo_transient.npz"))
|
|
|
|
with fits.open(os.path.join(C.OUT, "master-Luminance.fit")) as hd:
|
|
master = hd[0].data.astype(np.float32)
|
|
hdr = hd[0].header
|
|
wcs = WCS(hdr)
|
|
|
|
import sep
|
|
bkg = sep.Background(master, bw=64, bh=64, fw=3, fh=3)
|
|
msub = master - bkg.back()
|
|
|
|
n = len(rows)
|
|
fig = plt.figure(figsize=(15.0, 11.6))
|
|
gs = fig.add_gridspec(3, n, height_ratios=[1.0, 1.0, 1.55],
|
|
hspace=0.30, wspace=0.09,
|
|
left=0.055, right=0.985, top=0.845, bottom=0.075)
|
|
|
|
for i, r in enumerate(rows):
|
|
name, kind, phot = ident_for(r["ra_deg"])
|
|
x, y = int(round(r["x"])), int(round(r["y"]))
|
|
cut = msub[y - HALF:y + HALF + 1, x - HALF:x + HALF + 1]
|
|
ax = fig.add_subplot(gs[0, i])
|
|
vlo, vhi = stretch(cut, 20, 99.5)
|
|
ax.imshow(cut, origin="lower", cmap="gray", vmin=vlo, vmax=vhi)
|
|
ax.plot([HALF, HALF], [HALF + 6, HALF + 11], color="#f0c05a", lw=1.1)
|
|
ax.plot([HALF + 6, HALF + 11], [HALF, HALF], color="#f0c05a", lw=1.1)
|
|
ax.set_xticks([]); ax.set_yticks([])
|
|
for sp in ax.spines.values():
|
|
sp.set_color("#2f7d54" if kind.startswith("Cen") else "#b5484f")
|
|
sp.set_linewidth(1.4)
|
|
ax.set_title(f"{name}\n{r['ra_hms']} {r['dec_dms']}", fontsize=8.6,
|
|
color="#1b1e23", pad=5)
|
|
if i == 0:
|
|
ax.set_ylabel("this image\n60 min luminance", fontsize=9,
|
|
color="#3b3f45")
|
|
|
|
# archival plate
|
|
tag = f"{r['ra_deg']:.5f}{r['dec_deg']:+.5f}_mred.fits"
|
|
path = os.path.join(C.OUT, "_dss", tag)
|
|
ax2 = fig.add_subplot(gs[1, i])
|
|
if os.path.exists(path):
|
|
with fits.open(path) as h2:
|
|
d = np.asarray(h2[0].data, dtype=float)
|
|
v1, v2 = stretch(d, 20, 99.5)
|
|
ax2.imshow(d, origin="lower", cmap="gray", vmin=v1, vmax=v2)
|
|
m = d.shape[0] / 2.0
|
|
f = d.shape[0] / (2.0 * HALF + 1)
|
|
ax2.plot([m, m], [m + 6 * f, m + 11 * f], color="#f0c05a", lw=1.1)
|
|
ax2.plot([m + 6 * f, m + 11 * f], [m, m], color="#f0c05a", lw=1.1)
|
|
else:
|
|
ax2.text(0.5, 0.5, "no DSS cutout", ha="center", va="center",
|
|
fontsize=8, color="#8a9099", transform=ax2.transAxes)
|
|
ax2.set_xticks([]); ax2.set_yticks([])
|
|
for sp in ax2.spines.values():
|
|
sp.set_color("#9aa0a8")
|
|
sp.set_linewidth(0.9)
|
|
ax2.set_xlabel(f"{kind}\n{phot}\nG(this image) = {r['g_mag']:.2f}, "
|
|
f"{r['dist_cenA_arcmin']:.1f}' from nucleus",
|
|
fontsize=7.8, color="#4a4f57", labelpad=5)
|
|
if i == 0:
|
|
ax2.set_ylabel("archival\nDSS2 red, 1990s", fontsize=9,
|
|
color="#3b3f45")
|
|
|
|
# ---------------- magnitude completeness ----------------
|
|
axm = fig.add_subplot(gs[2, 0:2])
|
|
good = z["good"] & z["pointlike"]
|
|
mag = z["mag"][good]
|
|
matched = z["gaia_ok"][good]
|
|
bins = np.arange(14.0, 21.01, 0.25)
|
|
axm.hist([mag[matched], mag[~matched]], bins=bins, stacked=True,
|
|
color=["#3d7ba6", "#b5484f"],
|
|
label=["matched to Gaia DR3", "no Gaia counterpart"])
|
|
axm.axvline(20.5, color="#4a4f57", ls="--", lw=1.2)
|
|
axm.text(20.42, axm.get_ylim()[1] * 0.92, "Gaia cat.\nlimit G=20.5",
|
|
fontsize=7.6, color="#4a4f57", ha="right", va="top")
|
|
axm.set_xlabel("G magnitude (5 px aperture, master)", fontsize=9.5)
|
|
axm.set_ylabel("point-like sources per 0.25 mag", fontsize=9.5)
|
|
axm.set_title("where the reference catalogue runs out", fontsize=10.5,
|
|
color="#1b1e23", pad=7)
|
|
axm.legend(fontsize=8.2, frameon=False, loc="upper left")
|
|
axm.grid(alpha=0.22, lw=0.6)
|
|
for sp in ("top", "right"):
|
|
axm.spines[sp].set_visible(False)
|
|
|
|
# ---------------- morphology ----------------
|
|
axr = fig.add_subplot(gs[2, 2:4])
|
|
sel = z["good"]
|
|
axr.scatter(z["mag"][sel], z["r50"][sel] / z["r50_star"], s=2.5,
|
|
c="#9fb6c6", alpha=0.35, lw=0, label="all detections")
|
|
ci = z["cand_idx"]
|
|
axr.scatter(z["mag"][ci], z["r50"][ci] / z["r50_star"], s=9,
|
|
c="#e0a13c", alpha=0.8, lw=0,
|
|
label="no Gaia / SkyMapper match")
|
|
for r in rows:
|
|
axr.scatter([r["g_mag"]], [r["r50_over_psf"]], s=64,
|
|
facecolor="none",
|
|
edgecolor="#2f7d54" if ident_for(r["ra_deg"])[1]
|
|
.startswith("Cen") else "#b5484f", lw=1.5, zorder=5)
|
|
axr.axhline(1.0, color="#4a4f57", ls=":", lw=1.1)
|
|
axr.text(14.3, 1.03, "PSF (Gaia stars)", fontsize=7.6, color="#4a4f57")
|
|
axr.set_ylim(0.4, 3.2)
|
|
axr.set_xlim(14.0, 21.0)
|
|
axr.set_xlabel("G magnitude", fontsize=9.5)
|
|
axr.set_ylabel("half-light radius / PSF", fontsize=9.5)
|
|
axr.set_title("morphology cannot separate clusters from transients",
|
|
fontsize=10.5, color="#1b1e23", pad=7)
|
|
axr.legend(fontsize=8.0, frameon=False, loc="upper left",
|
|
markerscale=2.2)
|
|
axr.grid(alpha=0.22, lw=0.6)
|
|
for sp in ("top", "right"):
|
|
axr.spines[sp].set_visible(False)
|
|
|
|
# ---------------- field map ----------------
|
|
axf = fig.add_subplot(gs[2, 4:6])
|
|
sm = master[::6, ::6]
|
|
v1, v2 = np.percentile(sm[np.isfinite(sm)], [12, 99.75])
|
|
axf.imshow(np.arcsinh((sm - v1) / max(v2 - v1, 1e-6)), origin="lower",
|
|
cmap="gray", vmin=0, vmax=np.arcsinh(1.0))
|
|
for r in rows:
|
|
c = "#2f7d54" if ident_for(r["ra_deg"])[1].startswith("Cen") \
|
|
else "#b5484f"
|
|
axf.scatter([r["x"] / 6], [r["y"] / 6], s=90, facecolor="none",
|
|
edgecolor=c, lw=1.5)
|
|
axf.annotate(ident_for(r["ra_deg"])[0], (r["x"] / 6, r["y"] / 6),
|
|
textcoords="offset points", xytext=(9, 5), fontsize=7.4,
|
|
color=c)
|
|
axf.set_xticks([]); axf.set_yticks([])
|
|
axf.set_title("positions on the field (43' x 29')", fontsize=10.5,
|
|
color="#1b1e23", pad=7)
|
|
for sp in axf.spines.values():
|
|
sp.set_color("#9aa0a8")
|
|
|
|
fig.text(0.02, 0.962,
|
|
"Transient search in the 60 min luminance master of NGC 5128 "
|
|
"(Centaurus A), 2026-07-21", fontsize=14, weight="bold",
|
|
color="#1b1e23", ha="left")
|
|
fig.text(0.02, 0.936,
|
|
"7,519 point-like sources detected above SNR 10. 603 had no "
|
|
"Gaia DR3 counterpart; 567 also had none in SkyMapper DR2; 6 of "
|
|
"those survived the requirement of an independent detection in "
|
|
"at least 6 of the 12 subs.",
|
|
fontsize=9.4, color="#4a4f57", ha="left")
|
|
fig.text(0.02, 0.917,
|
|
"All 6 are identified below. Five are catalogued Centaurus A "
|
|
"globular clusters; the sixth is a resolved object with archival "
|
|
"detections from 1996 to 2020. No transient candidate remains.",
|
|
fontsize=9.4, color="#4a4f57", ha="left")
|
|
|
|
fig.savefig(OUTPNG, dpi=130, facecolor="white")
|
|
print(f"wrote {OUTPNG}")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|