""" gc-complete.py -- step 3. Artificial star tests. Why this exists: NGC 5128's light makes the noise rise steeply toward the nucleus, so the SAME cluster is harder to detect at r=2' than at r=20'. Without correcting for that, the measured radial density profile is the true profile MULTIPLIED by an unknown, radially-decreasing completeness -- which flattens or even inverts a real central concentration. So we measure the completeness directly: inject fake point sources of known magnitude at known positions into the real luminance image, re-run the identical detection and selection chain, and count how many come back, as a function of magnitude and projected radius. The PSF is empirical, built by median-stacking bright isolated stars from the image itself. Output: _gc_complete.npz (recovery grid), and the PSF stamp. """ import os, time, numpy as np, sep, warnings from astropy.io import fits from astropy.wcs import WCS import layout warnings.filterwarnings("ignore") S = layout.SESSION BW, FW, THRESH, MINAREA, APR = 16, 3, 3.0, 5, 5.0 ZP_L = 27.941 PIXSCALE = 0.5376 NUC_RA, NUC_DEC = 201.365063, -43.019113 NRUN = 6 # injection runs NPER = 2500 # fakes per run MAGS = (17.5, 22.5) HALF = 15 # PSF stamp half-size, px RNG = np.random.default_rng(20260721) def gk(fwhm=5.0): sig = fwhm / 2.3548 n = int(2 * round(3 * sig) + 1) r = np.arange(n) - n // 2 xx, yy = np.meshgrid(r, r) k = np.exp(-(xx ** 2 + yy ** 2) / (2 * sig ** 2)) return (k / k.sum()).astype(np.float32) def build_psf(data, cat): """median stack of bright, isolated, unsaturated stars""" m = (cat["base"] & cat["is_star"] & (cat["mag"] > 14.5) & (cat["mag"] < 17.0) & (cat["r_arcmin"] > 8) & (cat["fwhm"] < 5.6)) x, y = cat["x"][m], cat["y"][m] # isolation: no other detection within 20 px ax, ay = cat["x"], cat["y"] keep = [] for i in range(len(x)): dd = np.hypot(ax - x[i], ay - y[i]) if (dd < 20).sum() <= 1: keep.append(i) x, y = x[keep], y[keep] print(" PSF from %d isolated bright stars" % len(x)) stamps = [] for xi, yi in zip(x, y): ix, iy = int(round(xi)), int(round(yi)) if ix < HALF or iy < HALF or ix >= data.shape[1] - HALF or iy >= data.shape[0] - HALF: continue st = data[iy - HALF:iy + HALF + 1, ix - HALF:ix + HALF + 1].astype(np.float64) st = st - np.median(st[[0, -1], :]) s = st.sum() if s > 0: stamps.append(st / s) psf = np.median(np.array(stamps), axis=0) psf /= psf.sum() # normalise so that a source of total flux F injected as F*psf measures # F_ap through the r=5 px aperture; we need the aperture correction yy, xx = np.mgrid[-HALF:HALF + 1, -HALF:HALF + 1] apfrac = psf[np.hypot(xx, yy) <= APR].sum() print(" aperture fraction inside r=%.0f px: %.4f" % (APR, apfrac)) return psf.astype(np.float32), apfrac def main(): t0 = time.time() hdr = fits.getheader(layout.path("master-Luminance.fit")) w = WCS(hdr) ny, nx = hdr["NAXIS2"], hdr["NAXIS1"] cx, cy = [float(v) for v in w.all_world2pix(NUC_RA, NUC_DEC, 0)] cat = dict(np.load(layout.path("_gc_cat.npz"), allow_pickle=True)) FW_MIN, FW_MAX = float(cat["FW_MIN"]), float(cat["FW_MAX"]) base_img = np.ascontiguousarray( fits.getdata(layout.path("master-Luminance.fit")).astype(np.float32)) psf, apfrac = build_psf(base_img, cat) np.save(layout.path("_gc_psf.npy"), psf) sep.set_extract_pixstack(1000000) K = gk() rec_mag, rec_rad, rec_ok = [], [], [] for run in range(NRUN): img = base_img.copy() # positions: uniform over the frame, min 40 px apart from each other xs = RNG.uniform(HALF + 5, nx - HALF - 5, NPER) ys = RNG.uniform(HALF + 5, ny - HALF - 5, NPER) mags = RNG.uniform(*MAGS, NPER) # total flux such that the r=5px aperture measures the intended mag ftot = 10 ** (-0.4 * (mags - ZP_L)) / apfrac for xi, yi, f in zip(xs, ys, ftot): ix, iy = int(round(xi)), int(round(yi)) img[iy - HALF:iy + HALF + 1, ix - HALF:ix + HALF + 1] += (f * psf) b = sep.Background(img, bw=BW, bh=BW, fw=FW, fh=FW) ds = img - b.back() rm = b.rms() o = sep.extract(ds, THRESH, err=rm, minarea=MINAREA, filter_kernel=K, filter_type="matched", deblend_nthresh=32, deblend_cont=0.005, clean=True) fl, fe, afl = sep.sum_circle(ds, o["x"], o["y"], APR, err=rm, subpix=5) fl2, _, _ = sep.sum_circle(ds, o["x"], o["y"], 2 * APR, err=rm, subpix=5) rh, _ = sep.flux_radius(ds, o["x"], o["y"], np.full(len(o), 6 * APR), 0.5, normflux=fl2, subpix=5) egain = float(hdr.get("EGAIN", 0.2467)) fet = np.sqrt(fe ** 2 + np.clip(fl, 0, None) / (egain * 12.0)) snr = fl / np.clip(fet, 1e-9, None) omag = -2.5 * np.log10(np.clip(fl, 1e-9, None)) + ZP_L ofwhm = 2.0 * rh oelong = o["a"] / np.clip(o["b"], 1e-6, None) ok_sel = ((snr >= 5.0) & ((o["flag"].astype(int) & 0b1110) == 0) & (ofwhm > FW_MIN) & (ofwhm < FW_MAX) & (oelong < 2.0) & (omag > 17.5) & (omag < 21.5) & (fl > 0)) ox, oy = o["x"][ok_sel], o["y"][ok_sel] # match each injected star to the surviving detections, within 2 px for xi, yi, mg in zip(xs, ys, mags): dd = np.hypot(ox - xi, oy - yi) hit = (dd.min() < 2.0) if len(dd) else False rec_mag.append(mg) rec_rad.append(np.hypot(xi - cx, yi - cy) * PIXSCALE / 60.0) rec_ok.append(bool(hit)) print(" run %d/%d done (%.0f s)" % (run + 1, NRUN, time.time() - t0), flush=True) del img, b, ds, rm, o rec_mag = np.array(rec_mag); rec_rad = np.array(rec_rad); rec_ok = np.array(rec_ok) np.savez(layout.path("_gc_complete.npz"), mag=rec_mag, rad=rec_rad, ok=rec_ok, apfrac=apfrac) print("\ninjected %d, recovered %d (%.1f%%)" % (len(rec_ok), rec_ok.sum(), 100 * rec_ok.mean())) print("\ncompleteness grid (rows = mag, cols = radius arcmin):") rb = [0, 2, 4, 6, 9, 13, 18, 30] mb = np.arange(17.5, 22.6, 0.5) print(" " + "".join("%7s" % ("%d-%d" % (rb[i], rb[i + 1])) for i in range(len(rb) - 1))) for j in range(len(mb) - 1): row = "%4.1f " % mb[j] for i in range(len(rb) - 1): s = ((rec_mag >= mb[j]) & (rec_mag < mb[j + 1]) & (rec_rad >= rb[i]) & (rec_rad < rb[i + 1])) row += "%7s" % ("%.2f" % rec_ok[s].mean() if s.sum() > 20 else "-") print(row) if __name__ == "__main__": main()