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.
102 lines
4 KiB
Python
102 lines
4 KiB
Python
"""How deep did the luminance master actually go?
|
|
|
|
Calibrates instrumental magnitudes against Gaia G, then reports the faintest
|
|
star still detected at 5 sigma. That number decides which follow-up analyses
|
|
are worth attempting on this data and which are wishful thinking.
|
|
"""
|
|
import os
|
|
|
|
import numpy as np
|
|
import sep
|
|
from astropy import units as u
|
|
from astropy.coordinates import SkyCoord
|
|
from astropy.io import fits
|
|
from astropy.wcs import WCS
|
|
|
|
import layout
|
|
|
|
OUT = layout.SESSION
|
|
CACHE = layout.path("_gaia_deep.npz")
|
|
|
|
with fits.open(layout.path("master-Luminance.fit")) as hd:
|
|
img = hd[0].data.astype(np.float32)
|
|
hdr = hd[0].header
|
|
wcs = WCS(hdr, naxis=2)
|
|
ny, nx = img.shape
|
|
|
|
bkg = sep.Background(img, bw=64, bh=64, fw=3, fh=3)
|
|
sub = img - bkg.back()
|
|
objs = sep.extract(sub, 3.0, err=bkg.globalrms, minarea=6, deblend_cont=0.005)
|
|
objs = objs[(objs["flag"] == 0) & (objs["npix"] > 6)]
|
|
flux, fluxerr, _ = sep.sum_circle(sub, objs["x"], objs["y"], 5.0,
|
|
err=bkg.globalrms, subpix=5)
|
|
snr = flux / np.maximum(fluxerr, 1e-9)
|
|
keep = (flux > 0) & (snr > 3)
|
|
objs, flux, snr = objs[keep], flux[keep], snr[keep]
|
|
print(f"{len(objs)} sources detected at SNR > 3 (rms {bkg.globalrms:.2f} ADU)")
|
|
|
|
sky = wcs.pixel_to_world(objs["x"], objs["y"])
|
|
centre = wcs.pixel_to_world(nx / 2, ny / 2)
|
|
|
|
if os.path.exists(CACHE):
|
|
z = np.load(CACHE)
|
|
gra, gdec, gmag = z["ra"], z["dec"], z["g"]
|
|
else:
|
|
from astroquery.gaia import Gaia
|
|
Gaia.ROW_LIMIT = 60000
|
|
job = Gaia.launch_job_async(f"""
|
|
SELECT ra, dec, phot_g_mean_mag FROM gaiadr3.gaia_source
|
|
WHERE 1 = CONTAINS(POINT('ICRS', ra, dec),
|
|
CIRCLE('ICRS', {centre.ra.deg}, {centre.dec.deg}, 0.42))
|
|
AND phot_g_mean_mag IS NOT NULL AND phot_g_mean_mag < 20.5
|
|
""")
|
|
t = job.get_results()
|
|
gra = np.asarray(t["ra"], float)
|
|
gdec = np.asarray(t["dec"], float)
|
|
gmag = np.asarray(t["phot_g_mean_mag"], float)
|
|
np.savez_compressed(CACHE, ra=gra, dec=gdec, g=gmag)
|
|
print(f"{len(gmag)} Gaia sources in the field, G down to {gmag.max():.2f}")
|
|
|
|
gcoord = SkyCoord(gra * u.deg, gdec * u.deg)
|
|
idx, sep2d, _ = sky.match_to_catalog_sky(gcoord)
|
|
matched = sep2d.arcsec < 1.5
|
|
print(f"{matched.sum()} detections matched to Gaia within 1.5\"")
|
|
|
|
inst = -2.5 * np.log10(flux[matched])
|
|
gm = gmag[idx[matched]]
|
|
# Fit the zero point on well exposed, unsaturated stars only.
|
|
fit = (gm > 12) & (gm < 17) & (snr[matched] > 20)
|
|
zp = float(np.median(gm[fit] - inst[fit]))
|
|
scatter = float(np.std(gm[fit] - inst[fit] - 0.0))
|
|
print(f"zero point {zp:.3f} (G = inst + zp) from {fit.sum()} stars, "
|
|
f"scatter {scatter:.3f} mag")
|
|
|
|
mag_all = -2.5 * np.log10(flux) + zp
|
|
# SNR falls monotonically with magnitude, so read the SNR=5 crossing off a
|
|
# running median rather than requiring sources to land in a narrow SNR bin.
|
|
order = np.argsort(mag_all)
|
|
ms, ss = mag_all[order], snr[order]
|
|
win = max(11, len(ms) // 60)
|
|
run_m = np.array([np.median(ms[i:i + win]) for i in range(0, len(ms) - win, win // 2)])
|
|
run_s = np.array([np.median(ss[i:i + win]) for i in range(0, len(ss) - win, win // 2)])
|
|
below = np.where(run_s < 5.0)[0]
|
|
lim5 = float(run_m[below[0]]) if len(below) else float(run_m[-1])
|
|
print(f"limiting magnitude at SNR 5: G ~ {lim5:.2f} "
|
|
f"(SNR range {snr.min():.1f}-{snr.max():.0f})")
|
|
print(f"faintest detection: G ~ {mag_all.max():.2f} (SNR "
|
|
f"{snr[np.argmax(mag_all)]:.1f})")
|
|
|
|
unmatched = ~matched
|
|
print(f"{unmatched.sum()} detections with NO Gaia counterpart "
|
|
f"({100 * unmatched.mean():.1f}% of sources)")
|
|
r_gal = np.hypot(objs["x"] - nx / 2, objs["y"] - ny / 2) * 0.5376 / 60.0
|
|
near = unmatched & (r_gal < 12.0) & (mag_all > 18.0) & (mag_all < 22.0)
|
|
print(f" of those, {near.sum()} lie within 12' of the galaxy at "
|
|
f"G 18-22: the magnitude and radius range of Centaurus A's "
|
|
f"globular cluster system")
|
|
|
|
# Surface brightness of the sky, a fair summary of how much the moon cost.
|
|
pixarea = 0.5376 ** 2
|
|
sky_adu = float(np.median(bkg.back()))
|
|
print(f"sky background {sky_adu:.1f} ADU/px -> "
|
|
f"{zp - 2.5 * np.log10(max(sky_adu, 1e-6) / pixarea):.2f} mag/arcsec^2")
|