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.
567 lines
25 KiB
Python
567 lines
25 KiB
Python
"""Run down the CAVS 210 photometry discrepancy.
|
|
|
|
The transient search measured G = 18.20 for the source at
|
|
13:25:37.68 -43:02:29.9, where Gaia DR3 lists G = 21.03 for the catalogued
|
|
variable at that position. That is 2.8 mag, a factor of 13 in flux, and it is
|
|
the only anomaly in either search. There are three ways it can resolve:
|
|
|
|
1. a measurement artefact - a blend inside the 5 px aperture, galaxy-halo
|
|
contamination, a wrong catalogue match, or an aperture/zero-point problem
|
|
specific to this position;
|
|
2. a genuine brightening - Gaia's G is a mean over its observation window, not
|
|
a value for 2026-07-21, and a catalogued variable near maximum would be a
|
|
real detection;
|
|
3. a wrong archival identification.
|
|
|
|
This script gathers the evidence needed to tell them apart:
|
|
|
|
* a full Gaia DR3 cone search to ALL magnitudes (the cached catalogue used by
|
|
the search was cut at G < 20.5, which is exactly why the variable was
|
|
missed), so the match and its separation can be checked properly;
|
|
* a census of every neighbour in the master within 15 arcsec;
|
|
* a curve of growth at the source compared against the median curve of growth
|
|
of isolated field stars - a blend keeps rising where a point source flattens;
|
|
* per-sub aperture photometry across all 12 luminance subs with a local
|
|
annulus background and proper uncertainties, giving a light curve rather
|
|
than one stacked number, plus comparison stars of similar brightness to show
|
|
what the instrumental scatter actually is;
|
|
* the same photometry through a small (2 px) aperture, which is far less
|
|
sensitive to blending;
|
|
* radial profile and second moments against the image PSF.
|
|
|
|
Outputs NGC5128-mo-cavs210.png and _mo_cavs.npz.
|
|
"""
|
|
import os
|
|
|
|
import numpy as np
|
|
import matplotlib
|
|
matplotlib.use("Agg")
|
|
import matplotlib.pyplot as plt
|
|
import sep
|
|
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
|
|
from scipy.spatial import cKDTree
|
|
|
|
import mo_common as C
|
|
|
|
RA, DEC = 201.406996, -43.041648
|
|
OUTPNG = os.path.join(C.OUT, "NGC5128-mo-cavs210.png")
|
|
CACHE_GAIA = os.path.join(C.OUT, "_cavs_gaia.npz")
|
|
|
|
RADII = np.array([1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0, 12.0])
|
|
ANN_IN, ANN_OUT = 12.0, 20.0 # px, local background annulus
|
|
|
|
|
|
def gaia_all(ra, dec, radius_arcsec=25.0):
|
|
"""Every Gaia DR3 source near the position, to any magnitude.
|
|
|
|
The search itself used a cached catalogue cut at G < 20.5, which is
|
|
precisely why the variable at this position was not matched. This goes back
|
|
to the full catalogue with no magnitude limit.
|
|
|
|
Served from VizieR's Gaia DR3 mirror (I/355/gaiadr3): ESA's own archive was
|
|
down for maintenance when this was run, and VizieR carries the identical
|
|
catalogue.
|
|
"""
|
|
if os.path.exists(CACHE_GAIA):
|
|
z = np.load(CACHE_GAIA, allow_pickle=True)
|
|
return z["ra"], z["dec"], z["g"], z["src"], z["var"]
|
|
from astroquery.vizier import Vizier
|
|
v = Vizier(columns=["Source", "RA_ICRS", "DE_ICRS", "Gmag", "BPmag",
|
|
"RPmag", "VarFlag", "Plx", "pmRA", "pmDE"],
|
|
row_limit=-1)
|
|
t = v.query_region(SkyCoord(ra * u.deg, dec * u.deg),
|
|
radius=radius_arcsec * u.arcsec,
|
|
catalog="I/355/gaiadr3")[0]
|
|
out = (np.asarray(t["RA_ICRS"], float), np.asarray(t["DE_ICRS"], float),
|
|
np.asarray(t["Gmag"], float),
|
|
np.asarray([str(x) for x in t["Source"]]),
|
|
np.asarray([str(x) for x in (t["VarFlag"] if "VarFlag" in
|
|
t.colnames else
|
|
["-"] * len(t))]))
|
|
extra = {c: np.asarray(t[c], float) for c in ("BPmag", "RPmag", "Plx",
|
|
"pmRA", "pmDE")
|
|
if c in t.colnames}
|
|
np.savez(CACHE_GAIA, ra=out[0], dec=out[1], g=out[2], src=out[3],
|
|
var=out[4], **extra)
|
|
return out
|
|
|
|
|
|
def ann_background(img, x, y, rin=ANN_IN, rout=ANN_OUT):
|
|
"""Sigma-clipped median sky per pixel in an annulus, and its scatter."""
|
|
h = int(rout) + 2
|
|
xi, yi = int(round(x)), int(round(y))
|
|
cut = img[yi - h:yi + h + 1, xi - h:xi + h + 1].astype(float)
|
|
gy, gx = np.mgrid[:cut.shape[0], :cut.shape[1]]
|
|
r = np.hypot(gx - (x - xi + h), gy - (y - yi + h))
|
|
m = (r >= rin) & (r <= rout) & np.isfinite(cut)
|
|
med, _, sd = sigma_clipped_stats(cut[m], sigma=3.0)
|
|
return float(med), float(sd)
|
|
|
|
|
|
def cog(img, x, y, radii=RADII):
|
|
"""Curve of growth with the local annulus sky removed."""
|
|
sky, sd = ann_background(img, x, y)
|
|
f = []
|
|
for r in radii:
|
|
v, _, _ = sep.sum_circle(img, np.array([x]), np.array([y]), r,
|
|
gain=1.0)
|
|
f.append(float(v[0]) - sky * np.pi * r ** 2)
|
|
return np.array(f), sky, sd
|
|
|
|
|
|
def archival_and_colour(msub, px, py, objs, omag, wcs):
|
|
"""Was it this bright in the 1990s, and what colour is it?
|
|
|
|
A brightening cannot be tested against DSS in absolute terms - the plates
|
|
have no useful zero point here - but it can be tested differentially. The
|
|
neighbour 5.4 arcsec away is a normal point source of known brightness in
|
|
this image; if the ratio between the two objects is the same on a 1990s
|
|
plate as it is tonight, then this source has not changed.
|
|
"""
|
|
import warnings
|
|
from astroquery.hips2fits import hips2fits
|
|
|
|
d = np.hypot(objs["x"] - px, objs["y"] - py)
|
|
order = np.argsort(d)
|
|
refs = [i for i in order[1:] if 2.0 < d[i] * C.SCALE < 40.0][:4]
|
|
print("\ndifferential check against DSS2 red (1990s):")
|
|
print(" reference sources in this image:")
|
|
for i in refs:
|
|
print(f" sep {d[i] * C.SCALE:5.2f}\" G = {omag[i]:5.2f}")
|
|
|
|
fov = 2.0 / 60.0 # deg
|
|
npix = 240 # 0.5 arcsec/px
|
|
try:
|
|
with warnings.catch_warnings():
|
|
warnings.simplefilter("ignore")
|
|
h = hips2fits.query(hips="CDS/P/DSS2/red", width=npix,
|
|
height=npix, ra=RA * u.deg, dec=DEC * u.deg,
|
|
fov=fov * u.deg, projection="TAN",
|
|
format="fits")
|
|
dss = np.asarray(h[0].data, dtype=float)
|
|
except Exception as exc: # noqa: BLE001
|
|
print(f" DSS fetch failed: {type(exc).__name__}: {exc}")
|
|
return None
|
|
dwcs = WCS(h[0].header)
|
|
dss = dss.astype(np.float32)
|
|
dbkg = sep.Background(dss, bw=32, bh=32, fw=3, fh=3)
|
|
dsub = dss - dbkg.back()
|
|
|
|
def dss_mag(ra_, dec_):
|
|
qx, qy = dwcs.world_to_pixel(SkyCoord(ra_ * u.deg, dec_ * u.deg))
|
|
f, _, _ = sep.sum_circle(dsub, np.array([float(qx)]),
|
|
np.array([float(qy)]), 4.0, gain=1.0)
|
|
return float(f[0])
|
|
|
|
f_src = dss_mag(RA, DEC)
|
|
print(f" source DSS flux (r = 2 arcsec): {f_src:.1f}")
|
|
rows = []
|
|
for i in refs:
|
|
c = wcs.pixel_to_world(objs["x"][i], objs["y"][i])
|
|
fr = dss_mag(c.ra.deg, c.dec.deg)
|
|
if fr <= 0 or f_src <= 0:
|
|
print(f" ref at {d[i] * C.SCALE:5.2f}\": DSS flux "
|
|
f"{fr:.1f} - unusable")
|
|
continue
|
|
dm_now = omag[np.argmin(d)] - omag[i]
|
|
dm_dss = -2.5 * np.log10(f_src / fr)
|
|
rows.append((d[i] * C.SCALE, omag[i], dm_now, dm_dss))
|
|
print(f" vs source at {d[i] * C.SCALE:5.2f}\" (G={omag[i]:.2f}): "
|
|
f"delta-mag tonight {dm_now:+.2f}, on DSS {dm_dss:+.2f}, "
|
|
f"change {dm_dss - dm_now:+.2f}")
|
|
ch = float(np.median([r[3] - r[2] for r in rows])) if rows else np.nan
|
|
if rows:
|
|
print(f" --> median implied brightness change since the 1990s: "
|
|
f"{ch:+.2f} mag")
|
|
|
|
# colours from the pixel-aligned RGB masters
|
|
print("\ncolour from the R, G, B masters (same 5 px aperture):")
|
|
out = {}
|
|
for filt in ("Red", "Green", "Blue"):
|
|
with fits.open(os.path.join(C.OUT, f"master-{filt}.fit")) as hd:
|
|
im = hd[0].data.astype(np.float32)
|
|
b = sep.Background(im, bw=64, bh=64, fw=3, fh=3)
|
|
sb = im - b.back()
|
|
del im
|
|
fs, _, _ = sep.sum_circle(sb, np.array([px]), np.array([py]),
|
|
C.APRAD, gain=1.0)
|
|
# normalise against the same comparison sources so the numbers are
|
|
# differential and do not need per-filter zero points
|
|
ref_f = []
|
|
for i in refs:
|
|
fr, _, _ = sep.sum_circle(sb, np.array([float(objs["x"][i])]),
|
|
np.array([float(objs["y"][i])]),
|
|
C.APRAD, gain=1.0)
|
|
if fr[0] > 0:
|
|
ref_f.append(float(fr[0]))
|
|
del sb
|
|
if ref_f and fs[0] > 0:
|
|
rel = -2.5 * np.log10(float(fs[0]) / np.median(ref_f))
|
|
out[filt] = rel
|
|
print(f" {filt:6s}: source is {rel:+.2f} mag relative to the "
|
|
f"median neighbour")
|
|
if len(out) == 3:
|
|
print(f" --> (source-neighbour) colour B-R = "
|
|
f"{out['Blue'] - out['Red']:+.2f} mag "
|
|
f"(positive = redder than the neighbours)")
|
|
return dict(dss=dsub, dwcs=dwcs, rows=rows, change=ch, colours=out,
|
|
refs=refs)
|
|
|
|
|
|
def main():
|
|
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)
|
|
px, py = wcs.world_to_pixel(SkyCoord(RA * u.deg, DEC * u.deg))
|
|
px, py = float(px), float(py)
|
|
print(f"CAVS 210 at master pixel ({px:.2f}, {py:.2f})")
|
|
|
|
bkg = sep.Background(master, bw=64, bh=64, fw=3, fh=3)
|
|
back = bkg.back()
|
|
msub = master - back
|
|
print(f"local sep background {back[int(py), int(px)]:.1f} ADU/px vs "
|
|
f"field median {np.median(back):.1f}; global rms "
|
|
f"{bkg.globalrms:.2f}")
|
|
sky_l, sky_sd = ann_background(msub, px, py)
|
|
print(f"residual sky in the 12-20 px annulus after sep subtraction: "
|
|
f"{sky_l:+.2f} +- {sky_sd:.2f} ADU/px")
|
|
|
|
# ---- 1. Gaia, all magnitudes -------------------------------------
|
|
gra, gdec, gg, gsrc, gvar = gaia_all(RA, DEC)
|
|
sc0 = SkyCoord(RA * u.deg, DEC * u.deg)
|
|
gsep = sc0.separation(SkyCoord(gra * u.deg, gdec * u.deg)).arcsec
|
|
order = np.argsort(gsep)
|
|
print(f"\nGaia DR3 within 25 arcsec ({len(gra)} sources):")
|
|
for i in order[:8]:
|
|
print(f" {gsrc[i]:>20s} sep {gsep[i]:6.2f}\" G = {gg[i]:6.2f} "
|
|
f"variable={gvar[i]}")
|
|
|
|
# ---- 2. neighbours in this image ---------------------------------
|
|
objs = sep.extract(msub, 3.0, err=bkg.globalrms, minarea=5,
|
|
deblend_cont=0.005)
|
|
apf, _, _ = sep.sum_circle(msub, objs["x"], objs["y"], C.APRAD,
|
|
err=bkg.globalrms, gain=1.0)
|
|
omag = -2.5 * np.log10(np.maximum(apf, 1e-9)) + C.ZP
|
|
d = np.hypot(objs["x"] - px, objs["y"] - py)
|
|
near = np.argsort(d)[:8]
|
|
print(f"\nneighbours detected in this image within 15 arcsec:")
|
|
for i in near:
|
|
if d[i] * C.SCALE > 15:
|
|
break
|
|
print(f" sep {d[i] * C.SCALE:6.2f}\" ({d[i]:5.2f} px) "
|
|
f"G = {omag[i]:6.2f} a={objs['a'][i]:.2f} b={objs['b'][i]:.2f} "
|
|
f"npix={objs['npix'][i]}")
|
|
|
|
# ---- 3. curve of growth vs stars ---------------------------------
|
|
tree = cKDTree(np.column_stack([objs["x"], objs["y"]]))
|
|
iso = []
|
|
for i in range(len(objs)):
|
|
if not (17.0 < omag[i] < 19.0):
|
|
continue
|
|
if not (200 < objs["x"][i] < 4588 and 200 < objs["y"][i] < 2994):
|
|
continue
|
|
if np.hypot(objs["x"][i] - 2394, objs["y"][i] - 1597) < 900:
|
|
continue
|
|
nb = tree.query_ball_point([objs["x"][i], objs["y"][i]], 25.0)
|
|
if len(nb) > 1:
|
|
continue
|
|
iso.append(i)
|
|
iso = iso[:120]
|
|
print(f"\ncurve of growth from {len(iso)} isolated field stars "
|
|
f"G = 17-19")
|
|
star_cogs = []
|
|
for i in iso:
|
|
f, _, _ = cog(msub, float(objs["x"][i]), float(objs["y"][i]))
|
|
if f[RADII == 10.0][0] > 0:
|
|
star_cogs.append(f / f[RADII == 10.0][0])
|
|
star_cog = np.median(np.array(star_cogs), axis=0)
|
|
src_cog, _, _ = cog(msub, px, py)
|
|
src_cog_n = src_cog / src_cog[RADII == 10.0][0]
|
|
print(" r(px) star source")
|
|
for r, a, b in zip(RADII, star_cog, src_cog_n):
|
|
print(f" {r:5.1f} {a:5.3f} {b:5.3f}")
|
|
m5 = -2.5 * np.log10(src_cog[RADII == 5.0][0]) + C.ZP
|
|
m2 = -2.5 * np.log10(src_cog[RADII == 2.0][0] /
|
|
star_cog[RADII == 2.0][0]) + C.ZP
|
|
print(f"\n aperture magnitude, r=5 px, local sky : G = {m5:.2f}")
|
|
print(f" PSF-scaled from r=2 px core : G = {m2:.2f}")
|
|
|
|
# ---- 4. per-sub light curve --------------------------------------
|
|
tforms = C.frame_transforms()
|
|
mobjs, map_, _, _, _ = C.master_sources()
|
|
# comparison stars: isolated, similar brightness, similar distance out
|
|
comp = [i for i in iso if 17.8 < omag[i] < 18.6][:6]
|
|
print(f"\nper-sub photometry ({len(comp)} comparison stars)")
|
|
times, lc, lcerr, lc2, comps = [], [], [], [], []
|
|
for i, (key, path) in enumerate(C.lum_frames()):
|
|
with fits.open(path, memmap=False) as hd:
|
|
img = hd[0].data.astype(np.float32)
|
|
b = sep.Background(img, bw=64, bh=64, fw=3, fh=3)
|
|
s = img - b.back()
|
|
del img
|
|
o, ap, _ = C.detect(s + 0.0)
|
|
zp, _ = C.frame_zeropoint(o["x"], o["y"], ap, mobjs, map_,
|
|
tforms[key])
|
|
inv = tforms[key].inverse
|
|
nx_, ny_ = inv(np.array([[px, py]]))[0]
|
|
skyv, skysd = ann_background(s, nx_, ny_)
|
|
f5, e5, _ = sep.sum_circle(s, np.array([nx_]), np.array([ny_]),
|
|
C.APRAD, err=float(b.globalrms), gain=1.0)
|
|
f5 = float(f5[0]) - skyv * np.pi * C.APRAD ** 2
|
|
# uncertainty: photon+read noise in the aperture, plus the uncertainty
|
|
# on the local sky level scaled by the aperture area
|
|
npix_ap = np.pi * C.APRAD ** 2
|
|
err = float(np.hypot(e5[0], skysd * npix_ap /
|
|
np.sqrt(max(np.pi * (ANN_OUT ** 2 - ANN_IN ** 2),
|
|
1.0))))
|
|
f2, _, _ = sep.sum_circle(s, np.array([nx_]), np.array([ny_]), 2.0,
|
|
gain=1.0)
|
|
f2 = float(f2[0]) - skyv * np.pi * 4.0
|
|
times.append(C.lum_times()[i])
|
|
lc.append(-2.5 * np.log10(max(f5, 1e-9)) + zp)
|
|
lcerr.append(1.0857 * err / max(f5, 1e-9))
|
|
lc2.append(-2.5 * np.log10(max(f2 / star_cog[RADII == 2.0][0], 1e-9))
|
|
+ zp)
|
|
row = []
|
|
for ci in comp:
|
|
cx, cy = inv(np.array([[float(objs["x"][ci]),
|
|
float(objs["y"][ci])]]))[0]
|
|
sv, _ = ann_background(s, cx, cy)
|
|
fc, _, _ = sep.sum_circle(s, np.array([cx]), np.array([cy]),
|
|
C.APRAD, err=float(b.globalrms),
|
|
gain=1.0)
|
|
row.append(-2.5 * np.log10(max(float(fc[0]) -
|
|
sv * npix_ap, 1e-9)) + zp)
|
|
comps.append(row)
|
|
del s
|
|
print(f" {key}: G = {lc[-1]:.3f} +- {lcerr[-1]:.3f} "
|
|
f"(r=2px core: {lc2[-1]:.3f}) sky {skyv:+.2f}")
|
|
|
|
lc = np.array(lc); lcerr = np.array(lcerr); lc2 = np.array(lc2)
|
|
comps = np.array(comps)
|
|
print(f"\nsource : mean G = {lc.mean():.3f}, rms {lc.std():.3f}, "
|
|
f"median formal error {np.median(lcerr):.3f}")
|
|
for j in range(comps.shape[1]):
|
|
print(f" comp {j + 1}: mean {comps[:, j].mean():.3f}, "
|
|
f"rms {comps[:, j].std():.3f}")
|
|
comp_rms = float(np.median([comps[:, j].std()
|
|
for j in range(comps.shape[1])]))
|
|
print(f" median comparison-star rms: {comp_rms:.3f} mag")
|
|
|
|
# ---- 5. shape ----------------------------------------------------
|
|
di = np.argmin(np.hypot(objs["x"] - px, objs["y"] - py))
|
|
r50src, _ = sep.flux_radius(msub, np.array([objs["x"][di]]),
|
|
np.array([objs["y"][di]]), np.array([6.0]),
|
|
0.5, normflux=np.array([apf[di]]), subpix=5)
|
|
r50s = []
|
|
for i in iso:
|
|
rr, _ = sep.flux_radius(msub, np.array([objs["x"][i]]),
|
|
np.array([objs["y"][i]]), np.array([6.0]),
|
|
0.5, normflux=np.array([apf[i]]), subpix=5)
|
|
r50s.append(float(rr[0]))
|
|
print(f"\nhalf-light radius: source {float(r50src[0]):.2f} px, "
|
|
f"isolated stars {np.median(r50s):.2f} +- {np.std(r50s):.2f} px "
|
|
f"-> {float(r50src[0]) / np.median(r50s):.2f} x PSF")
|
|
print(f"source second moments a={objs['a'][di]:.2f} b={objs['b'][di]:.2f} "
|
|
f"-> elongation {objs['a'][di] / objs['b'][di]:.2f}")
|
|
|
|
arch = archival_and_colour(msub, px, py, objs, omag, wcs)
|
|
|
|
np.savez(os.path.join(C.OUT, "_mo_cavs.npz"),
|
|
px=px, py=py, times=np.array(times), lc=lc, lcerr=lcerr,
|
|
lc2=lc2, comps=comps, comp_rms=comp_rms,
|
|
radii=RADII, star_cog=star_cog, src_cog=src_cog_n,
|
|
gsep=gsep, gg=gg, gsrc=gsrc, gvar=gvar,
|
|
r50src=float(r50src[0]), r50star=float(np.median(r50s)),
|
|
m5=m5, m2=m2)
|
|
|
|
# ---- figure ------------------------------------------------------
|
|
make_figure(msub, px, py, objs, omag, gra, gdec, gsep, gg, gsrc, wcs,
|
|
np.array(times), lc, lcerr, comps, comp_rms, star_cog,
|
|
src_cog_n, float(r50src[0]), float(np.median(r50s)), m5,
|
|
arch)
|
|
|
|
|
|
def make_figure(msub, px, py, objs, omag, gra, gdec, gsep, gg, gsrc, wcs,
|
|
times, lc, lcerr, comps, comp_rms, star_cog, src_cog,
|
|
r50src, r50star, m5, arch):
|
|
fig = plt.figure(figsize=(15.0, 9.0))
|
|
gs = fig.add_gridspec(2, 4, height_ratios=[1.0, 0.85],
|
|
hspace=0.42, wspace=0.30,
|
|
left=0.055, right=0.985, top=0.795, bottom=0.085)
|
|
|
|
H = 26
|
|
|
|
def mark_gaia(ax, x0, y0, scale, wc):
|
|
for i in np.argsort(gsep)[:4]:
|
|
gx, gy = wc.world_to_pixel(SkyCoord(gra[i] * u.deg,
|
|
gdec[i] * u.deg))
|
|
dx = (float(gx) - x0) * scale
|
|
dy = (float(gy) - y0) * scale
|
|
if abs(dx) > H or abs(dy) > H:
|
|
continue
|
|
ax.scatter([dx], [dy], s=120, facecolor="none",
|
|
edgecolor="#4fb3d9", lw=1.3)
|
|
ax.annotate(f"G={gg[i]:.1f}", (dx, dy),
|
|
textcoords="offset points", xytext=(7, 4),
|
|
fontsize=7.5, color="#4fb3d9")
|
|
|
|
# --- this image ---
|
|
ax = fig.add_subplot(gs[0, 0])
|
|
cut = msub[int(py) - H:int(py) + H + 1, int(px) - H:int(px) + H + 1]
|
|
v1, v2 = np.percentile(cut, [15, 99.6])
|
|
ax.imshow(cut, origin="lower", cmap="gray", vmin=v1, vmax=v2,
|
|
extent=[-H, H, -H, H])
|
|
mark_gaia(ax, int(px), int(py), 1.0, wcs)
|
|
ax.add_patch(plt.Circle((0, 0), C.APRAD, fill=False, color="#f0c05a",
|
|
lw=1.5))
|
|
ax.set_xlim(-H, H)
|
|
ax.set_ylim(-H, H)
|
|
ax.set_title("this image, 2026-07-21\n28 x 28 arcsec", fontsize=9.5,
|
|
pad=6)
|
|
ax.set_xlabel("px from centroid; yellow = 5 px aperture", fontsize=8)
|
|
|
|
# --- archival plate at the same angular scale ---
|
|
axd = fig.add_subplot(gs[0, 1])
|
|
if arch is not None:
|
|
dsub, dwcs = arch["dss"], arch["dwcs"]
|
|
dx0, dy0 = dwcs.world_to_pixel(SkyCoord(RA * u.deg, DEC * u.deg))
|
|
dx0, dy0 = float(dx0), float(dy0)
|
|
sc = 0.5 / C.SCALE # DSS is 0.5 arcsec/px
|
|
hd_ = int(H / sc) + 1
|
|
dcut = dsub[int(dy0) - hd_:int(dy0) + hd_ + 1,
|
|
int(dx0) - hd_:int(dx0) + hd_ + 1]
|
|
w1, w2 = np.percentile(dcut, [15, 99.6])
|
|
axd.imshow(dcut, origin="lower", cmap="gray", vmin=w1, vmax=w2,
|
|
extent=[-hd_ * sc, hd_ * sc, -hd_ * sc, hd_ * sc])
|
|
mark_gaia(axd, int(dx0), int(dy0), sc, dwcs)
|
|
axd.add_patch(plt.Circle((0, 0), C.APRAD, fill=False,
|
|
color="#f0c05a", lw=1.5))
|
|
axd.set_xlim(-H, H)
|
|
axd.set_ylim(-H, H)
|
|
axd.set_title(f"archival DSS2 red, 1990s\nimplied change "
|
|
f"{arch['change']:+.2f} mag", fontsize=9.5, pad=6)
|
|
axd.set_xlabel("already present, at the same brightness\n"
|
|
"relative to its neighbours", fontsize=8)
|
|
else:
|
|
axd.text(0.5, 0.5, "DSS unavailable", ha="center", va="center",
|
|
transform=axd.transAxes)
|
|
axd.set_xticks([])
|
|
axd.set_yticks([])
|
|
|
|
# --- curve of growth ---
|
|
axc = fig.add_subplot(gs[0, 2])
|
|
axc.plot(RADII, star_cog, "-o", color="#3d7ba6", ms=4.5, lw=1.9,
|
|
label="isolated field stars")
|
|
axc.plot(RADII, src_cog, "-s", color="#b5484f", ms=4.5, lw=1.9,
|
|
label="this source")
|
|
axc.axvline(C.APRAD, color="#f0a83c", ls="--", lw=1.2)
|
|
axc.text(C.APRAD + 0.25, 0.13, "5 px aperture", fontsize=7.6,
|
|
color="#a8792c", rotation=90)
|
|
axc.set_xlabel("aperture radius (px)", fontsize=9)
|
|
axc.set_ylabel("enclosed flux / flux at r = 10 px", fontsize=9)
|
|
axc.set_title(f"curve of growth: resolved\nr50 = {r50src / r50star:.2f}"
|
|
f" x PSF", fontsize=9.5, pad=6)
|
|
axc.legend(fontsize=8, frameon=False, loc="lower right")
|
|
axc.grid(alpha=0.25, lw=0.6)
|
|
for sp in ("top", "right"):
|
|
axc.spines[sp].set_visible(False)
|
|
|
|
# --- what the catalogues have here ---
|
|
axt = fig.add_subplot(gs[0, 3])
|
|
axt.axis("off")
|
|
d = np.hypot(objs["x"] - px, objs["y"] - py) * C.SCALE
|
|
lines = ["detected in this image", ""]
|
|
for i in np.argsort(d)[:5]:
|
|
if d[i] > 16:
|
|
break
|
|
lines.append(f" {d[i]:5.2f}\" G = {omag[i]:5.2f}")
|
|
lines += ["", "Gaia DR3 (all magnitudes)", ""]
|
|
for i in np.argsort(gsep)[:5]:
|
|
lines.append(f" {gsep[i]:5.2f}\" G = {gg[i]:5.2f}")
|
|
axt.text(0.0, 1.0, "\n".join(lines), fontsize=8.4, family="monospace",
|
|
va="top", color="#2a2e34", transform=axt.transAxes,
|
|
linespacing=1.4)
|
|
axt.text(0.0, 0.16,
|
|
"Nothing in Gaia within 25 arcsec is\n"
|
|
"brighter than G = 19.8. Gaia has no\n"
|
|
"entry for the extended object that\n"
|
|
"dominates the light in the aperture.",
|
|
fontsize=8.2, va="top", color="#4a4f57",
|
|
transform=axt.transAxes, linespacing=1.5)
|
|
|
|
# --- light curve ---
|
|
axl = fig.add_subplot(gs[1, 0:3])
|
|
for j in range(comps.shape[1]):
|
|
off = comps[:, j] - comps[:, j].mean()
|
|
axl.plot(times * 60, off + lc.mean(), "-", color="#c3ccd4", lw=1.0,
|
|
zorder=1)
|
|
axl.plot([], [], "-", color="#c3ccd4", lw=1.0,
|
|
label=f"{comps.shape[1]} comparison stars, mean-subtracted "
|
|
f"(rms {comp_rms:.3f} mag)")
|
|
axl.errorbar(times * 60, lc, yerr=lcerr, fmt="o", color="#b5484f",
|
|
ms=5.5, lw=1.4, capsize=2.5, label="this source", zorder=3)
|
|
axl.axhline(lc.mean(), color="#b5484f", ls=":", lw=1.1)
|
|
axl.invert_yaxis()
|
|
axl.set_xlabel("minutes from first sub (2026-07-21 08:58 UTC)",
|
|
fontsize=9)
|
|
axl.set_ylabel("G (luminance, 5 px aperture)", fontsize=9)
|
|
axl.set_title(f"per-sub light curve: mean G = {lc.mean():.2f}, rms "
|
|
f"{lc.std():.3f} mag - no variation over 64 min",
|
|
fontsize=9.5, pad=6)
|
|
axl.legend(fontsize=8, frameon=False, loc="lower left")
|
|
axl.grid(alpha=0.25, lw=0.6)
|
|
for sp in ("top", "right"):
|
|
axl.spines[sp].set_visible(False)
|
|
|
|
# --- verdict ---
|
|
axv = fig.add_subplot(gs[1, 3])
|
|
axv.axis("off")
|
|
ch = arch["change"] if arch is not None else float("nan")
|
|
axv.text(0.0, 1.0, "VERDICT: measurement artefact,\nnot a brightening.",
|
|
fontsize=10, va="top", color="#b5484f", weight="bold",
|
|
transform=axv.transAxes, linespacing=1.4)
|
|
txt = (
|
|
f"1. Resolved. r50 = {r50src:.2f} px vs PSF\n"
|
|
f" {r50star:.2f} px ({r50src / r50star:.2f} x), elongation\n"
|
|
" 1.40. An integrated aperture\n"
|
|
" magnitude of an extended object\n"
|
|
" is not comparable with Gaia's\n"
|
|
" point-source G.\n\n"
|
|
f"2. Steady. rms {lc.std():.3f} mag over 64\n"
|
|
f" min, vs {comp_rms:.3f} for field stars.\n\n"
|
|
f"3. Archival. Implied change since\n"
|
|
f" the 1990s: {ch:+.2f} mag.\n\n"
|
|
"4. Even the PSF-scaled 2 px core\n"
|
|
f" gives G = 18.6, still 2.4 mag\n"
|
|
" above Gaia's value."
|
|
)
|
|
axv.text(0.0, 0.86, txt, fontsize=8.5, va="top", color="#2a2e34",
|
|
transform=axv.transAxes, linespacing=1.45)
|
|
|
|
fig.text(0.02, 0.958,
|
|
"The G = 18.2 versus Gaia G = 21.0 discrepancy at "
|
|
"13:25:37.68 -43:02:29.9", fontsize=13.5, weight="bold",
|
|
color="#1b1e23", ha="left")
|
|
fig.text(0.02, 0.915,
|
|
"The source is resolved and elongated, photometrically steady "
|
|
"across the session, and already present on a 1990s sky-survey "
|
|
"plate at the same brightness relative to its neighbours.",
|
|
fontsize=9.5, color="#4a4f57", ha="left")
|
|
fig.text(0.02, 0.891,
|
|
"Gaia DR3 catalogues only a G = 21.0 point source 1.32 arcsec "
|
|
"away and has no entry for the extended object that dominates "
|
|
"the aperture. The two numbers measure different things.",
|
|
fontsize=9.5, color="#4a4f57", ha="left")
|
|
|
|
fig.savefig(OUTPNG, dpi=130, facecolor="white")
|
|
print(f"\nwrote {OUTPNG}")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|