"""Laser lines and photon pulses in archival ultraviolet data of nearby stars.""" import argparse import glob import json import math import os import re import time import urllib.parse import urllib.request from concurrent.futures import ThreadPoolExecutor import numpy as np from setitools import core from setitools.core import quiet_warnings # -------------------------------------------------------------------------------------------------------------- # archive pointings carry catalogue positions of any epoch since 1950: match proper-motion tracks, nearest star within the aperture def uv_nearby(plx): """SIMBAD stars above plx mas, planets excluded: main_id, ra, dec, pmra, pmdec, plx_value, sp_type, rvz_radvel, otype.""" return core.simbad("SELECT main_id, ra, dec, pmra, pmdec, plx_value, sp_type, rvz_radvel, otype FROM basic " "WHERE plx_value > %g AND otype NOT IN ('Pl', 'Pl?') AND ra IS NOT NULL" % plx) def uv_track_dist(star, ra, dec, mjd): """Arcsec from pointings (deg, MJD) to a star's proper-motion track from 1950 to the observation or 2000, whichever is later.""" ra0, de0 = float(star["ra"]), float(star["dec"]) x = ((np.asarray(ra, float) - ra0 + 180.0) % 360.0 - 180.0) * math.cos(math.radians(de0)) * 3600.0 y = (np.asarray(dec, float) - de0) * 3600.0 mx, my = float(star.get("pmra") or 0) / 1000.0, float(star.get("pmdec") or 0) / 1000.0 t1 = np.maximum((np.asarray(mjd, float) - 51544.5) / 365.25, 0.0) v2 = mx * mx + my * my t = np.clip((x * mx + y * my) / v2, -50.0, t1) if v2 > 0 else np.zeros_like(x) return np.hypot(x - mx * t, y - my * t) def uv_owner(dists, tol): """Nearest star to a pointing from {name: track distance} within tol arcsec, or None.""" if not dists: return None name = min(dists, key=dists.get) return name if dists[name] <= tol else None def uv_groups(lams, tols): """Indices of ascending wavelengths grouped where neighbours lie within the lower one's tolerance.""" if not len(lams): return [] out, cur = [], [0] for i in range(1, len(lams)): if lams[i] - lams[i - 1] < tols[i - 1]: cur.append(i) else: out.append(cur) cur = [i] return out + [cur] def uv_forced(items, lam_of, z_of, thr=5.0): """Stars whose forced measurement exceeds thr: lam_of(star) the wavelength in the star's frame, z_of(path, lam) the excess in noise units or None.""" hit = set() for star, path in items: if star in hit: continue z = z_of(path, lam_of(star)) if z is not None and z > thr: hit.add(star) return hit def uv_resolve(names): """{target name: SIMBAD main_id} for names SIMBAD knows, a spaced component letter ('HD 48915 B') also tried joined.""" var = {} for n in names: n = " ".join(str(n).split()) for v in (n, re.sub(r"(\d) ([A-Ca-c])$", r"\1\2", n)): var.setdefault(v, n) out = {} keys = sorted(var) for k in range(0, len(keys), 50): q = ("SELECT i.id, b.main_id FROM ident AS i JOIN basic AS b ON i.oidref = b.oid WHERE i.id IN (%s)" % ",".join("'%s'" % x.replace("'", "''") for x in keys[k:k + 50])) for r in core.simbad(q): out.setdefault(var.get(" ".join(str(r["id"]).split()), ""), " ".join(str(r["main_id"]).split())) out.pop("", None) return out # -------------------------------------------------------------------------------------------------------------- starcat_DOC = r"""Laser lines in StarCAT HST/STIS echelle spectra (Ayres 2010): unresolved 6-sigma peaks persistent in two exposures of a setting, absent from other stars in either frame. python -I seti.py ultraviolet starcat fetch python -I seti.py ultraviolet starcat search""" starcat_ROOT = os.environ.get("STARCAT_ROOT", core.data("starcat")) starcat_HERE = core.RESULTS starcat_BASE = "https://archive.stsci.edu/hlsps/starcat/fits/" starcat_RES = {"E140M": 45800.0, "E230M": 30000.0, "E140H": 114000.0, "E230H": 114000.0} starcat_AIRGLOW = (1215.67, 1302.17, 1304.86, 1306.03, 1355.60, 1358.51, 1199.55, 1200.22, 1200.71) # chromospheric and transition-region lines whose self-reversed cores leave narrow peaks: C III, Si III, N V, C II, Si IV, C IV, He II, Mg II starcat_CHROMOSPHERIC = (1175.71, 1206.50, 1238.82, 1242.80, 1334.53, 1335.71, 1393.76, 1402.77, 1548.20, 1550.77, 1640.4, 2795.53, 2796.35, 2802.70, 2803.53) starcat_Z_MIN = 6.0 starcat_C_KMS = 299792.458 starcat_ENVELOPE_A = 1.0 starcat_ENVELOPE_PCT = 90 def starcat_get(u): return core.get(u, timeout=180, retries=3, soft=True, pause=2) def starcat_fetch(): os.makedirs(os.path.join(starcat_ROOT, "fits"), exist_ok=True) top = starcat_get(starcat_BASE).decode("utf-8", "replace") objs = [o for o in re.findall(r'href="([^"?/][^"]*)/"', top)] jobs = [] def listing(o): t = starcat_get(starcat_BASE + o + "/") return [(o, f) for f in re.findall(r'href="(h_o[0-9a-z]+_spc\.fits)"', t.decode("utf-8", "replace"))] if t else [] with ThreadPoolExecutor(6) as ex: for js in ex.map(listing, objs): jobs += js print("%d objects, %d single exposures" % (len(objs), len(jobs)), flush=True) def one(j): o, f = j d = os.path.join(starcat_ROOT, "fits", o) p = os.path.join(d, f) if os.path.exists(p): return 1 b = starcat_get(starcat_BASE + o + "/" + f) if not b or not b.startswith(b"SIMPLE"): return 0 os.makedirs(d, exist_ok=True) open(p + ".tmp", "wb").write(b) os.replace(p + ".tmp", p) return 1 with ThreadPoolExecutor(6) as ex: n = sum(ex.map(one, jobs)) print("%d exposures on disk" % n) def starcat_load(path): """(wavelength, flux, error) of an exposure's good samples and its metadata; None for a setting outside RES.""" from astropy.io import fits with fits.open(path) as h: hd = h[0].header d = h[1].data[0] w, f, e, dq = (np.asarray(d[k], float) for k in ("WAVE", "FLUX", "ERROR", "DQ")) mode = str(hd.get("OPT_ELEM", "")).strip() if mode not in starcat_RES: return None meta = dict(obj=str(hd.get("CATNAME", "")).strip(), mode=mode, cenwave=hd.get("CENWAVE"), mjd=hd.get("EXPSTART"), ra=hd.get("RA_TARG"), dec=hd.get("DEC_TARG"), file=os.path.basename(path), exptime=hd.get("EXPTIME")) good = (e > 0) & np.isfinite(f) & ((dq.astype(int) & (4 | 8 | 16 | 32 | 64 | 128 | 256)) == 0) return w[good], f[good], e[good], meta def starcat_features(path): x = starcat_load(path) if x is None: return None w, f, e, meta = x if len(w) < 1000: return dict(meta, feats=[]) return dict(meta, feats=starcat_detect(w, f, e, meta["mode"])) def starcat_detect(w, f, e, mode): """Unresolved peaks: 6 sigma above the running upper envelope, 3 sigma above the brightest samples within 2.5 A, width 0.5-1.5 instrumental.""" from scipy.optimize import curve_fit step = float(np.median(np.diff(w))) # upper envelope: a running median between dense absorption lines makes the windows read as peaks from scipy.ndimage import percentile_filter win = max(int(starcat_ENVELOPE_A / step) | 1, 31) cont = percentile_filter(f, starcat_ENVELOPE_PCT, size=win, mode="nearest") z = (f - cont) / e # and above the brightest neighbours within 2.5 A, which line-dense windows reach from scipy.ndimage import maximum_filter1d half = max(int(2.5 / step), 30) left = maximum_filter1d(np.concatenate([np.full(half, -np.inf), f[:-half]]), size=half, origin=(half - 1) // 2) right = maximum_filter1d(np.concatenate([f[half:], np.full(half, -np.inf)]), size=half, origin=-(half // 2)) top = np.maximum(left, right) out = [] cand = np.where((z >= starcat_Z_MIN) & (z >= np.roll(z, 1)) & (z >= np.roll(z, -1)))[0] for i in cand: if i < 8 or i > len(w) - 9: continue lam = w[i] if any(abs(lam - a) < 0.6 for a in starcat_AIRGLOW) or any(abs(lam - a) < 1.2 for a in starcat_CHROMOSPHERIC): continue if not f[i] - top[i] > 3 * e[i]: continue inst = lam / starcat_RES[mode] sl = slice(i - 8, i + 9) x, y, s = w[sl], (f - cont)[sl], e[sl] try: p, _ = curve_fit(lambda x, a, m, sg: a * np.exp(-0.5 * ((x - m) / sg) ** 2), x, y, p0=[y[8], lam, inst / 2.355], sigma=s, maxfev=2000) except Exception: continue fwhm = 2.355 * abs(p[2]) if 0.5 * inst <= fwhm <= 1.5 * inst and abs(p[1] - lam) < inst: out.append(dict(lam=float(p[1]), z=float(z[i]), fwhm=float(fwhm), inst=float(inst), flux=float(p[0]))) return out def starcat_resolve(objs): """SIMBAD main_id, parallax, spectral type and RV of StarCAT objects {name: (ra, dec)}: by identifier, else by position among stars.""" def num(x): return float(x) if x not in (None, "") else None out = {} ident = {k: re.sub(r"^(HD|HR|HIP|GJ|GL|BD|CD|LHS|LTT|WD)(?=[0-9+-])", r"\1 ", k.replace("-", " "), flags=re.I) for k in objs} ids = sorted(set(ident.values())) for j in range(0, len(ids), 50): q = ("SELECT i.id, b.main_id, b.plx_value, b.sp_type, b.rvz_radvel FROM ident AS i JOIN basic AS b ON " "i.oidref = b.oid WHERE i.id IN (%s)" % ",".join("'%s'" % x.replace("'", "''") for x in ids[j:j + 50])) for r in core.simbad(q): key = " ".join(str(r["id"]).split()).lower() for k, v in ident.items(): if key == " ".join(v.split()).lower(): out[k] = dict(main_id=" ".join(r["main_id"].split()), plx=num(r.get("plx_value")), sp=r.get("sp_type") or "", rv=num(r.get("rvz_radvel"))) rest = {k: v for k, v in objs.items() if k not in out and v[0] is not None} m = core.xmatch([dict(source_id=k, ra=v[0], dec=v[1]) for k, v in rest.items()], "simbad", 10.0) if rest else [] best = {} for r in m: if r.get("otype") in ("Pl", "Pl?") or "Planet" in (r.get("main_type") or ""): continue if r["source_id"] not in best or float(r["angDist"]) < float(best[r["source_id"]]["angDist"]): best[r["source_id"]] = r for k, r in best.items(): out[k] = dict(main_id=" ".join(r["main_id"].split()), plx=num(r.get("plx")), sp=r.get("sp_type") or "", rv=num(r.get("radvel"))) return out def starcat_forced_z(w, f, e, lam): """Largest excess over the upper envelope in noise units within two samples of lam; None without coverage.""" from scipy.ndimage import percentile_filter j = int(np.searchsorted(w, lam)) if j < 3 or j > len(w) - 4: return None step = float(np.median(np.diff(w[max(j - 50, 0):j + 50]))) if abs(w[j] - lam) > 2 * step: return None lo, hi = max(j - 600, 0), min(j + 600, len(w)) cont = percentile_filter(f[lo:hi], starcat_ENVELOPE_PCT, size=max(int(starcat_ENVELOPE_A / step) | 1, 31), mode="nearest") k = j - lo return float(np.max((f[lo:hi][k - 2:k + 3] - cont[k - 2:k + 3]) / e[lo:hi][k - 2:k + 3])) def starcat__forced(job): """Forced measurements of candidate features in one exposure: (object, [(z at rest wavelength, z at observed wavelength in the same setting)]).""" quiet_warnings() path, rv, feats = job x = starcat_load(path) if x is None: return None w, f, e, meta = x out = [] for p in feats: zr = starcat_forced_z(w, f, e, p["rest"] * (1 + rv / starcat_C_KMS)) same = meta["mode"] == p["mode"] and meta["cenwave"] in p["cenwaves"] out.append((zr, starcat_forced_z(w, f, e, p["lam"]) if same else None)) return meta["obj"], out def starcat_search(): from concurrent.futures import ProcessPoolExecutor paths = [os.path.join(dp, f) for dp, _, fs in os.walk(os.path.join(starcat_ROOT, "fits")) for f in fs if f.endswith("_spc.fits")] procs = int(os.environ.get("UV_PROCS", "2")) with ProcessPoolExecutor(procs) as ex: res = [r for r in ex.map(starcat_features, paths, chunksize=8) if r] print("%d exposures of %d objects; %d narrow peaks at %d sigma or more" % ( len(res), len({r["obj"] for r in res}), sum(len(r["feats"]) for r in res), starcat_Z_MIN), flush=True) # persistent features per object and setting by = {} for r in res: by.setdefault((r["obj"], r["mode"]), []).append(r) pers = [] for (obj, mode), rs in by.items(): if len(rs) < 2: continue allf = sorted(((x, r["file"], r["cenwave"]) for r in rs for x in r["feats"]), key=lambda t: t[0]["lam"]) for grp in uv_groups([x["lam"] for x, _, _ in allf], [x["inst"] / 3 for x, _, _ in allf]): files = {allf[j][1] for j in grp} if len(files) >= 2: xs = [allf[j][0] for j in grp] pers.append(dict(obj=obj, mode=mode, lam=float(np.mean([y["lam"] for y in xs])), n=len(files), of=len(rs), z=float(max(y["z"] for y in xs)), inst=xs[0]["inst"], cenwaves=sorted({allf[j][2] for j in grp}), ra=rs[0]["ra"], dec=rs[0]["dec"])) print("%d persistent features (two or more exposures in one setting) in %d objects" % ( len(pers), len({p["obj"] for p in pers})), flush=True) names = starcat_resolve({r["obj"]: (r["ra"], r["dec"]) for r in res}) rv = {k: v["rv"] or 0.0 for k, v in names.items()} for p in pers: p["rest"] = p["lam"] / (1 + rv.get(p["obj"], 0.0) / starcat_C_KMS) p["main_id"], p["sp"] = names.get(p["obj"], {}).get("main_id", p["obj"]), names.get(p["obj"], {}).get("sp", "") # other stars' peaks: rest frame (stellar), observed frame (instrument or sky) pk = [(x["lam"], r["obj"]) for r in res for x in r["feats"]] obsw = np.array([lam for lam, _ in pk]) rest = np.array([lam / (1 + rv.get(o, 0.0) / starcat_C_KMS) for lam, o in pk]) objs = np.array([o for _, o in pk]) for p in pers: near_r = np.abs(rest - p["rest"]) / p["rest"] * starcat_C_KMS < 15.0 near_o = np.abs(obsw - p["lam"]) / p["lam"] * starcat_C_KMS < 15.0 p["stars_rest"] = len(set(objs[near_r]) - {p["obj"]}) p["stars_obs"] = len(set(objs[near_o]) - {p["obj"]}) left = [p for p in pers if p["stars_rest"] < 2 and p["stars_obs"] < 2] # forced measurements in other stars below detection: rest wavelength (stellar), observed wavelength in the same setting (instrument) obj_of = {r["file"]: r["obj"] for r in res} jobs = [(q, rv.get(obj_of.get(os.path.basename(q)), 0.0), left) for q in paths if os.path.basename(q) in obj_of] hits_r, hits_o = [set() for _ in left], [set() for _ in left] with ProcessPoolExecutor(procs) as ex: for out in ex.map(starcat__forced, jobs, chunksize=8): if not out: continue obj, zs = out for k, (zr, zo) in enumerate(zs): if obj != left[k]["obj"]: if zr is not None and zr > 5.0: hits_r[k].add(obj) if zo is not None and zo > 5.0: hits_o[k].add(obj) for k, p in enumerate(left): p["forced_rest"], p["forced_obs"] = sorted(hits_r[k]), sorted(hits_o[k]) unique = sorted([p for p in left if len(p["forced_rest"]) < 2 and len(p["forced_obs"]) < 2], key=lambda p: -p["z"]) print("%d persistent features in at most one other star in either frame; %d of them above 5 sigma in two or more " "other stars when measured there; %d left:" % (len(left), len(left) - len(unique), len(unique))) for p in unique[:60]: print(" %-18s %s %.3f A (rest %.3f), in %d of %d exposures, peak %.1f sigma, other stars %d/%d, measured %d/%d" % (p["obj"], p["mode"], p["lam"], p["rest"], p["n"], p["of"], p["z"], p["stars_rest"], p["stars_obs"], len(p["forced_rest"]), len(p["forced_obs"]))) json.dump(dict(exposures=len(res), objects=len({r["obj"] for r in res}), peaks=sum(len(r["feats"]) for r in res), persistent=pers), open(os.path.join(starcat_HERE, "uv_laser.json"), "w"), indent=1) def starcat_main(): p = argparse.ArgumentParser() p.add_argument("mode", choices=("fetch", "search")) a = p.parse_args() {"fetch": starcat_fetch, "search": starcat_search}[a.mode]() def starcat_cli(): quiet_warnings() starcat_main() # -------------------------------------------------------------------------------------------------------------- cos_DOC = r"""Laser lines in HST/COS FUV spectra of stars within 10 pc: unresolved peaks at p 1e-9, persistent across FP-POS, absent from other stars and their summed spectra. python -I seti.py ultraviolet cos spectra (STIS_ROOT event files; COS_PROCS processes) python -I seti.py ultraviolet cos search python -I seti.py ultraviolet cos identify (NIST ASD lines, cached beside STIS_ROOT)""" cos_ROOT = os.environ.get("STIS_ROOT", core.data("stis_tag")) cos_HERE = core.RESULTS cos_C_KMS = 299792.458 cos_BIN_PX = 2 cos_LSF_PX = 6.5 cos_AIRGLOW = ((1215.67, 2.0), (1302.17, 0.6), (1304.86, 0.6), (1306.03, 0.6), (1355.60, 0.6), (1358.51, 0.6), (1199.55, 0.6), (1200.22, 0.6), (1200.71, 0.6)) cos_STRONG = (1175.71, 1206.50, 1238.82, 1242.80, 1334.53, 1335.71, 1393.76, 1402.77, 1548.20, 1550.77, 1640.4) cos_P_MAX = 1e-9 cos_DQ_BAD = 2 | 8 | 16 | 32 | 64 | 128 | 512 | 1024 | 2048 | 4096 | 8192 | 16384 cos_APERTURE_TOL = 3.0 # arcsec; the primary science aperture is 2.5 arcsec across def cos_owned(rs): """Spectra whose pointing holds a star within 10 pc, each given that star as 'star'.""" inv = json.load(open(os.path.join(cos_ROOT, "inventory.json"))) tags = {t["file"]: t for t in inv["tags"]} out = [] for r in rs: tg = tags.get(r["file"]) o = inv["obs"].get(tg["parent"], inv["obs"].get(tg["obsid"], {})) if tg else {} star = uv_owner(o.get("dist", {}), cos_APERTURE_TOL) if star: out.append(dict(r, star=star)) return out cos_CAL_TARGETS = ("DARK", "BIAS", "DEUTERIUM", "WAVE", "WAVEHITM", "WAVELINE", "FLAT", "CCDFLAT", "NONE") def cos_left(p): """True for a persistent feature no two other stars show in either frame and absent from their summed spectra.""" return p["stars_rest"] < 2 and p["stars_obs"] < 2 and p.get("template_p", 1.0) > 1e-5 def cos_science(hd): """False for calibration exposures (darks, flats, lamps) filed under a science visit.""" et = str(hd.get("EXPTYPE", "EXTERNAL/SCI")).strip().upper() tn = str(hd.get("TARGNAME", "")).strip().upper() return et in ("EXTERNAL/SCI", "") and not any(tn.startswith(c) for c in cos_CAL_TARGETS) def cos_spectrum(path): """Heliocentric count spectrum of one COS FUV corrtag segment: on-band and scaled background counts per two-pixel bin, with metadata.""" from astropy.io import fits with fits.open(path) as h: hd = h[0].header if str(hd.get("DETECTOR", "")).strip() != "FUV": return None if not cos_science(hd): return None seg = str(hd.get("SEGMENT", "")).strip()[-1:] e = h["EVENTS"] eh, d = e.header, e.data if d is None or len(d) < 100: return None cols = d.columns.names t = np.asarray(d["TIME"], float) xf = np.asarray(d["XFULL"], float) yf = np.asarray(d["YFULL"], float) wl = np.asarray(d["WAVELENGTH"], float) good = (np.asarray(d["DQ"], int) & cos_DQ_BAD) == 0 if "DQ" in cols else np.ones(len(t), bool) if "PHA" in cols: pha = np.asarray(d["PHA"], int) good &= (pha >= int(eh.get("PHALOWR" + seg, 2))) & (pha <= int(eh.get("PHAUPPR" + seg, 23))) good &= wl > 0 for i in range(1, 40): if "LMP_ON%d" % i in eh: on = float(eh["LMP_ON%d" % i]) off = float(eh.get("LMPOFF%d" % i, on + 15.0)) good &= ~((t >= on - 2.0) & (t <= off + 2.0)) loc, hgt = eh.get("SP_LOC_" + seg), eh.get("SP_HGT_" + seg) if loc is None or hgt is None: return None slp = float(eh.get("SP_SLP_" + seg, 0.0) or 0.0) best = None for xr in (0.0, 8192.0): inb = np.abs(yf - (float(loc) + slp * (xf - xr))) <= float(hgt) / 2 if best is None or (good & inb).sum() > (good & best).sum(): best = inb on_band = good & best bk = np.zeros(len(t), bool) bh = 0.0 for k in (1, 2): b, bhh = eh.get("B_BKG%d_%s" % (k, seg)), eh.get("B_HGT%d_%s" % (k, seg)) if b is not None and bhh: bk |= np.abs(yf - float(b)) <= float(bhh) / 2 bh += float(bhh) bk &= good vh = float(eh.get("V_HELIO", 0.0) or 0.0) meta = dict(file=os.path.basename(path), grating=str(hd.get("OPT_ELEM", "")).strip(), cenwave=hd.get("CENWAVE"), fppos=hd.get("FPPOS"), segment=seg, target=str(hd.get("TARGNAME", "")).strip(), mjd=float(eh.get("EXPSTART", 0) or 0), exptime=float(eh.get("EXPTIME", 0) or 0), v_helio=vh) if on_band.sum() < 200: return None edges = np.arange(0, 16384 + cos_BIN_PX, cos_BIN_PX) c_on = np.histogram(xf[on_band], edges)[0].astype(float) c_bk = np.histogram(xf[bk], edges)[0].astype(float) * (float(hgt) / bh if bh else 0.0) # dispersion: cubic wavelength against XFULL per segment from the band's events sel = on_band if sel.sum() > 50: cf = np.polyfit(xf[sel], wl[sel], 3) else: return None centres = 0.5 * (edges[:-1] + edges[1:]) lam = np.polyval(cf, centres) lo, hi = np.percentile(xf[on_band], [0.5, 99.5]) keep = (centres >= lo) & (centres <= hi) return dict(meta, lam=lam[keep].tolist(), on=c_on[keep].tolist(), bk=c_bk[keep].tolist(), disp=float(abs(np.polyval(np.polyder(cf), 8192.0)))) def cos_spectra(): from concurrent.futures import ProcessPoolExecutor inv = json.load(open(os.path.join(cos_ROOT, "inventory.json"))) obs = inv["obs"] files = [(os.path.join(cos_ROOT, "events", tg["file"]), obs.get(tg["parent"], obs.get(tg["obsid"], {})).get("stars", [])) for tg in inv["tags"] if tg["file"].endswith(("_corrtag_a.fits", "_corrtag_b.fits"))] files = [(p, s) for p, s in files if os.path.exists(p)] out = [] with ProcessPoolExecutor(int(os.environ.get("COS_PROCS", "4"))) as ex: for (p, s), r in zip(files, ex.map(cos__safe_spectrum, [p for p, _ in files], chunksize=4)): if r: r["stars"] = [" ".join(x.split()) for x in s] out.append(r) json.dump(out, open(os.path.join(cos_ROOT, "cos_spectra.json"), "w")) print("%d COS FUV segment spectra of %d targets" % (len(out), len({r["target"] for r in out}))) def cos__safe_spectrum(p): try: return cos_spectrum(p) except Exception as e: print(os.path.basename(p), "failed:", str(e)[:80], flush=True) return None def cos_peaks(r, sign, p_max=cos_P_MAX): """Narrow peaks of one spectrum: (heliocentric wavelength, significance, width in pixels, net counts).""" from scipy.ndimage import percentile_filter from scipy.optimize import curve_fit from scipy.stats import poisson lam = np.array(r["lam"]) * (1 + sign * r["v_helio"] / cos_C_KMS) on, bk = np.array(r["on"]), np.array(r["bk"]) if len(lam) < 200: return [] s3 = np.convolve(on, np.ones(3), "same") b3 = np.convolve(bk, np.ones(3), "same") step = abs(float(np.median(np.diff(lam)))) win = max(int(1.0 / step) * 2 + 1, 31) env = percentile_filter(s3, 90, size=win, mode="nearest") mu = np.maximum(env, b3 + 0.5) p = poisson.sf(s3 - 1, mu) cand = np.where((p < p_max) & (s3 >= np.roll(s3, 1)) & (s3 >= np.roll(s3, -1)))[0] out = [] net = on - bk for i in cand: if i < 8 or i > len(lam) - 9: continue l0 = lam[i] if any(abs(l0 - a) < w for a, w in cos_AIRGLOW) or any(abs(l0 - a) < 1.2 for a in cos_STRONG): continue x = np.arange(-8, 9) * cos_BIN_PX y = net[i - 8:i + 9] - np.median(net[max(i - 40, 0):i + 40]) try: q, _ = curve_fit(lambda x, a, m, sg: a * np.exp(-0.5 * ((x - m) / sg) ** 2), x, y, p0=[max(y[8], 1.0), 0.0, cos_LSF_PX / 2.355], maxfev=2000) except Exception: continue fw = 2.355 * abs(q[2]) if 0.5 * cos_LSF_PX <= fw <= 1.5 * cos_LSF_PX and abs(q[1]) < cos_LSF_PX: out.append(dict(lam=float(np.interp(i + q[1] / cos_BIN_PX, np.arange(len(lam)), lam)), p=float(max(p[i], 1e-300)), fwhm_px=float(fw), net=float(q[0]))) return out def cos_search(): rs = cos_owned(json.load(open(os.path.join(cos_ROOT, "cos_spectra.json")))) # CalCOS HELCORR already makes WAVELENGTH heliocentric best_sign = 0 feats, loose = [], [] for r in rs: for f in cos_peaks(r, best_sign, 1e-5): g = dict(f, file=r["file"], grating=r["grating"], fppos=r["fppos"], target=r["target"], disp=r["disp"], star=r["star"]) loose.append(g) if f["p"] < cos_P_MAX: feats.append(g) print("%d spectra of %d stars; %d narrow peaks at p < %.0e" % (len(rs), len({r["star"] for r in rs}), len(feats), cos_P_MAX), flush=True) by = {} for f in feats: by.setdefault((f["star"], f["grating"]), []).append(f) pers = [] for (star, grating), fs in by.items(): fs.sort(key=lambda f: f["lam"]) used = set() for k, f in enumerate(fs): if k in used: continue tol = cos_BIN_PX * f["disp"] * 3 / 3.0 grp = [j for j, g in enumerate(fs) if abs(g["lam"] - f["lam"]) < tol] files = {fs[j]["file"][:9] for j in grp} # persistence across FP-POS: a fixed wavelength moves on the detector if len(files) >= 2 and len({fs[j]["fppos"] for j in grp}) >= 2: used.update(grp) pers.append(dict(star=star, grating=grating, lam=float(np.mean([fs[j]["lam"] for j in grp])), n=len(files), p=min(fs[j]["p"] for j in grp), fppos=sorted({fs[j]["fppos"] for j in grp}), disp=f["disp"])) print("%d persistent features (two or more exposures) in %d stars" % (len(pers), len({p["star"] for p in pers})), flush=True) names = sorted({p["star"] for p in pers} | {g["star"] for g in loose}) rv = {} for n in names: try: q = core.simbad("SELECT b.rvz_radvel FROM basic b JOIN ident i ON i.oidref = b.oid WHERE i.id = '%s'" % n.replace("'", "''")) rv[n] = float(q[0]["rvz_radvel"]) if q and q[0]["rvz_radvel"] not in ("", None) else 0.0 except Exception: rv[n] = 0.0 for p in pers: p["rest"] = p["lam"] / (1 + rv.get(p["star"], 0.0) / cos_C_KMS) # other stars' peaks at p < 1e-5 within 15 km/s: rest frame (stellar), observed frame (instrument or sky) rest = np.array([g["lam"] / (1 + rv.get(g["star"], 0.0) / cos_C_KMS) for g in loose]) obsw = np.array([g["lam"] for g in loose]) star = np.array([g["star"] for g in loose]) for p in pers: near_r = np.abs(rest - p["rest"]) / p["rest"] * cos_C_KMS < 15.0 near_o = np.abs(obsw - p["lam"]) / p["lam"] * cos_C_KMS < 15.0 p["stars_rest"] = len(set(star[near_r]) - {p["star"]}) p["stars_obs"] = len(set(star[near_o]) - {p["star"]}) # deep test: other stars summed in their rest frames, a peak at Poisson p < 1e-5 per resolution element over the running 90th percentile from scipy.stats import poisson grid = np.arange(1100.0, 1800.0, 0.01) stack = {} for r in rs: st = r["star"] lam = np.array(r["lam"]) * (1 + best_sign * r["v_helio"] / cos_C_KMS) / (1 + rv.get(st, 0.0) / cos_C_KMS) on = np.array(r["on"]) idx = np.clip(np.searchsorted(grid, lam), 0, len(grid) - 1) h = stack.setdefault((r["grating"], st), np.zeros(len(grid))) np.add.at(h, idx, on) for p in pers: tot = np.zeros(len(grid)) for (g, st), h in stack.items(): if g == p["grating"] and st != p["star"]: tot += h res_el = 3 * cos_BIN_PX * p["disp"] k = max(int(round(res_el / 0.01)), 1) s_k = np.convolve(tot, np.ones(k), "same") j = int(np.searchsorted(grid, p["rest"])) lo, hi = max(j - 100, 0), min(j + 101, len(grid)) env = np.percentile(s_k[lo:hi], 90) i0 = int(np.argmax(s_k[max(j - k, 0):j + k + 1])) + max(j - k, 0) p["template_counts"] = float(s_k[i0]) p["template_p"] = float(poisson.sf(s_k[i0] - 1, max(env, 0.5))) unique = sorted([p for p in pers if cos_left(p)], key=lambda p: p["p"]) print("%d persistent features also in the other stars' summed spectra (template p < 1e-5)" % sum(1 for p in pers if p["template_p"] <= 1e-5)) from collections import Counter print("persistent features by star: " + ", ".join("%s %d" % kv for kv in Counter(p["star"] for p in pers).most_common())) print("%d persistent features in at most one other star in either frame and absent from the summed spectra:" % len(unique)) for p in unique[:80]: print(" %-22s %s %.3f A (rest %.3f), in %d exposures (FP-POS %s), p %.1e, other stars %d/%d, other stars' " "sum %.0f counts (p %.1e)" % (p["star"][:22], p["grating"], p["lam"], p["rest"], p["n"], ",".join(str(x) for x in p["fppos"]), p["p"], p["stars_rest"], p["stars_obs"], p["template_counts"], p["template_p"])) json.dump(dict(sign=best_sign, spectra=len(rs), stars=len({r["star"] for r in rs}), peaks=len(feats), persistent=pers), open(os.path.join(cos_HERE, "cos_laser.json"), "w"), indent=1) cos_IONS = ("H I", "He II", "C I", "C II", "C III", "C IV", "N I", "N II", "N III", "N IV", "N V", "O I", "O II", "O III", "O IV", "O V", "Si I", "Si II", "Si III", "Si IV", "S I", "S II", "S III", "S IV", "Fe I", "Fe II", "Fe III", "Fe XII", "Fe XXI", "Mg I", "Mg II", "Al I", "Al II", "Al III", "Ca II", "Cl I", "Cl II", "P I", "P II", "Ni II", "Cr II", "Mn II", "Co II", "Zn II", "Ti II", "Ne V") def cos_identify(): """Persistent features within 15 km/s of NIST chromospheric and transition-region lines (intensity listed or A above 1e7 /s, low levels for neutral and singly ionised).""" lines = [] for ion in cos_IONS: lines += core.nist_lines(ion, 1100, 1800, cache_dir=os.path.join(os.path.dirname(cos_ROOT), "nist")) def num(x): try: return float(str(x).strip("[]()+x? ")) except ValueError: return float("nan") strong = [] for ln in lines: stage = ln["ion"].split()[1] ei, aki = num(ln["ei"]), num(ln["aki"]) if (ln["intens"].strip() or aki > 1e7) and (stage not in ("I", "II") or ei < 30000): strong.append(ln) strong.sort(key=lambda ln: ln["lam"]) lam = np.array([ln["lam"] for ln in strong]) d = json.load(open(os.path.join(cos_HERE, "cos_laser.json"))) left = [p for p in d["persistent"] if cos_left(p)] rng = np.random.default_rng(1) r = rng.uniform(1150, 1750, 20000) k = np.searchsorted(lam, r).clip(1, len(lam) - 1) dv = np.minimum(np.abs(lam[k] - r), np.abs(lam[k - 1] - r)) / r * cos_C_KMS print("%d NIST lines, %d strong or low-excitation; a random wavelength lies within 2 and 15 km/s of one with " "probability %.2f and %.2f" % (len(lines), len(strong), np.mean(dv < 2), np.mean(dv < 15))) unid = [] for p in sorted(left, key=lambda p: p["p"]): dv = (lam - p["rest"]) / p["rest"] * cos_C_KMS j = int(np.argmin(np.abs(dv))) p["id"] = "%s %.3f (%+.0f km/s)" % (strong[j]["ion"], strong[j]["lam"], dv[j]) if abs(dv[j]) < 15 else None if not p["id"]: unid.append(p) print(" %-14s %s rest %.3f, %d exposures, p %.0e: %s" % (p["star"][:14], p["grating"], p["rest"], p["n"], p["p"], p["id"] or "no line within 15 km/s")) print("%d of %d persistent features unique to one star have no line within 15 km/s" % (len(unid), len(left))) json.dump(dict(d, identified=left), open(os.path.join(cos_HERE, "cos_laser.json"), "w"), indent=1) def cos_main(): p = argparse.ArgumentParser() p.add_argument("mode", choices=("spectra", "search", "identify")) a = p.parse_args() {"spectra": cos_spectra, "search": cos_search, "identify": cos_identify}[a.mode]() def cos_cli(): quiet_warnings() cos_main() # -------------------------------------------------------------------------------------------------------------- iue_DOC = r"""Laser lines in IUE high-dispersion spectra of stars within 10 pc: unresolved 6-sigma peaks persistent in two images of one camera and aperture, absent from other stars and the other aperture. python -I seti.py ultraviolet iue inventory python -I seti.py ultraviolet iue fetch [--threads 4] python -I seti.py ultraviolet iue search""" iue_ROOT = os.environ.get("IUE_ROOT", core.data("iue")) iue_HERE = core.RESULTS iue_DL = "https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:IUE/url/missions/iue/data/%s/%05d/%s.mxhi.gz" iue_R = 10000.0 iue_Z_MIN = float(os.environ.get("IUE_ZMIN", "6.0")) iue_N_MIN = int(os.environ.get("IUE_NMIN", "2")) iue_C_KMS = 299792.458 iue_MATCH_ARCSEC = 30.0 iue_APERTURE_TOL = 10.0 # arcsec; the large aperture is 10 by 20 arcsec def iue_owned(inv): """One inventory entry per image: the nearest star within the aperture tolerance.""" by = {} for o in inv: by.setdefault(o["obs_id"], []).append(o) out = [] for v in by.values(): # archived coordinates err by up to 30 arcsec: the target name's star first named = [o for o in v if o.get("named")] o = min(named or v, key=lambda o: o["sep"]) if named or o["sep"] <= iue_APERTURE_TOL: out.append(o) return out def iue_inventory(): os.makedirs(iue_ROOT, exist_ok=True) stars = uv_nearby(100) every, page = [], 1 while True: d = hst_pulses_invoke({"service": "Mast.Caom.Filtered", "format": "json", "pagesize": 50000, "page": page, "params": {"columns": "obsid,obs_id,instrument_name,target_name,t_min,s_ra,s_dec", "filters": [{"paramName": "obs_collection", "values": ["IUE"]}, {"paramName": "filters", "values": ["HIGH DISP"]}]}}) rows = d.get("data", []) every += rows if len(rows) < 50000: break page += 1 print("%d IUE high-dispersion observations" % len(every), flush=True) ora = np.array([float(x["s_ra"]) if x.get("s_ra") is not None else np.nan for x in every]) ode = np.array([float(x["s_dec"]) if x.get("s_dec") is not None else np.nan for x in every]) tmin = np.array([float(x["t_min"]) if x.get("t_min") is not None else 46000.0 for x in every]) out = [] for s in stars: sep = uv_track_dist(s, ora, ode, tmin) for j in np.where(sep < iue_MATCH_ARCSEC)[0]: x = every[j] out.append(dict(star=" ".join(s["main_id"].split()), rv=s.get("rvz_radvel"), obs_id=x["obs_id"], camera=x["instrument_name"], target=x["target_name"], mjd=float(tmin[j]), sep=float(sep[j]))) names = uv_resolve({o["target"] for o in out}) for o in out: o["named"] = names.get(" ".join(str(o["target"]).split())) == o["star"] json.dump(out, open(os.path.join(iue_ROOT, "inventory.json"), "w"), indent=1) print("%d high-dispersion images of %d objects within 10 pc" % (len(out), len({o["star"] for o in out}))) def iue_url(obs_id): cam = "".join(c for c in obs_id if c.isalpha()).lower() num = int("".join(c for c in obs_id if c.isdigit())) return iue_DL % (cam, num // 1000 * 1000, obs_id.lower()) def iue_fetch(threads): from concurrent.futures import ThreadPoolExecutor inv = json.load(open(os.path.join(iue_ROOT, "inventory.json"))) os.makedirs(os.path.join(iue_ROOT, "mxhi"), exist_ok=True) def one(o): p = os.path.join(iue_ROOT, "mxhi", o + ".mxhi.gz") if os.path.exists(p): return 1 b = core.get(iue_url(o), timeout=180, retries=3, soft=True, pause=2) if not b or b[:2] != b"\x1f\x8b" or len(b) < 5000: return 0 open(p + ".tmp", "wb").write(b) os.replace(p + ".tmp", p) return 1 ids = sorted({o["obs_id"] for o in inv}) with ThreadPoolExecutor(threads) as ex: n = sum(ex.map(one, ids)) print("%d of %d MXHI files on disk" % (n, len(ids))) def iue_orders(path): """[(wavelength, net, noise, abs_cal)] of an image's orders at QUALITY 0; ABS_CAL converts flux numbers to erg/cm2/s/A.""" from astropy.io import fits out = [] with fits.open(path) as h: d = h[1].data names = [c.upper() for c in d.columns.names] for row in d: n = int(row["NPOINTS"]) if n < 100: continue w = float(row["WAVELENGTH"]) + float(row["DELTAW"]) * np.arange(n) f = np.asarray(row["NET"][:n], float) e = np.asarray(row["NOISE"][:n], float) q = np.asarray(row["QUALITY"][:n], int) cal = np.asarray(row["ABS_CAL"][:n], float) if "ABS_CAL" in names else np.full(n, np.nan) ok = (q == 0) & np.isfinite(f) & np.isfinite(e) & (e > 0) if ok.sum() < 80: continue out.append((w[ok], f[ok], e[ok], cal[ok])) return out def iue_detect(w, f, e): """Unresolved peaks of one order: 6 sigma above a 5 A running upper envelope and 3 sigma above the brightest samples within 2.5 A.""" from scipy.ndimage import maximum_filter1d, percentile_filter from scipy.optimize import curve_fit out = [] step = float(np.median(np.diff(w))) env = percentile_filter(f, 90, size=max(int(5.0 / step) | 1, 61), mode="nearest") z = (f - env) / e half = max(int(2.5 / step), 10) left = maximum_filter1d(np.concatenate([np.full(half, -np.inf), f[:-half]]), size=half, origin=(half - 1) // 2) right = maximum_filter1d(np.concatenate([f[half:], np.full(half, -np.inf)]), size=half, origin=-(half // 2)) top = np.maximum(left, right) for i in np.where((z >= iue_Z_MIN) & (z >= np.roll(z, 1)) & (z >= np.roll(z, -1)))[0]: if i < 6 or i > len(w) - 7 or not f[i] - top[i] > 3 * e[i]: continue lam = w[i] if any(abs(lam - a) < 0.6 for a in starcat_AIRGLOW) or any(abs(lam - a) < 1.2 for a in starcat_CHROMOSPHERIC): continue inst = lam / iue_R x, y = w[i - 6:i + 7], (f - env)[i - 6:i + 7] try: p, _ = curve_fit(lambda x, a, m, sg: a * np.exp(-0.5 * ((x - m) / sg) ** 2), x, y, p0=[y[6], lam, inst / 2.355], sigma=e[i - 6:i + 7], maxfev=2000) except Exception: continue fw = 2.355 * abs(p[2]) if 0.5 * inst <= fw <= 1.5 * inst and abs(p[1] - lam) < inst: out.append(dict(lam=float(p[1]), z=float(z[i]), inst=float(inst))) return out def iue_features(path): return [x for w, f, e, _ in iue_orders(path) for x in iue_detect(w, f, e)] def iue_aperture(path): """LARGE or SMALL: the aperture whose spectrum an MXHI file holds.""" from astropy.io import fits with fits.open(path) as h: return str(h[0].header.get("APERTURE", "")).strip().upper() or "?" def iue_forced_z(path, lam): """Largest excess over the upper envelope in noise units within two pixels of lam; None without coverage.""" from scipy.ndimage import percentile_filter best = None for w, f, e, _ in iue_orders(path): if w[0] + 1 < lam < w[-1] - 1: j = int(np.searchsorted(w, lam)) step = float(np.median(np.diff(w))) env = percentile_filter(f, 90, size=max(int(5.0 / step) | 1, 61), mode="nearest") v = float(np.max((f[j - 2:j + 3] - env[j - 2:j + 3]) / e[j - 2:j + 3])) best = v if best is None else max(best, v) return best def iue__feat(job): quiet_warnings() p = os.path.join(iue_ROOT, "mxhi", job["obs_id"] + ".mxhi.gz") if not os.path.exists(p): return job, None try: return job, (iue_features(p), iue_aperture(p)) except Exception: return job, None def iue_search(): from concurrent.futures import ProcessPoolExecutor inv = iue_owned(json.load(open(os.path.join(iue_ROOT, "inventory.json")))) feats, n_read, apert, rvs = [], 0, {}, {} with ProcessPoolExecutor(int(os.environ.get("IUE_PROCS", "4"))) as ex: for job, res in ex.map(iue__feat, inv, chunksize=8): if res is None: continue fs, ap = res n_read += 1 apert[(job["star"], job["obs_id"])] = (job["camera"], ap) rv = float(job["rv"]) if job.get("rv") not in (None, "") else 0.0 rvs[job["star"]] = rv for f in fs: feats.append(dict(f, star=job["star"], camera=job["camera"], aperture=ap, obs_id=job["obs_id"], rest=f["lam"] / (1 + rv / iue_C_KMS))) print("%d images of %d objects read; %d narrow peaks at %.0f sigma" % (n_read, len({j["star"] for j in inv}), len(feats), iue_Z_MIN), flush=True) # persistence within one camera and aperture, each with its own detector position by = {} for f in feats: by.setdefault((f["star"], f["camera"], f["aperture"]), []).append(f) pers = [] for (star, cam, ap), fs in by.items(): fs.sort(key=lambda f: f["lam"]) for grp in uv_groups([f["lam"] for f in fs], [f["inst"] / 3 for f in fs]): imgs = {fs[j]["obs_id"] for j in grp} if len(imgs) >= iue_N_MIN: pers.append(dict(star=star, camera=cam, aperture=ap, lam=float(np.mean([fs[j]["lam"] for j in grp])), rest=float(np.mean([fs[j]["rest"] for j in grp])), n=len(imgs), of=sum(1 for (s, o), (c, a) in apert.items() if s == star and c == cam and a == ap), z=max(fs[j]["z"] for j in grp), inst=fs[grp[0]]["inst"])) print("%d persistent features in %d stars" % (len(pers), len({p["star"] for p in pers})), flush=True) rest = np.array([f["rest"] for f in feats]) obsw = np.array([f["lam"] for f in feats]) star = np.array([f["star"] for f in feats]) for p in pers: near_r = np.abs(rest - p["rest"]) / p["rest"] * iue_C_KMS < 15.0 near_o = np.abs(obsw - p["lam"]) / p["lam"] * iue_C_KMS < 15.0 p["stars_rest"] = len(set(star[near_r]) - {p["star"]}) p["stars_obs"] = len(set(star[near_o]) - {p["star"]}) left = sorted([p for p in pers if p["stars_rest"] < 2 and p["stars_obs"] < 2], key=lambda p: -p["z"]) # a detector-position feature is absent through the other aperture for p in left: other = [o for (s, o), (c, a) in apert.items() if s == p["star"] and c == p["camera"] and a != p["aperture"]] zs = [z for z in (iue_forced_z(os.path.join(iue_ROOT, "mxhi", o + ".mxhi.gz"), p["lam"]) for o in other) if z is not None] p["other_aperture_z"] = [round(z, 1) for z in zs] p["aperture_fixed"] = bool(zs) and max(zs) < 2.0 # forced measurements in other stars below detection: rest wavelength, same camera (stellar); observed wavelength, same camera and aperture (detector) for p in left: imgs = [(s, os.path.join(iue_ROOT, "mxhi", o + ".mxhi.gz"), a) for (s, o), (c, a) in apert.items() if s != p["star"] and c == p["camera"]] p["forced_rest"] = sorted(uv_forced([(s, q) for s, q, _ in imgs], lambda s: p["rest"] * (1 + rvs.get(s, 0.0) / iue_C_KMS), iue_forced_z)) p["forced_obs"] = sorted(uv_forced([(s, q) for s, q, a in imgs if a == p["aperture"]], lambda s: p["lam"], iue_forced_z)) kept = [p for p in left if not p["aperture_fixed"] and len(p["forced_rest"]) < 2 and len(p["forced_obs"]) < 2] print("%d persistent features in at most one other star in either frame; %d of them absent from the same star's " "images through the other aperture; %d above 5 sigma in two or more other stars when measured there; " "%d left:" % (len(left), sum(p["aperture_fixed"] for p in left), sum(1 for p in left if not p["aperture_fixed"] and p not in kept), len(kept))) for p in kept[:60]: print(" %-20s %s %s %.3f A (rest %.3f), %d of %d images, peak %.1f sigma, other stars %d/%d, other aperture %s" % (p["star"][:20], p["camera"], p["aperture"], p["lam"], p["rest"], p["n"], p["of"], p["z"], p["stars_rest"], p["stars_obs"], p["other_aperture_z"] or "none")) json.dump(dict(images=n_read, peaks=len(feats), persistent=pers), open(os.path.join(iue_HERE, "iue_laser.json"), "w"), indent=1) def iue_main(): p = argparse.ArgumentParser() p.add_argument("mode", choices=("inventory", "fetch", "search")) p.add_argument("--threads", type=int, default=4) a = p.parse_args() {"inventory": iue_inventory, "fetch": lambda: iue_fetch(a.threads), "search": iue_search}[a.mode]() def iue_cli(): quiet_warnings() iue_main() # -------------------------------------------------------------------------------------------------------------- fuse_DOC = r"""Laser lines in FUSE spectra (905-1187 A) of stars within 20 pc: point-source peaks at p 1e-6 recurring in a second channel, away from stellar and airglow lines and other stars' peaks. python -I seti.py ultraviolet fuse inventory [--plx 50] python -I seti.py ultraviolet fuse fetch python -I seti.py ultraviolet fuse search""" fuse_ROOT = os.environ.get("FUSE_ROOT", core.data("fuse")) fuse_HERE = core.RESULTS fuse_MAST = "https://mast.stsci.edu/api/v0/invoke" fuse_DOWNLOAD = "https://mast.stsci.edu/api/v0.1/Download/file?uri=" def fuse_invoke(req): data = urllib.parse.urlencode({"request": json.dumps(req)}).encode() for attempt in range(3): try: return json.loads(urllib.request.urlopen(urllib.request.Request( fuse_MAST, data=data, headers=dict(core.UA, **{"Content-Type": "application/x-www-form-urlencoded"})), timeout=240).read()) except Exception: time.sleep(5 * (attempt + 1)) return {"data": []} def fuse_inventory(plx): os.makedirs(fuse_ROOT, exist_ok=True) stars = uv_nearby(plx) every, page = [], 1 while True: d = fuse_invoke({"service": "Mast.Caom.Filtered", "format": "json", "pagesize": 5000, "page": page, "params": {"columns": "obsid,obs_id,t_exptime,target_name,s_ra,s_dec,t_min", "filters": [{"paramName": "obs_collection", "values": ["FUSE"]}]}}) rows = d.get("data", []) every += rows if len(rows) < 5000: break page += 1 print("%d FUSE observations; %d SIMBAD objects with parallax above %g mas" % (len(every), len(stars), plx), flush=True) ora = np.array([float(x["s_ra"]) if x.get("s_ra") is not None else np.nan for x in every]) ode = np.array([float(x["s_dec"]) if x.get("s_dec") is not None else np.nan for x in every]) tmin = np.array([float(x["t_min"]) if x.get("t_min") is not None else 52500.0 for x in every]) obs = {} for s in stars: dist = uv_track_dist(s, ora, ode, tmin) for j in np.where(dist < 108.0)[0]: x = every[j] o = obs.setdefault(str(x["obsid"]), dict(x, stars=[])) o["stars"].append(dict(main_id=s["main_id"], plx=s["plx_value"], sp=s["sp_type"], otype=s["otype"], rv=s["rvz_radvel"], dist=round(float(dist[j]), 2))) files = [] ids = list(obs) for k in range(0, len(ids), 50): p = fuse_invoke({"service": "Mast.Caom.Products", "format": "json", "pagesize": 20000, "params": {"obsid": ",".join(ids[k:k + 50])}}) for x in p.get("data", []): fn = x["productFilename"] if fn.endswith("00000all4ttagfcal.fit.gz") or fn.endswith("00000all4histfcal.fit.gz"): files.append(dict(file=fn, uri=x["dataURI"], size=int(x.get("size") or 0), obsid=str(x.get("obsID")))) json.dump(dict(obs=obs, files=files), open(os.path.join(fuse_ROOT, "inventory.json"), "w"), indent=1) names = {s["main_id"] for o in obs.values() for s in o["stars"]} print("%d FUSE observations of %d objects within %.0f pc; %d combined files, %.0f MB" % (len(obs), len(names), 1000.0 / plx, len(files), sum(f["size"] for f in files) / 1e6)) def fuse_fetch(): inv = json.load(open(os.path.join(fuse_ROOT, "inventory.json"))) d = os.path.join(fuse_ROOT, "all") os.makedirs(d, exist_ok=True) def one(f): p = os.path.join(d, f["file"]) if os.path.exists(p): return 1 b = core.get(fuse_DOWNLOAD + f["uri"], timeout=300, retries=3, soft=True, pause=5) if not b: return 0 open(p + ".tmp", "wb").write(b) os.replace(p + ".tmp", p) return 1 with ThreadPoolExecutor(4) as ex: n = sum(ex.map(one, inv["files"])) print("%d of %d combined files on disk" % (n, len(inv["files"]))) fuse_C_KMS = 299792.458 fuse_LSF_SIG = 0.026 # A, point-source line spread function (FWHM about 0.06 A) fuse_THR = 5.0 fuse_TOL_KMS = 20.0 # km/s; channel wavelength zero points differ by up to this fuse_APERTURE_TOL = 15.0 # arcsec; the large square aperture is 30 arcsec across def fuse_owner(o): """Inventory entry of the star in an observation's aperture, or None.""" st = {" ".join(s["main_id"].split()): s for s in o.get("stars") or []} name = uv_owner({k: s["dist"] for k, s in st.items()}, fuse_APERTURE_TOL) return dict(st[name], main_id=name) if name else None def fuse_on_target(path): """False for observations that put the star at the reference point, outside every aperture.""" from astropy.io import fits return str(fits.getheader(path, 0).get("APERTURE", "")).strip().upper() != "RFPT" def fuse_channels(path): """{channel: (wave, flux, error)} of a combined file, pixels without an error dropped.""" from astropy.io import fits out = {} with fits.open(path) as h: for e in h[1:]: d = e.data if d is None or len(d) < 500: continue w, f, s = (np.asarray(d[c], float) for c in ("WAVE", "FLUX", "ERROR")) ok = np.isfinite(w) & np.isfinite(f) & np.isfinite(s) & (s > 0) if ok.sum() > 500: out[e.name.strip()] = (w[ok], f[ok], s[ok]) return out def fuse_matched(w, f, s): """Point-source line amplitude over its error at every pixel, above a 2 A running 90th percentile.""" from scipy.ndimage import percentile_filter step = float(np.median(np.diff(w))) sp = fuse_LSF_SIG / step x = np.arange(-int(4 * sp), int(4 * sp) + 1) k = np.exp(-0.5 * (x / sp) ** 2) r = f - percentile_filter(f, 90, size=int(2.0 / step) | 1, mode="nearest") iv = 1.0 / s ** 2 num = np.convolve(r * iv, k[::-1], mode="same") den = np.convolve(iv, (k ** 2)[::-1], mode="same") a = num / np.maximum(den, 1e-300) err = 1.0 / np.sqrt(np.maximum(den, 1e-300)) return a / err, r, err def fuse_line_limit(w, f, s): """Median THR-sigma flux of an unresolved line over the channel, erg cm-2 s-1.""" _, _, err = fuse_matched(w, f, s) return float(np.median(fuse_THR * err) * np.sqrt(2 * np.pi) * fuse_LSF_SIG) def fuse_channel_peaks(w, f, s, thr): """[(wavelength, z, fitted sigma A, amplitude)] of point-source peaks above thr.""" from scipy.optimize import curve_fit z, r, _ = fuse_matched(w, f, s) out = [] step = float(np.median(np.diff(w))) edge = int(1.0 / step) cand = np.nonzero(z > thr)[0] cand = cand[(cand > edge) & (cand < len(w) - edge)] used = set() for i in cand: if i in used: continue j = i + int(np.argmax(z[i:i + 12])) used.update(range(j - 12, j + 13)) sl = slice(max(j - 20, 0), j + 21) try: p, _ = curve_fit(lambda t, a, m, sg, b: a * np.exp(-0.5 * ((t - m) / sg) ** 2) + b, w[sl], r[sl], p0=(r[j], w[j], fuse_LSF_SIG, 0.0), sigma=s[sl], maxfev=3000) except Exception: continue a, m, sg, b = p sg = abs(sg) if a > 0 and abs(m - w[j]) < 0.05 and 0.5 * fuse_LSF_SIG <= sg <= 2.0 * fuse_LSF_SIG: out.append((float(m), float(z[j]), float(sg), float(a))) return out def fuse_search(): inv = json.load(open(os.path.join(fuse_ROOT, "inventory.json"))) rv = {} found, limits, n_read = [], {}, 0 for f in inv["files"]: p = os.path.join(fuse_ROOT, "all", f["file"]) own = fuse_owner(inv["obs"].get(f["obsid"]) or {}) if own is None or not os.path.exists(p) or not fuse_on_target(p): continue n_read += 1 star = own["main_id"] rv[star] = float(own["rv"]) if own.get("rv") not in (None, "") else 0.0 ch = fuse_channels(p) lim = [fuse_line_limit(*v) for v in ch.values()] if lim: plx = float(own.get("plx") or "nan") limits.setdefault(star, []).append(dict(file=f["file"], flux=float(np.median(lim)), plx=plx)) pk = {c: fuse_channel_peaks(*v, thr=fuse_THR) for c, v in ch.items()} low = {c: fuse_channel_peaks(*v, thr=3.0) for c, v in ch.items()} for c, lst in pk.items(): for lam, z, sg, a in lst: cover = [d for d, (w, _, _) in ch.items() if d != c and w[0] + 1 < lam < w[-1] - 1] hits = [d for d in cover if any(abs(l2 - lam) / lam * fuse_C_KMS < fuse_TOL_KMS for l2, _, _, _ in low[d])] found.append(dict(file=f["file"], star=star, rv=rv.get(star, 0.0), channel=c, lam=lam, z=z, sig=sg, amp=a, cover=cover, confirmed=hits)) multi = [x for x in found if x["confirmed"]] for x in multi: x["rest"] = x["lam"] / (1 + x["rv"] / fuse_C_KMS) # stars with a confirmed peak at the same wavelength, stellar and heliocentric frames for x in multi: x["stars_rest"] = len({y["star"] for y in multi if y["star"] != x["star"] and abs(y["rest"] - x["rest"]) / x["rest"] * fuse_C_KMS < fuse_TOL_KMS}) x["stars_obs"] = len({y["star"] for y in multi if y["star"] != x["star"] and abs(y["lam"] - x["lam"]) / x["lam"] * fuse_C_KMS < fuse_TOL_KMS}) left = [x for x in multi if x["stars_rest"] < 2 and x["stars_obs"] < 2] lam, strong = fuse_nist_strong() rng = np.random.default_rng(1) r = rng.uniform(910, 1185, 20000) k = np.searchsorted(lam, r).clip(1, len(lam) - 1) chance = float(np.mean(np.minimum(np.abs(lam[k] - r), np.abs(lam[k - 1] - r)) / r * fuse_C_KMS < fuse_TOL_KMS)) for x in left: dv = (lam - x["rest"]) / x["rest"] * fuse_C_KMS j = int(np.argmin(np.abs(dv))) x["id"] = "%s %.3f (%+.0f km/s)" % (strong[j]["ion"], strong[j]["lam"], dv[j]) if abs(dv[j]) < fuse_TOL_KMS else None dh = (lam - x["lam"]) / x["lam"] * fuse_C_KMS j = int(np.argmin(np.abs(dh))) if x["id"] is None and abs(dh[j]) < fuse_TOL_KMS and strong[j]["ion"] in fuse_AIRGLOW: x["id"] = "airglow %s %.3f (%+.0f km/s)" % (strong[j]["ion"], strong[j]["lam"], dh[j]) best = {} for star, ls in limits.items(): x = min(ls, key=lambda x: x["flux"]) d_cm = 1000.0 / x["plx"] * 3.0857e18 if x["plx"] > 0 else float("nan") best[star] = dict(flux=x["flux"], plx=x["plx"], lum=4 * np.pi * d_cm ** 2 * x["flux"], file=x["file"]) fl = np.array([b["flux"] for b in best.values()]) lu = np.array([b["lum"] for b in best.values() if np.isfinite(b["lum"])]) print("%d objects: median %g-sigma unresolved-line limit %.1e erg cm-2 s-1, isotropic %.1e erg/s (%.1e W)" % (len(best), fuse_THR, np.median(fl), np.median(lu), np.median(lu) * 1e-7)) json.dump(dict(files=n_read, objects=len(limits), found=found, left=left, chance=chance, limits=best), open(os.path.join(fuse_HERE, "fuse_laser.json"), "w"), indent=1) print("%d observations with a star in the aperture; %d point-source peaks above %g; %d confirmed in a second " "channel; %d not shared with two or more other stars; a random wavelength lies within %.0f km/s of a strong " "NIST line with probability %.2f" % (n_read, len(found), fuse_THR, len(multi), len(left), fuse_TOL_KMS, chance)) seen = set() for x in sorted(left, key=lambda x: -x["z"]): key = (x["star"], round(x["rest"], 1)) if key in seen: continue seen.add(key) print(" %-22s %s %.3f A (rest %.3f), z %.1f, sigma %.3f A, confirmed in %s of %s, other stars %d/%d: %s" % (x["star"], x["channel"], x["lam"], x["rest"], x["z"], x["sig"], x["confirmed"], x["cover"], x["stars_rest"], x["stars_obs"], x["id"] or "no line within %.0f km/s" % fuse_TOL_KMS)) print("%d distinct features have no strong line within %.0f km/s" % (len({(x["star"], round(x["rest"], 1)) for x in left if not x["id"]}), fuse_TOL_KMS)) fuse_IONS = ("H I", "He I", "He II", "C I", "C II", "C III", "C IV", "N I", "N II", "N III", "N IV", "N V", "O I", "O II", "O III", "O IV", "O V", "O VI", "Ne V", "Ne VI", "Si II", "Si III", "Si IV", "P II", "P IV", "P V", "S II", "S III", "S IV", "S V", "S VI", "Ar I", "Ar II", "Fe II", "Fe III", "Fe XVIII", "Fe XIX") fuse_AIRGLOW = ("H I", "O I", "N I", "N II", "He I") def fuse_nist_strong(): """NIST's strong FUV lines of chromospheric, transition-region, coronal and airglow ions (intensity listed or A above 1e7 /s, low levels for neutral and singly ionised).""" lines = [] for ion in fuse_IONS: lines += core.nist_lines(ion, 900, 1190, cache_dir=os.path.join(os.path.dirname(fuse_ROOT), "nist")) def num(x): try: return float(str(x).strip("[]()+x? ")) except ValueError: return float("nan") strong = sorted((ln for ln in lines if (ln["intens"].strip() or num(ln["aki"]) > 1e7) and (ln["ion"].split()[1] not in ("I", "II") or num(ln["ei"]) < 30000)), key=lambda ln: ln["lam"]) return np.array([ln["lam"] for ln in strong]), strong def fuse_main(): p = argparse.ArgumentParser() p.add_argument("mode", choices=("inventory", "fetch", "search")) p.add_argument("--plx", type=float, default=50.0) a = p.parse_args() if a.mode == "inventory": fuse_inventory(a.plx) elif a.mode == "fetch": fuse_fetch() else: fuse_search() def fuse_cli(): quiet_warnings() fuse_main() # -------------------------------------------------------------------------------------------------------------- stis_DOC = r"""Laser lines in post-StarCAT HST/STIS echelle x1d spectra (2004 July on) of stars whose tracks lie within 3 arcsec, through StarCAT's detector and tests. python -I seti.py ultraviolet stis inventory python -I seti.py ultraviolet stis fetch [--threads 4] python -I seti.py ultraviolet stis search (UV_PROCS processes) python -I seti.py ultraviolet stis identify (the features left, against NIST's lines)""" stis_ROOT = os.environ.get("STIS_X1D_ROOT", core.data("stis_x1d")) stis_HERE = core.RESULTS stis_T0 = 53200.0 # 2004 July, after StarCAT's last exposures stis_TOL = 3.0 # arcsec; the echelle apertures are 0.2 arcsec across stis_TOL_NAMED = 15.0 # arcsec from the track of the star the target name resolves to stis_NOT_STARS = ("QSO", "G", "GiG", "GiC", "GiP", "BiC", "Sy1", "Sy2", "Sy*", "AGN", "Bla", "BLL", "LIN", "SBG", "rG", "PN", "SNR", "HII", "Cl*", "GlC", "OpC", "ISM", "MoC", "Pl", "Pl?", "LSB", "EmG", "H2G", "bCG", "IG") def stis_inventory(): os.makedirs(stis_ROOT, exist_ok=True) every, page = [], 1 while True: d = hst_pulses_invoke({"service": "Mast.Caom.Filtered", "format": "json", "pagesize": 50000, "page": page, "params": {"columns": "obsid,obs_id,instrument_name,filters,target_name,s_ra,s_dec,t_min", "filters": [{"paramName": "obs_collection", "values": ["HST"]}, {"paramName": "instrument_name", "values": ["STIS/FUV-MAMA", "STIS/NUV-MAMA"]}, {"paramName": "filters", "values": list(starcat_RES)}, {"paramName": "t_min", "values": [{"min": stis_T0, "max": 99999}]}]}}) rows = d.get("data", []) every += rows if len(rows) < 50000: break page += 1 def finite(x): try: return math.isfinite(float(x)) except (TypeError, ValueError): return False every = [r for r in every if finite(r.get("s_ra")) and finite(r.get("s_dec"))] pts = {} for r in every: pts.setdefault((round(float(r["s_ra"]), 5), round(float(r["s_dec"]), 5)), []).append(r) rows = [dict(source_id="%d" % i, ra=k[0], dec=k[1]) for i, k in enumerate(pts)] keys = list(pts) m = core.xmatch(rows, "simbad", 30.0) cands = {} for x in m: if (x.get("otype") or "").strip() in stis_NOT_STARS: continue cands.setdefault(int(x["source_id"]), []).append(x) # a resolved target name credits its star within stis_TOL_NAMED, since binary orbits carry stars off their tracks tnames = sorted({str(r["target_name"]) for r in every if r.get("target_name")}) ident = {t: re.sub(r"^(HD|HR|HIP|GJ|GL|BD|CD|LHS|LTT|WD)(?=[0-9+-])", r"\1 ", t.replace("-", " "), flags=re.I) for t in tnames} named = uv_resolve(sorted(set(ident.values()))) target_star = {t: named.get(" ".join(v.split())) for t, v in ident.items()} stars = sorted({s for s in target_star.values() if s}) info = {} for j in range(0, len(stars), 50): q = ("SELECT main_id, ra, dec, pmra, pmdec, plx_value AS plx, sp_type, rvz_radvel AS radvel, otype FROM basic " "WHERE main_id IN (%s)" % ",".join("'%s'" % s.replace("'", "''") for s in stars[j:j + 50])) for x in core.simbad(q): if (x.get("otype") or "").strip() not in stis_NOT_STARS and x.get("ra") is not None: info[" ".join(str(x["main_id"]).split())] = x obs = [] for i, k in enumerate(keys): best = None for x in cands.get(i, []) + [info[s] for s in {target_star.get(str(r.get("target_name"))) for r in pts[k]} if s in info]: nm = " ".join(str(x["main_id"]).split()) in {target_star.get(str(r.get("target_name"))) for r in pts[k]} for r in pts[k]: dist = float(uv_track_dist(dict(ra=x["ra"], dec=x["dec"], pmra=x.get("pmra") or 0, pmdec=x.get("pmdec") or 0), [float(r["s_ra"])], [float(r["s_dec"])], [float(r["t_min"])])[0]) ok = dist <= (stis_TOL_NAMED if nm else stis_TOL) if ok and (best is None or (not nm, dist) < (not best[2], best[0])): best = (dist, x, nm) if best is None: continue x = best[1] for r in pts[k]: obs.append(dict(r, star=" ".join(str(x["main_id"]).split()), otype=x.get("otype"), sp=x.get("sp_type"), rv=x.get("radvel"), plx=x.get("plx"), dist=round(best[0], 2), named=best[2])) ids = sorted({str(o["obsid"]) for o in obs}) files = [] for k in range(0, len(ids), 50): p = hst_pulses_invoke({"service": "Mast.Caom.Products", "format": "json", "pagesize": 20000, "params": {"obsid": ",".join(ids[k:k + 50])}}) for x in p.get("data", []): if x["productFilename"].endswith("_x1d.fits"): files.append(dict(file=x["productFilename"], uri=x["dataURI"], size=int(x.get("size") or 0), obsid=str(x.get("obsID")), parent=str(x.get("parent_obsid")))) files = list({f["file"]: f for f in files}.values()) json.dump(dict(obs=obs, files=files), open(os.path.join(stis_ROOT, "inventory.json"), "w"), indent=1) print("%d STIS echelle observations since 2004 July; %d of %d stars; %d x1d files, %.1f GB" % (len(every), len(obs), len({o["star"] for o in obs}), len(files), sum(f["size"] for f in files) / 1e9)) def stis_fetch(threads): inv = json.load(open(os.path.join(stis_ROOT, "inventory.json"))) d = os.path.join(stis_ROOT, "x1d") os.makedirs(d, exist_ok=True) def one(f): p = os.path.join(d, f["file"]) if os.path.exists(p): return 1 b = core.get(hst_pulses_DOWNLOAD + f["uri"], timeout=300, retries=3, soft=True, pause=5) if not b or not b.startswith(b"SIMPLE"): return 0 open(p + ".tmp", "wb").write(b) os.replace(p + ".tmp", p) return 1 with ThreadPoolExecutor(threads) as ex: n = sum(ex.map(one, inv["files"])) print("%d of %d x1d files on disk" % (n, len(inv["files"]))) def stis_orders(path): """[(wavelength, flux, error)] of each echelle order at good pixels, with the exposure's metadata.""" from astropy.io import fits out = [] with fits.open(path) as h: hd = h[0].header meta = dict(mode=str(hd.get("OPT_ELEM", "")).strip(), cenwave=hd.get("CENWAVE"), file=os.path.basename(path), rootname=str(hd.get("ROOTNAME", "")).strip().lower(), target=str(hd.get("TARGNAME", "")).strip()) for e in h[1:]: if e.name != "SCI" or e.data is None: continue meta["mjd"] = e.header.get("EXPSTART", meta.get("mjd")) for row in e.data: w, f, s, dq = (np.asarray(row[c], float) for c in ("WAVELENGTH", "FLUX", "ERROR", "DQ")) good = (s > 0) & np.isfinite(f) & ((dq.astype(int) & (4 | 8 | 16 | 32 | 64 | 128 | 256)) == 0) if good.sum() > 200: out.append((w[good], f[good], s[good])) return out, meta def stis__feat(job): quiet_warnings() path, star = job try: od, meta = stis_orders(path) except Exception: return None if meta["mode"] not in starcat_RES: return None return dict(meta, star=star, feats=[x for w, f, e in od for x in starcat_detect(w, f, e, meta["mode"])]) def stis_forced(job): """(star, [(z at each rest wavelength, z at each observed wavelength in the same grating and central wavelength)]) for one exposure.""" quiet_warnings() path, star, rv, feats = job try: od, meta = stis_orders(path) except Exception: return None out = [] for p in feats: zr = zo = None for w, f, e in od: a = starcat_forced_z(w, f, e, p["rest"] * (1 + rv / starcat_C_KMS)) zr = a if zr is None or (a is not None and a > zr) else zr if meta["mode"] == p["mode"] and meta["cenwave"] in p["cenwaves"]: b = starcat_forced_z(w, f, e, p["lam"]) zo = b if zo is None or (b is not None and b > zo) else zo out.append((zr, zo)) return star, out def stis_search(): from concurrent.futures import ProcessPoolExecutor inv = json.load(open(os.path.join(stis_ROOT, "inventory.json"))) star_of = {str(o["obsid"]): o for o in inv["obs"]} jobs = [] for f in inv["files"]: o = star_of.get(f["parent"]) or star_of.get(f["obsid"]) p = os.path.join(stis_ROOT, "x1d", f["file"]) if o and os.path.exists(p): jobs.append((p, o["star"])) rv = {o["star"]: float(o["rv"]) if o.get("rv") not in (None, "") else 0.0 for o in inv["obs"]} procs = int(os.environ.get("UV_PROCS", "2")) with ProcessPoolExecutor(procs) as ex: res = [r for r in ex.map(stis__feat, jobs, chunksize=4) if r] print("%d exposures of %d stars; %d narrow peaks at %d sigma or more" % ( len(res), len({r["star"] for r in res}), sum(len(r["feats"]) for r in res), starcat_Z_MIN), flush=True) by = {} for r in res: by.setdefault((r["star"], r["mode"]), []).append(r) pers = [] for (star, mode), rs in by.items(): if len({r["rootname"] for r in rs}) < 2: continue allf = sorted(((x, r["rootname"], r["cenwave"]) for r in rs for x in r["feats"]), key=lambda t: t[0]["lam"]) for grp in uv_groups([x["lam"] for x, _, _ in allf], [x["inst"] / 3 for x, _, _ in allf]): files = {allf[j][1] for j in grp} if len(files) >= 2: xs = [allf[j][0] for j in grp] lam = float(np.mean([y["lam"] for y in xs])) pers.append(dict(star=star, mode=mode, lam=lam, rest=lam / (1 + rv.get(star, 0.0) / starcat_C_KMS), n=len(files), of=len({r["rootname"] for r in rs}), z=float(max(y["z"] for y in xs)), inst=xs[0]["inst"], cenwaves=sorted({allf[j][2] for j in grp}))) print("%d persistent features (two or more exposures in one grating) in %d stars" % (len(pers), len({p["star"] for p in pers})), flush=True) pk = [(x["lam"], r["star"]) for r in res for x in r["feats"]] obsw = np.array([lam for lam, _ in pk]) rest = np.array([lam / (1 + rv.get(s, 0.0) / starcat_C_KMS) for lam, s in pk]) who = np.array([s for _, s in pk]) for p in pers: p["stars_rest"] = len(set(who[np.abs(rest - p["rest"]) / p["rest"] * starcat_C_KMS < 15.0]) - {p["star"]}) p["stars_obs"] = len(set(who[np.abs(obsw - p["lam"]) / p["lam"] * starcat_C_KMS < 15.0]) - {p["star"]}) left = [p for p in pers if p["stars_rest"] < 2 and p["stars_obs"] < 2] hits_r, hits_o = [set() for _ in left], [set() for _ in left] with ProcessPoolExecutor(procs) as ex: for out in ex.map(stis_forced, [(p, s, rv.get(s, 0.0), left) for p, s in jobs], chunksize=4): if not out: continue s, zs = out for k, (zr, zo) in enumerate(zs): if s != left[k]["star"]: if zr is not None and zr > 5.0: hits_r[k].add(s) if zo is not None and zo > 5.0: hits_o[k].add(s) for k, p in enumerate(left): p["forced_rest"], p["forced_obs"] = sorted(hits_r[k]), sorted(hits_o[k]) unique = sorted([p for p in left if len(p["forced_rest"]) < 2 and len(p["forced_obs"]) < 2], key=lambda p: -p["z"]) print("%d persistent features in at most one other star in either frame; %d of them above 5 sigma in two or more " "other stars when measured there; %d left:" % (len(left), len(left) - len(unique), len(unique))) for p in unique[:60]: print(" %-24s %s %.3f A (rest %.3f), in %d of %d exposures, peak %.1f sigma, other stars %d/%d" % (p["star"][:24], p["mode"], p["lam"], p["rest"], p["n"], p["of"], p["z"], p["stars_rest"], p["stars_obs"])) json.dump(dict(exposures=len(res), objects=len({r["star"] for r in res}), peaks=sum(len(r["feats"]) for r in res), persistent=pers), open(os.path.join(stis_HERE, "stis_laser.json"), "w"), indent=1) def stis_identify(): """Surviving STIS features against NIST lines by cos identify's rule, with the chance a random wavelength lies as close.""" lines = [] for ion in cos_IONS: lines += core.nist_lines(ion, 1140, 3150, cache_dir=os.path.join(os.path.dirname(stis_ROOT), "nist")) def num(x): try: return float(str(x).strip("[]()+x? ")) except ValueError: return float("nan") strong = sorted((ln for ln in lines if (ln["intens"].strip() or num(ln["aki"]) > 1e7) and (ln["ion"].split()[1] not in ("I", "II") or num(ln["ei"]) < 30000)), key=lambda ln: ln["lam"]) lam = np.array([ln["lam"] for ln in strong]) path = os.path.join(stis_HERE, "stis_laser.json") d = json.load(open(path)) left = [p for p in d["persistent"] if "forced_rest" in p and len(p["forced_rest"]) < 2 and len(p["forced_obs"]) < 2] r = np.random.default_rng(1).uniform(1150, 3100, 20000) k = np.searchsorted(lam, r).clip(1, len(lam) - 1) near = np.minimum(np.abs(lam[k] - r), np.abs(lam[k - 1] - r)) / r * cos_C_KMS print("%d NIST lines, %d strong or low-excitation; a random wavelength lies within 15 km/s of one with probability " "%.2f (FUV %.2f, NUV %.2f)" % (len(lines), len(strong), np.mean(near < 15), np.mean(near[r < 1700] < 15), np.mean(near[r >= 1700] < 15))) unid = [] for p in sorted(left, key=lambda p: -p["z"]): dv = (lam - p["rest"]) / p["rest"] * cos_C_KMS j = int(np.argmin(np.abs(dv))) p["id"] = "%s %.3f (%+.0f km/s)" % (strong[j]["ion"], strong[j]["lam"], dv[j]) if abs(dv[j]) < 15 else None if not p["id"]: unid.append(p) by = {} for p in left: by.setdefault(p["star"], [0, 0])[0 if p["id"] else 1] += 1 print("%d features left in %d stars; %d without a line within 15 km/s" % (len(left), len(by), len(unid))) for s, (a, b) in sorted(by.items(), key=lambda kv: -sum(kv[1])): print(" %-24s identified %d, unidentified %d" % (s[:24], a, b)) for p in unid[:40]: print(" unidentified %-20s %s rest %.3f, in %d of %d exposures, %.1f sigma" % (p["star"][:20], p["mode"], p["rest"], p["n"], p["of"], p["z"])) json.dump(dict(d, identified=left), open(path, "w"), indent=1) def stis_main(): p = argparse.ArgumentParser() p.add_argument("mode", choices=("inventory", "fetch", "search", "identify")) p.add_argument("--threads", type=int, default=4) a = p.parse_args() if a.mode == "inventory": stis_inventory() elif a.mode == "fetch": stis_fetch(a.threads) elif a.mode == "identify": stis_identify() else: stis_search() def stis_cli(): quiet_warnings() stis_main() # -------------------------------------------------------------------------------------------------------------- limits_DOC = r"""Sensitivity of the ultraviolet laser searches: nominal unresolved-line thresholds per exposure and injection recovery, a star's limit m90 times its second-best exposure. python -I seti.py ultraviolet limits starcat|iue|fuse|stis [--n 300] [--procs 4] python -I seti.py ultraviolet limits cos [--n 300] [--procs 4] (where the COS spectra are, STIS_ROOT) python -I seti.py ultraviolet limits table (per-star limits, distances and powers) python -I seti.py ultraviolet limits record (the four searches and their limits, coverage store)""" limits_HERE = core.RESULTS limits_MULT = (1.0, 1.5, 2.0, 3.0, 5.0, 8.0) limits_PER = 6 # lines injected per exposure and multiple, 6 instrumental widths apart limits_ROOT2PI = math.sqrt(2 * math.pi) def limits_m_at(frac, mult, rec): """Multiple recovered with probability frac, interpolated in the recovery curve.""" rec = np.maximum.accumulate(np.asarray(rec, float)) if rec[-1] < frac: return float("inf") return float(np.interp(frac, rec, mult)) if rec[0] < frac else float(mult[0]) def limits_masked(lam, masks, widths): return any(abs(lam - a) < w for a, w in zip(masks, widths)) # ------------------------------------------------------------------------------------------------------- StarCAT def limits_starcat_one(args): path, seed = args x = starcat_load(path) if x is None or len(x[0]) < 1000: return None w, f, e, meta = x R = starcat_RES[meta["mode"]] sig = w / R / 2.3548 nom = starcat_Z_MIN * e * limits_ROOT2PI * sig free = np.array([not limits_masked(lam, starcat_AIRGLOW + starcat_CHROMOSPHERIC, [0.6] * len(starcat_AIRGLOW) + [1.2] * len(starcat_CHROMOSPHERIC)) for lam in w]) out = dict(file=meta["file"], obj=meta["obj"], mode=meta["mode"], nominal=float(np.median(nom[free])), lam_med=float(np.median(w)), ra=meta["ra"], dec=meta["dec"]) if seed is None: return out rng = np.random.default_rng(seed) lo, hi = w[0] + 10, w[-1] - 10 rec = [] for m in limits_MULT: inj = np.zeros_like(f) lams = [] for _ in range(limits_PER * 20): l0 = rng.uniform(lo, hi) j = int(np.searchsorted(w, l0)) if j <= 0 or j >= len(w) or not free[j] or any(abs(l0 - q) < 6 * l0 / R for q in lams): continue lams.append(l0) inj += m * nom[j] / (limits_ROOT2PI * sig[j]) * np.exp(-0.5 * ((w - l0) / sig[j]) ** 2) if len(lams) == limits_PER: break found = starcat_detect(w, f + inj, e, meta["mode"]) rec.append(sum(1 for l0 in lams if any(abs(q["lam"] - l0) < 0.5 * l0 / R for q in found)) / max(len(lams), 1)) out["rec"] = rec return out def limits_stis_one(args): """Nominal threshold and recovery of one post-2004 STIS exposure, one line per echelle order and multiple.""" path, seed, star = args try: od, meta = stis_orders(path) except Exception: return None if meta["mode"] not in starcat_RES or not od: return None R = starcat_RES[meta["mode"]] noms, lams = [], [] rng = np.random.default_rng(seed) if seed is not None else None hits, tries = np.zeros(len(limits_MULT)), np.zeros(len(limits_MULT)) masks = starcat_AIRGLOW + starcat_CHROMOSPHERIC widths = [0.6] * len(starcat_AIRGLOW) + [1.2] * len(starcat_CHROMOSPHERIC) for w, f, e in od: sig = w / R / 2.3548 nom = starcat_Z_MIN * e * limits_ROOT2PI * sig free = np.array([not limits_masked(lam, masks, widths) for lam in w]) noms += list(nom[free]) lams += list(w[free]) if rng is None or len(w) < 200: continue edge = min(3.0, 0.25 * (w[-1] - w[0])) # an E140H order spans about 4 A if edge <= 0: continue for k, m in enumerate(limits_MULT): l0 = rng.uniform(w[0] + edge, w[-1] - edge) j = int(np.searchsorted(w, l0)) if j <= 0 or j >= len(w) or not free[j]: continue inj = m * nom[j] / (limits_ROOT2PI * sig[j]) * np.exp(-0.5 * ((w - l0) / sig[j]) ** 2) found = starcat_detect(w, f + inj, e, meta["mode"]) tries[k] += 1 hits[k] += any(abs(q["lam"] - l0) < 0.5 * l0 / R for q in found) if not noms: return None out = dict(file=meta["file"], star=star, mode=meta["mode"], nominal=float(np.median(noms)), lam_med=float(np.median(lams))) if rng is not None: out["rec"] = (hits / np.maximum(tries, 1)).tolist() return out # ------------------------------------------------------------------------------------------------------------ IUE def limits_iue_one(args): path, seed = args try: od = iue_orders(path) except Exception: return None if not od: return None R = iue_R noms, lam_all = [], [] rng = np.random.default_rng(seed) if seed is not None else None hits = np.zeros(len(limits_MULT)) tries = np.zeros(len(limits_MULT)) for w, f, e, cal in od: okc = np.isfinite(cal) & (np.abs(f) > 5 * e) if okc.sum() < 10: continue conv = np.interp(w, w[okc], (cal / f)[okc]) # erg/cm2/s/A per flux number sig = w / R / 2.3548 nom = iue_Z_MIN * e * np.abs(conv) * limits_ROOT2PI * sig free = np.array([not limits_masked(lam, starcat_AIRGLOW + starcat_CHROMOSPHERIC, [0.6] * len(starcat_AIRGLOW) + [1.2] * len(starcat_CHROMOSPHERIC)) for lam in w]) noms += list(nom[free]) lam_all += list(w[free]) if rng is None or len(w) < 200: continue for k, m in enumerate(limits_MULT): l0 = rng.uniform(w[0] + 2, w[-1] - 2) j = int(np.searchsorted(w, l0)) if j <= 0 or j >= len(w) or not free[j]: continue amp = m * iue_Z_MIN * e[j] # peak in flux numbers found = iue_detect(w, f + amp * np.exp(-0.5 * ((w - l0) / sig[j]) ** 2), e) tries[k] += 1 hits[k] += any(abs(q["lam"] - l0) < 0.5 * l0 / R for q in found) if not noms: return None out = dict(file=os.path.basename(path), nominal=float(np.median(noms)), lam_med=float(np.median(lam_all))) if rng is not None: out["rec"] = (hits / np.maximum(tries, 1)).tolist() out["tries"] = tries.tolist() return out # ----------------------------------------------------------------------------------------------------------- FUSE def limits_fuse_threshold(ch, errs, l0): """Unresolved-line flux at l0 reaching THR in the most sensitive channel and 3 sigma in the next; None under two channels.""" cover = [c for c, (w, _, _) in ch.items() if w[0] + 1 < l0 < w[-1] - 1] if len(cover) < 2: return None, cover e = sorted(float(errs[c][min(int(np.searchsorted(ch[c][0], l0)), len(ch[c][0]) - 1)]) for c in cover) return max(fuse_THR * e[0], 3.0 * e[1]) * limits_ROOT2PI * fuse_LSF_SIG, cover def limits_fuse_one(args): path, seed = args if not fuse_on_target(path): return None ch = fuse_channels(path) if not ch: return None errs = {c: fuse_matched(*v)[2] for c, v in ch.items()} grid = np.arange(910.0, 1187.0, 0.5) noms = [x for x in (limits_fuse_threshold(ch, errs, lam)[0] for lam in grid) if x is not None and np.isfinite(x)] if not noms: return None out = dict(file=os.path.basename(path), nominal=float(np.median(noms)), lam_med=1040.0) if seed is None: return out rng = np.random.default_rng(seed) rec = [] for m in limits_MULT: hit = n = 0 for _ in range(limits_PER): c0 = list(ch)[rng.integers(len(ch))] w0 = ch[c0][0] l0 = rng.uniform(w0[0] + 2, w0[-1] - 2) F, cover = limits_fuse_threshold(ch, errs, l0) if F is None or not np.isfinite(F): continue F = m * F det5, det3 = [], [] for c in cover: w, f, s = ch[c] g = F / (limits_ROOT2PI * fuse_LSF_SIG) * np.exp(-0.5 * ((w - l0) / fuse_LSF_SIG) ** 2) if any(abs(x[0] - l0) / l0 * fuse_C_KMS < fuse_TOL_KMS for x in fuse_channel_peaks(w, f + g, s, thr=fuse_THR)): det5.append(c) if any(abs(x[0] - l0) / l0 * fuse_C_KMS < fuse_TOL_KMS for x in fuse_channel_peaks(w, f + g, s, thr=3.0)): det3.append(c) n += 1 hit += any(any(d != c for d in det3) for c in det5) rec.append(hit / max(n, 1)) out["rec"] = rec return out # ------------------------------------------------------------------------------------------------------------ COS def limits_cos_one(args): r, seed, sens = args from scipy.stats import poisson lam = np.array(r["lam"]) on, bk = np.array(r["on"]), np.array(r["bk"]) if len(lam) < 200 or r["exptime"] <= 0: return None from scipy.ndimage import percentile_filter s3 = np.convolve(on, np.ones(3), "same") b3 = np.convolve(bk, np.ones(3), "same") step = abs(float(np.median(np.diff(lam)))) env = percentile_filter(s3, 90, size=max(int(1.0 / step) * 2 + 1, 31), mode="nearest") mu = np.maximum(env, b3 + 0.5) # 3-bin counts for p < P_MAX, scaled by the 72 per cent of the line in 6 pixels thr = poisson.isf(cos_P_MAX, mu) + 1 sig_px = cos_LSF_PX / 2.3548 frac = math.erf(3.0 / (sig_px * math.sqrt(2))) nline = np.maximum(thr - mu, 1.0) / frac key = "%s/%s" % (r["grating"], r["cenwave"]) s = sens.get(key) if s is None: return None S = np.interp(lam, s["lam"], s["sens"], left=np.nan, right=np.nan) # counts/s/pixel per erg/cm2/s/A k = r["disp"] / (r["exptime"] * S) # erg/cm2/s per count free = np.array([not (any(abs(x - a) < w for a, w in cos_AIRGLOW) or any(abs(x - a) < 1.2 for a in cos_STRONG)) for x in lam]) & np.isfinite(k) if free.sum() < 100: return None nom = nline * k out = dict(file=r["file"], star=r["star"], grating=r["grating"], fppos=r["fppos"], nominal=float(np.median(nom[free])), lam_med=float(np.median(lam[free]))) if seed is None: return out rng = np.random.default_rng(seed) rec = [] px = np.arange(len(lam)) * cos_BIN_PX for m in limits_MULT: hit = n = 0 for _ in range(limits_PER): j = int(rng.integers(20, len(lam) - 20)) if not free[j]: continue prof = np.exp(-0.5 * ((px - px[j]) / sig_px) ** 2) prof = prof / prof.sum() add = rng.poisson(m * nline[j] * prof) rr = dict(r, on=(on + add).tolist()) found = cos_peaks(rr, 0, cos_P_MAX) n += 1 hit += any(abs(q["lam"] - lam[j]) < 3 * r["disp"] * cos_BIN_PX for q in found) rec.append(hit / max(n, 1)) out["rec"] = rec return out def limits_cos_sensitivity(): """NET/FLUX per pixel against wavelength per grating and central wavelength, from one x1d each.""" import urllib.request from astropy.io import fits root = os.environ.get("STIS_ROOT", core.data("stis_tag")) rs = cos_owned(json.load(open(os.path.join(root, "cos_spectra.json")))) path = os.path.join(root, "cos_sens.json") if os.path.exists(path): return json.load(open(path)), rs pick = {} for r in rs: pick.setdefault("%s/%s" % (r["grating"], r["cenwave"]), r["file"][:9]) os.makedirs(os.path.join(root, "x1d"), exist_ok=True) sens = {} for key, rootname in pick.items(): p = os.path.join(root, "x1d", rootname + "_x1d.fits") if not os.path.exists(p): u = "https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:HST/product/%s_x1d.fits" % rootname try: b = urllib.request.urlopen(urllib.request.Request(u, headers=core.UA), timeout=300).read() except Exception: continue if not b.startswith(b"SIMPLE"): continue open(p, "wb").write(b) lam, sv = [], [] with fits.open(p) as h: for row in h[1].data: w, fl, rate = (np.asarray(row[c], float) for c in ("WAVELENGTH", "FLUX", "NET")) ok = (fl > 0) & (rate > 0) & np.isfinite(fl) & np.isfinite(rate) lam += list(w[ok]) sv += list((rate / fl)[ok]) if lam: o = np.argsort(lam) lam, sv = np.array(lam)[o], np.array(sv)[o] # 51-pixel running median of the smooth ratio from scipy.ndimage import median_filter sens[key] = dict(lam=lam[::10].tolist(), sens=median_filter(sv, 51)[::10].tolist()) json.dump(sens, open(path, "w")) return sens, rs # ---------------------------------------------------------------------------------------------------------- driver def limits_run(which, n, procs): from concurrent.futures import ProcessPoolExecutor rng = np.random.default_rng(7) if which == "starcat": paths = sorted(os.path.join(dp, f) for dp, _, fs in os.walk(os.path.join(starcat_ROOT, "fits")) for f in fs if f.endswith("_spc.fits")) fn, items = limits_starcat_one, paths elif which == "iue": items = sorted(os.path.join(iue_ROOT, "mxhi", f) for f in os.listdir(os.path.join(iue_ROOT, "mxhi")) if f.endswith(".gz")) fn = limits_iue_one elif which == "fuse": items = sorted(os.path.join(fuse_ROOT, "all", f) for f in os.listdir(os.path.join(fuse_ROOT, "all"))) fn = limits_fuse_one elif which == "stis": inv = json.load(open(os.path.join(stis_ROOT, "inventory.json"))) star_of = {str(o["obsid"]): o["star"] for o in inv["obs"]} items = sorted((os.path.join(stis_ROOT, "x1d", f["file"]), star_of.get(f["parent"]) or star_of.get(f["obsid"])) for f in inv["files"] if os.path.exists(os.path.join(stis_ROOT, "x1d", f["file"]))) fn = limits_stis_one else: sens, rs = limits_cos_sensitivity() items = rs fn = limits_cos_one inj = set(rng.choice(len(items), size=min(n, len(items)), replace=False).tolist()) args = [(it, (int(rng.integers(1 << 30)) if i in inj else None)) for i, it in enumerate(items)] if which == "cos": args = [(it, s, sens) for it, s in args] elif which == "stis": args = [(it[0], s, it[1]) for it, s in args] with ProcessPoolExecutor(procs) as ex: res = [r for r in ex.map(fn, args, chunksize=4) if r] recs = np.array([r["rec"] for r in res if "rec" in r]) curve = recs.mean(axis=0).tolist() if len(recs) else [] # plateau: findable fraction of the band; the limit is 90 per cent of it plateau = float(max(curve)) if curve else None out = dict(set=which, mult=list(limits_MULT), recovery=curve, n_injected=int(len(recs)), plateau=plateau, m50=limits_m_at(0.5 * plateau, limits_MULT, curve) if curve else None, m90=limits_m_at(0.9 * plateau, limits_MULT, curve) if curve else None, exposures=res) json.dump(out, open(os.path.join(limits_HERE, "uv_limits_%s.json" % which), "w")) print("%s: %d exposures, median nominal limit %.2e erg cm-2 s-1; recovery at %s times nominal: %s; plateau %.2f; " "half and 90 per cent of it at %.2f and %.2f times nominal" % (which, len(res), np.median([r["nominal"] for r in res]), "/".join("%g" % m for m in limits_MULT), " ".join("%.2f" % x for x in curve), plateau or float("nan"), out["m50"] or float("nan"), out["m90"] or float("nan"))) limits_PC_M = 3.0857e16 limits_D_TX = 10.0 # transmitter aperture, m def limits_star_of_file(): """{file: SIMBAD main_id} for IUE images and FUSE observations with a star in the aperture.""" out = {} for o in iue_owned(json.load(open(os.path.join(iue_ROOT, "inventory.json")))): out[o["obs_id"].lower() + ".mxhi.gz"] = o["star"] fi = json.load(open(os.path.join(fuse_ROOT, "inventory.json"))) for f in fi["files"]: own = fuse_owner(fi["obs"].get(f["obsid"]) or {}) if own: out[f["file"]] = own["main_id"] return out def limits_table(): """Per-star limits: m90 times the second most sensitive exposure's nominal limit (best observation for FUSE), with distance, isotropic luminosity and 10 m diffraction-limited transmitter power.""" sets = {w: json.load(open(os.path.join(limits_HERE, "uv_limits_%s.json" % w))) for w in ("starcat", "iue", "cos", "fuse", "stis") if os.path.exists(os.path.join(limits_HERE, "uv_limits_%s.json" % w))} fmap = limits_star_of_file() plx, sp = {}, {} if "starcat" in sets: objs = {e["obj"]: (float(e["ra"]), float(e["dec"])) for e in sets["starcat"]["exposures"] if e.get("ra") is not None} for k, v in starcat_resolve(objs).items(): plx[k] = v["plx"] sp[k] = "%s; %s" % (v["main_id"], v["sp"]) names = set() for w in ("iue", "cos", "fuse", "stis"): for e in sets.get(w, {}).get("exposures", []): n = e.get("star") or fmap.get(e["file"]) if n: names.add(n) names = sorted(names) for k in range(0, len(names), 50): q = ("SELECT i.id, b.plx_value, b.sp_type FROM ident AS i JOIN basic AS b ON i.oidref = b.oid WHERE i.id IN (%s)" % ",".join("'%s'" % n.replace("'", "''") for n in names[k:k + 50])) for r in core.simbad(q): n = " ".join(str(r["id"]).split()) plx[n] = float(r["plx_value"]) if r.get("plx_value") not in (None, "") else None sp[n] = r.get("sp_type") or "" rows = [] for which, s in sets.items(): m90 = s["m90"] by = {} for e in s["exposures"]: star = e.get("obj") if which == "starcat" else e.get("star") or fmap.get(e["file"]) if not star: continue setting = e.get("mode") or e.get("grating") or e["file"][:3] by.setdefault((star, setting), []).append(e) best = {} for (star, setting), es in by.items(): v = sorted(x["nominal"] for x in es) lim = v[0] if which == "fuse" else (v[1] if len(v) >= 2 else None) if lim is None: continue lam = float(np.median([x["lam_med"] for x in es])) if star not in best or lim < best[star][0]: best[star] = (lim, setting, lam, len(es)) for star, (lim, setting, lam, n) in best.items(): F = lim * m90 # erg cm-2 s-1 p = plx.get(star) d = 1000.0 / p * limits_PC_M if p and p > 0 else None L = 4 * math.pi * d ** 2 * F * 1e-3 if d else None # W (erg cm-2 s-1 = 1e-3 W m-2) P = F * 1e-3 * d ** 2 * (lam * 1e-10 / limits_D_TX) ** 2 if d else None rows.append(dict(set=which, star=star, setting=setting, n=n, lam=lam, flux=F, plx=p, sp=sp.get(star, ""), L_iso_W=L, P_10m_W=P)) json.dump(dict(m90={w: s["m90"] for w, s in sets.items()}, rows=rows), open(os.path.join(limits_HERE, "uv_limits_table.json"), "w"), indent=1) for which in sets: rs = [r for r in rows if r["set"] == which] near = [r for r in rs if r["plx"] and r["plx"] > 100] print("%s: m90 %.2f; %d stars with two exposures, median %.1e erg cm-2 s-1; %d within 10 pc, median %.1e W " "isotropic, %.1e W through 10 m" % (which, sets[which]["m90"], len(rs), np.median([r["flux"] for r in rs]), len(near), np.median([r["L_iso_W"] for r in near]) if near else np.nan, np.median([r["P_10m_W"] for r in near]) if near else np.nan)) def limits_left(which, d): """Candidates a search's result file leaves after every test.""" if which in ("starcat", "stis"): return [p for p in d["persistent"] if "forced_rest" in p and len(p["forced_rest"]) < 2 and len(p["forced_obs"]) < 2] if which == "iue": return [p for p in d["persistent"] if "aperture_fixed" in p and not p["aperture_fixed"] and len(p["forced_rest"]) < 2 and len(p["forced_obs"]) < 2] if which == "cos": return [p for p in d["persistent"] if cos_left(p)] return d["left"] limits_SEARCHES = ( ("uv_laser_starcat", "starcat", "uv_laser.json", "Ultraviolet stellar spectra, STIS", "MAST HLSP StarCAT", "single exposures (Ayres 2010)", "every StarCAT exposure in E140M, E140H, E230M and E230H"), ("cos_laser_lines", "cos", "cos_laser.json", "Ultraviolet stellar spectra, COS", "MAST HST COS FUV corrtag", "CalCOS photon events", "COS FUV exposures with a SIMBAD star within 10 pc in the aperture"), ("iue_laser_lines", "iue", "iue_laser_6sig.json", "Ultraviolet stellar spectra, IUE", "MAST IUE NEWSIPS MXHI", "high dispersion", "images whose target is a SIMBAD star within 10 pc, by name or by its proper-motion track"), ("fuse_laser_lines", "fuse", "fuse_laser.json", "Ultraviolet stellar spectra, FUSE", "MAST FUSE CalFUSE combined channel spectra", "all4ttagfcal and all4histfcal", "observations with a SIMBAD star within 20 pc in the aperture, reference-point pointings left out"), ("uv_laser_stis", "stis", "stis_laser.json", "Ultraviolet stellar spectra, STIS after 2004", "MAST HST STIS x1d", "pipeline x1d, 2004 July on", "E140M, E140H, E230M and E230H exposures within 3 arcsec of a SIMBAD star's " "proper-motion track")) limits_VETTED = { "starcat": "steady features of Arcturus, in every exposure covering them from 1998 to 2003, emission of the Mira R " "Leo (the strongest within 4 km/s of Fe II 2626.27) and C III] 1908.73 of the symbiotic EG And", "cos": "chromospheric lines of eps Eri within 2 km/s of Cr II 1330.93, C I 1468.41 and Si I 1769.79, steady from " "2015 to 2025", "iue": "the last were fixed to one aperture's camera position or lay 440 km/s redward of Si II 1808 and 1817 A in " "the brightest G and K dwarfs, absent from STIS", "fuse": "of 9 single-channel peaks 4 lie on Ly beta, C III and O VI and 5 reach 5.0-5.5 sigma with nothing in the " "overlapping channel", "stis": "all lie in twelve stars whose ultraviolet spectra are crowded with natural lines: cool supergiants and giants " "(alpha Ori, alpha Her, gamma Cru, beta Gem), T Tauri stars (TW Hya, RU Lup), eta Car, the Ap star GY And, " "beta Pic, two stars catalogued as [KT2007] 1181 and 1184, and Procyon, whose features are resolved peaks " "of its spectrum steady in every exposure"} def limits_record(): """Record each search's row and per-star limits.""" t = json.load(open(os.path.join(limits_HERE, "uv_limits_table.json"))) for sid, which, fn, approach, archive, release, selection in limits_SEARCHES: d = json.load(open(os.path.join(limits_HERE, fn))) s = json.load(open(os.path.join(limits_HERE, "uv_limits_%s.json" % which))) rows = [r for r in t["rows"] if r["set"] == which] near = [r for r in rows if r["plx"] and r["plx"] > 100] left = limits_left(which, d) if which == "fuse": null = "%d single-channel peaks, %d confirmed in a second channel" % ( len(d["found"]), sum(1 for x in d["found"] if x["confirmed"])) else: null = "%d peaks, %d persistent, %d left after the cross-star, aperture and forced tests" % ( d["peaks"], len(d["persistent"]), len(left)) core.record_search( sid, approach, archive=archive, release=release, selection=selection, code=os.path.abspath(__file__), entry=which, n_targets=len(rows), null_summary=null, ledger="Ultraviolet stellar spectra", params=dict(plateau=s["plateau"], m50=s["m50"], m90=s["m90"], injected=s["n_injected"]), sensitivity="median per-star limit %.1e erg cm-2 s-1; %d stars within 10 pc, median %.1e W isotropic and " "%.1e W through a 10 m aperture" % (np.median([r["flux"] for r in rows]), len(near), np.median([r["L_iso_W"] for r in near]) if near else np.nan, np.median([r["P_10m_W"] for r in near]) if near else np.nan), result="%d features pass every test; %s" % (len(left), limits_VETTED[which])) core.clear_targets(sid) core.record_targets(sid, [ dict(target="simbad:%s" % (r["sp"].split(";")[0] if which == "starcat" and r.get("sp") else r["star"]), upper_limit=r["flux"], unit="erg cm-2 s-1, unresolved line at 90 per cent of the recovery plateau", note="%s, %.0f A; isotropic %.1e W; 10 m beam %.1e W" % (r["setting"], r["lam"], r["L_iso_W"] or np.nan, r["P_10m_W"] or np.nan)) for r in rows]) print("%s: %d stars, %d left" % (sid, len(rows), len(left))) def limits_main(): p = argparse.ArgumentParser() p.add_argument("which", choices=("starcat", "iue", "fuse", "cos", "stis", "table", "record")) p.add_argument("--n", type=int, default=300) p.add_argument("--procs", type=int, default=4) a = p.parse_args() if a.which == "table": limits_table() elif a.which == "record": limits_record() else: limits_run(a.which, a.n, a.procs) def limits_cli(): quiet_warnings() limits_main() # -------------------------------------------------------------------------------------------------------------- hst_pulses_DOC = r"""Ultraviolet laser pulses in STIS TIME-TAG and COS corrtag events of stars within 10 pc: same-bin photon clusters in 1 ms, 32 ms or 0.25 s beyond the local rate, Poisson FAP with trials. python -I seti.py ultraviolet hst_pulses inventory python -I seti.py ultraviolet hst_pulses fetch python -I seti.py ultraviolet hst_pulses search""" hst_pulses_ROOT = os.environ.get("STIS_ROOT", core.data("stis_tag")) hst_pulses_HERE = core.RESULTS hst_pulses_MAST = "https://mast.stsci.edu/api/v0/invoke" hst_pulses_DOWNLOAD = "https://mast.stsci.edu/api/v0.1/Download/file?uri=" hst_pulses_INSTRUMENTS = ["STIS/FUV-MAMA", "STIS/NUV-MAMA", "COS/FUV", "COS/NUV"] hst_pulses_MIN_SPAN = 10.0 hst_pulses_LOCAL = 5.0 hst_pulses_SPOT = 30.0 def hst_pulses_invoke(req): data = urllib.parse.urlencode({"request": json.dumps(req)}).encode() for attempt in range(3): try: return json.loads(urllib.request.urlopen(urllib.request.Request( hst_pulses_MAST, data=data, headers=dict(core.UA, **{"Content-Type": "application/x-www-form-urlencoded"})), timeout=240).read()) except Exception: time.sleep(5 * (attempt + 1)) return {"data": []} def hst_pulses_inventory(): os.makedirs(hst_pulses_ROOT, exist_ok=True) stars = uv_nearby(100) # one filtered MAST query, matched to the stars' tracks locally every, page = [], 1 while True: d = hst_pulses_invoke({"service": "Mast.Caom.Filtered", "format": "json", "pagesize": 50000, "page": page, "params": {"columns": "obsid,obs_id,instrument_name,t_exptime,target_name,s_ra,s_dec,t_min", "filters": [{"paramName": "obs_collection", "values": ["HST"]}, {"paramName": "instrument_name", "values": hst_pulses_INSTRUMENTS}]}}) rows = d.get("data", []) every += rows if len(rows) < 50000: break page += 1 print("%d HST observations with %s" % (len(every), ", ".join(hst_pulses_INSTRUMENTS)), flush=True) ora = np.array([float(x["s_ra"]) if x.get("s_ra") is not None else np.nan for x in every]) ode = np.array([float(x["s_dec"]) if x.get("s_dec") is not None else np.nan for x in every]) tmin = np.array([float(x["t_min"]) if x.get("t_min") is not None else 55000.0 for x in every]) allobs = {} res = [] for s in stars: dist = uv_track_dist(s, ora, ode, tmin) near = np.where(dist < 108.0)[0] res.append((s, near)) for j in near: x = every[j] o = allobs.setdefault(str(x["obsid"]), dict(x, stars=[], dist={})) o["stars"].append(s["main_id"]) o["dist"][" ".join(s["main_id"].split())] = round(float(dist[j]), 2) ids = list(allobs) tags = [] for k in range(0, len(ids), 50): p = hst_pulses_invoke({"service": "Mast.Caom.Products", "format": "json", "pagesize": 20000, "params": {"obsid": ",".join(ids[k:k + 50])}}) for x in p.get("data", []): fn = x["productFilename"] if fn.endswith(("_tag.fits", "_corrtag_a.fits", "_corrtag_b.fits", "_corrtag.fits")): tags.append(dict(file=fn, uri=x["dataURI"], size=int(x.get("size") or 0), obsid=str(x.get("obsID")), parent=str(x.get("parent_obsid")))) uniq = {t["file"]: t for t in tags} json.dump(dict(obs=allobs, tags=list(uniq.values())), open(os.path.join(hst_pulses_ROOT, "inventory.json"), "w"), indent=1) print("%d HST STIS and COS MAMA observations near %d objects; %d event files, %.1f GB" % (len(allobs), sum(1 for s, near in res if len(near)), len(uniq), sum(t["size"] for t in uniq.values()) / 1e9)) def hst_pulses_fetch(): inv = json.load(open(os.path.join(hst_pulses_ROOT, "inventory.json"))) os.makedirs(os.path.join(hst_pulses_ROOT, "events"), exist_ok=True) def one(t): p = os.path.join(hst_pulses_ROOT, "events", t["file"]) if os.path.exists(p): return 1 b = core.get(hst_pulses_DOWNLOAD + t["uri"], timeout=600, retries=3, soft=True, pause=5) if not b or not b.startswith(b"SIMPLE"): return 0 open(p + ".tmp", "wb").write(b) os.replace(p + ".tmp", p) return 1 with ThreadPoolExecutor(4) as ex: n = sum(ex.map(one, inv["tags"])) print("%d of %d event files on disk" % (n, len(inv["tags"]))) hst_pulses_CAL_TARGETS = ("DARK", "BIAS", "DEUTERIUM", "WAVE", "WAVEHITM", "WAVELINE", "FLAT", "CCDFLAT", "NONE") def hst_pulses_science(hd): """False for calibration exposures (darks, flats, lamps) filed under a science visit.""" et = str(hd.get("EXPTYPE", "EXTERNAL/SCI")).strip().upper() tn = str(hd.get("TARGNAME", "")).strip().upper() return et in ("EXTERNAL/SCI", "") and not any(tn.startswith(c) for c in hst_pulses_CAL_TARGETS) def hst_pulses_read_events(path): """Times (s), detector bin indices and test windows from a STIS _tag or COS _corrtag file.""" from astropy.io import fits out = [] with fits.open(path) as h: hd = h[0].header inst = str(hd.get("INSTRUME", "")).strip() if not hst_pulses_science(hd): return out for ext in h: if ext.name != "EVENTS" or ext.data is None or len(ext.data) == 0: continue d = ext.data cols = [c.upper() for c in d.columns.names] t = np.asarray(d["TIME"], float) if inst == "STIS": x, y = np.asarray(d["AXIS1"], float), np.asarray(d["AXIS2"], float) bx, by = 8, 8 windows = (0.001, 0.032, 0.25) # merge MAMA double counts: events within 2 pixels and 250 us of the previous o = np.argsort(t, kind="stable") t, x, y = t[o], x[o], y[o] dup = np.zeros(len(t), bool) for lag in range(1, 7): dup[lag:] |= (t[lag:] - t[:-lag] < 250e-6) & (np.abs(x[lag:] - x[:-lag]) <= 2) & \ (np.abs(y[lag:] - y[:-lag]) <= 2) t, x, y = t[~dup], x[~dup], y[~dup] else: x = np.asarray(d["XCORR"] if "XCORR" in cols else d["RAWX"], float) y = np.asarray(d["YCORR"] if "YCORR" in cols else d["RAWY"], float) fuv = str(hd.get("DETECTOR", "")).strip() == "FUV" bx, by = (16, 4) if fuv else (8, 8) good = np.ones(len(t), bool) if "DQ" in cols: # COS DQ bits (Data Handbook table 2.19): 2 hot spot, 8 calibration, 16/1024 low response, 32 background, 64 burst, 128 bounds, 512 PHA, 2048 bad time, 4096 low PHA, 8192 gain sag, 16384 dark edge good &= (np.asarray(d["DQ"], int) & (2 | 8 | 16 | 32 | 64 | 128 | 512 | 1024 | 2048 | 4096 | 8192 | 16384)) == 0 if fuv and "PHA" in cols: pha = np.asarray(d["PHA"], int) good &= (pha >= int(ext.header.get("PHALOWRA", 2))) & (pha <= int(ext.header.get("PHAUPPRA", 23))) if "WAVELENGTH" in cols: # no wavelength: digitizer-edge rows outside the calibrated area good &= np.asarray(d["WAVELENGTH"], float) > 0 # aperture light only: events within the SP_LOC, SP_HGT extraction band in YFULL if "YFULL" in cols: yf = np.asarray(d["YFULL"], float) xf = np.asarray(d["XFULL"], float) band = np.zeros(len(t), bool) for s in "ABC": loc, hgt = ext.header.get("SP_LOC_" + s), ext.header.get("SP_HGT_" + s) if loc is None or hgt is None or float(loc) < 0: continue if fuv and str(hd.get("SEGMENT", "")).strip()[-1:] != s: continue slp = float(ext.header.get("SP_SLP_" + s, 0.0) or 0.0) half = float(hgt) / 2 + 2 # slope reference column unrecorded: the trace is whichever of edge and centre puts more events in the band best = None for xr in (0.0, 8192.0 if fuv else 512.0): inb = np.abs(yf - (float(loc) + slp * (xf - xr))) <= half if best is None or inb.sum() > best.sum(): best = inb band |= best good &= band t, x, y = t[good], x[good], y[good] # TAGFLASH lamp intervals (LMP_ONn, LMPOFFn) cut with 2 s margins eh = ext.header flashes = [(float(eh["LMP_ON%d" % i]), float(eh.get("LMPOFF%d" % i, eh["LMP_ON%d" % i] + 15.0))) for i in range(1, 40) if "LMP_ON%d" % i in eh] if flashes: keep = np.ones(len(t), bool) for on, off in flashes: keep &= ~((t >= on - 2.0) & (t <= off + 2.0)) t, x, y = t[keep], x[keep], y[keep] windows = (0.032, 0.25) ok = np.isfinite(t) & np.isfinite(x) & np.isfinite(y) t, x, y = t[ok], x[ok], y[ok] b = (np.floor(x / bx).astype(np.int64) << 20) + np.floor(y / by).astype(np.int64) out.append(dict(inst=inst, target=str(hd.get("TARGNAME", "")).strip(), t=t, b=b, windows=windows, det=str(hd.get("DETECTOR", "")).strip() + str(hd.get("SEGMENT", "")).strip(), bx=bx, by=by, exptime=float(hd.get("TEXPTIME", 0) or (t.max() - t.min() if len(t) else 0)))) return out def hst_pulses_scan(ev): """Largest same-bin cluster per window, its bin's expectation and the false-alarm probability.""" from scipy.stats import poisson t, b = ev["t"], ev["b"] if len(t) < 50: return None t0 = float(t.min()) span = float(t.max()) - t0 if span < hst_pulses_MIN_SPAN: return None # one sort key by bin then time, bins over 2 x SPOT seconds apart key_scale = 10.0 ** math.ceil(math.log10(span + 2 * hst_pulses_SPOT + 10.0)) ub, inv, cnt = np.unique(b, return_inverse=True, return_counts=True) key = inv.astype(np.float64) * key_scale + (t - t0) o = np.argsort(key) key, inv_s, ts = key[o], inv[o], t[o] tall = np.sort(t) # expectation: the larger of the mean rate scaled by the detector's LOCAL rate (flares, bursts) and the spot's own SPOT-second rate (airglow, hot pixels) near = np.searchsorted(tall, ts + hst_pulses_LOCAL, side="right") - np.searchsorted(tall, ts - hst_pulses_LOCAL, side="left") scale = near / (2 * hst_pulses_LOCAL * len(t) / span) rate = cnt[inv_s] / span * np.maximum(scale, 1e-3) own = np.searchsorted(key, key + hst_pulses_SPOT, side="right") - np.searchsorted(key, key - hst_pulses_SPOT, side="left") dur = np.minimum(ts + hst_pulses_SPOT, t0 + span) - np.maximum(ts - hst_pulses_SPOT, t0) res = {} for w in ev["windows"]: # COS 32 ms ticks: a one-tick window holds one tick k = np.searchsorted(key, key + w * (1 - 1e-6), side="right") - np.arange(len(key)) mu = np.maximum(rate * w, (own - k) / np.maximum(dur - w, 1.0) * w) cand = np.where(k >= 2)[0] if len(cand) == 0: res[str(w)] = dict(k=1, mu=0.0, p=1.0, fap=1.0, t=None, detector=None, mean_detector=len(t) / span * w) continue # probability of the other k - 1 events at this bin's rate p = poisson.sf(k[cand] - 2, mu[cand]) j = cand[int(np.argmin(p))] pmin = float(p.min()) det = int(np.searchsorted(tall, ts[j] + w, side="right") - np.searchsorted(tall, ts[j], side="left")) fap = -math.expm1(len(t) * math.log1p(-min(pmin, 1 - 1e-16))) bj = int(ub[inv_s[j]]) res[str(w)] = dict(k=int(k[j]), mu=float(mu[j]), p=pmin, fap=fap, t=float(ts[j]), detector=det, mean_detector=len(t) / span * w, site=[ev["det"], (bj >> 20) * ev["bx"], (bj & 0xFFFFF) * ev["by"]]) return dict(n=int(len(t)), span=span, windows=res) def hst_pulses_scan_file(job): tg, stars = job p = os.path.join(hst_pulses_ROOT, "events", tg["file"]) if not os.path.exists(p): return [] try: evs = hst_pulses_read_events(p) except Exception as e: print(tg["file"], "unreadable:", str(e)[:60], flush=True) return [] out = [] for k, ev in enumerate(evs): r = hst_pulses_scan(ev) if r: out.append(dict(file=tg["file"], ext=k, inst=ev["inst"], target=ev["target"], stars=stars, **r)) return out def hst_pulses_search(): from concurrent.futures import ProcessPoolExecutor inv = json.load(open(os.path.join(hst_pulses_ROOT, "inventory.json"))) obs = inv["obs"] jobs = [(tg, obs.get(tg["parent"], obs.get(tg["obsid"], {})).get("stars", [])) for tg in inv["tags"]] with ProcessPoolExecutor(int(os.environ.get("STIS_PROCS", "2"))) as ex: rows = [r for rs in ex.map(hst_pulses_scan_file, jobs, chunksize=4) for r in rs] faps = [min(w["fap"] for w in r["windows"].values()) for r in rows] print("%d event lists (%d STIS, %d COS) of %d HST targets; %.1f million photons over %.1f days" % (len(rows), sum(r["inst"] == "STIS" for r in rows), sum(r["inst"] == "COS" for r in rows), len({r["target"] for r in rows}), sum(r["n"] for r in rows) / 1e6, sum(r["span"] for r in rows) / 86400)) for thr in (1e-2, 1e-3, 1e-5): print(" lists with a same-spot cluster at false-alarm probability below %.0e: %d" % (thr, sum(f < thr for f in faps))) # explanations: a flare raises the whole detector; a site hosting two targets' strongest clusters is a hot pixel top = [(r, f, min(r["windows"].items(), key=lambda kv: kv[1]["fap"])) for r, f in zip(rows, faps) if f < 1e-2] sites = {} for r, f, w in top: sites.setdefault(tuple(w[1]["site"]), set()).add(r["target"]) for r, f, w in top: x = w[1] rest, m = (x["detector"] or 0) - x["k"], x["mean_detector"] if rest > 3 * m + 5 * math.sqrt(m): r["why"] = "flare: %d other events detector-wide against a mean of %.1f" % (rest, m) elif len(sites[tuple(x["site"])]) > 1: r["why"] = "hot pixel: the same detector site (%s x %d y %d) hosts the strongest cluster of %s" % ( x["site"][0], x["site"][1], x["site"][2], ", ".join(sorted(sites[tuple(x["site"])]))) else: r["why"] = "" for r, f, w in sorted(top, key=lambda z: z[1])[:20]: print(" %s [%d] %s %-14s %s: %d photons over %.0f s; %s s window: %d in one spot (expected %.1e), FAP %.1e, at " "t %.4f with %d events detector-wide (mean %.2f)%s" % (r["file"], r["ext"], r["inst"], r["target"][:14], ",".join(s.strip() for s in r["stars"])[:30], r["n"], r["span"], w[0], w[1]["k"], w[1]["mu"], w[1]["fap"], w[1]["t"] or 0, w[1]["detector"] or 0, w[1]["mean_detector"], "; " + r["why"] if r["why"] else "")) left = [z for z in top if not z[0]["why"]] # smallest FAPs against the lists and windows tried trials = sum(len(r["windows"]) for r in rows) print("unexplained lists below 1e-3: %d (about %.1f expected from %d windows); below 1e-4: %d" % (sum(f < 1e-3 for r, f, w in left), trials * 1e-3, trials, sum(f < 1e-4 for r, f, w in left))) json.dump(rows, open(os.path.join(hst_pulses_HERE, "stis_pulses.json"), "w"), indent=1, default=str) def hst_pulses_record(): """Record per SIMBAD object from stis_pulses.json: lists, photons, hours and the smallest FAP left after flare and hot-pixel explanations.""" rows = json.load(open(os.path.join(hst_pulses_HERE, "stis_pulses.json"))) per = {} for r in rows: f = min(w["fap"] for w in r["windows"].values()) for s in r["stars"] or [r["target"]]: s = " ".join(s.split()) x = per.setdefault(s, dict(lists=0, n=0, span=0.0, fap=1.0, inst=set(), why=None)) x["lists"] += 1 x["n"] += r["n"] x["span"] += r["span"] x["inst"].add(r["inst"]) if f < x["fap"]: x["fap"], x["why"] = f, r.get("why") or None out = [dict(target="simbad:" + s, stat=x["n"], p=x["fap"], unit="photons; p is the smallest list false-alarm probability", vet=x["why"], note="%d %s event lists over %.1f h" % (x["lists"], "+".join(sorted(x["inst"])), x["span"] / 3600)) for s, x in per.items()] core.record_search("hst_photon_pulses", "HST photon-event pulses", archive="MAST HST STIS TIME-TAG and COS corrtag", release="1997-2024", selection="SIMBAD objects with parallax above 100 mas within 0.03 deg", code=os.path.abspath(__file__), entry="hst_pulses", params=dict(windows=[0.001, 0.032, 0.25], local=hst_pulses_LOCAL, spot=hst_pulses_SPOT, min_span=hst_pulses_MIN_SPAN), n_targets=len(out), result="%d unexplained lists below FAP 1e-3 against %.1f expected; %d below 1e-4" % (sum(1 for r in rows if not r.get("why") and min(w["fap"] for w in r["windows"].values()) < 1e-3), 1e-3 * sum(len(r["windows"]) for r in rows), sum(1 for r in rows if not r.get("why") and min(w["fap"] for w in r["windows"].values()) < 1e-4)), ledger="GALEX time-tagged photons") core.record_targets("hst_photon_pulses", out) print("recorded %d objects" % len(out)) def hst_pulses_main(): p = argparse.ArgumentParser() p.add_argument("mode", choices=("inventory", "fetch", "search", "record")) a = p.parse_args() {"inventory": hst_pulses_inventory, "fetch": hst_pulses_fetch, "search": hst_pulses_search, "record": hst_pulses_record}[a.mode]() def hst_pulses_cli(): quiet_warnings() hst_pulses_main() # -------------------------------------------------------------------------------------------------------------- galex_pulses_DOC = r"""Ultraviolet pulses from stars within 10 pc in GALEX time-tagged photons: largest counts in 5, 20 and 100 ms against the local rate, with control apertures and detector-wide bursts. python -I seti.py ultraviolet galex_pulses fetch [--band fuv|nuv] python -I seti.py ultraviolet galex_pulses search [--band fuv|nuv] python -I seti.py ultraviolet galex_pulses record [--band fuv|nuv]""" galex_pulses_ROOT = os.environ.get("GALEX_ROOT", core.data("galex")) galex_pulses_HERE = core.RESULTS galex_pulses_R_SRC = 8.0 galex_pulses_R_CTRL = 60.0 galex_pulses_N_CTRL = 8 galex_pulses_WINDOWS = (0.005, 0.02, 0.1) galex_pulses_LOCAL = 10.0 # band: photon-list suffix, directory of the cut apertures, results file galex_pulses_BANDS = {"fuv": ("fd", "cut", "galex_pulses.json"), "nuv": ("nd", "cut_nuv", "galex_pulses_nuv.json")} def galex_pulses_year_of(t): """GALEX time (seconds since 1980 January 6, the GPS epoch) as a decimal year.""" return 1980.0 + (t + 5 * 86400.0) / 3.15576e7 def galex_pulses_targets(): """Single-pointing eclipses and stars within 0.55 deg of each boresight, cached in ROOT/eclipse_stars.json.""" cache = os.path.join(galex_pulses_ROOT, "eclipse_stars.json") if os.path.exists(cache): return {int(e): v for e, v in json.load(open(cache)).items()} import pyarrow.parquet as pq b = pq.read_table(os.path.join(galex_pulses_ROOT, "boresight.parquet")).to_pandas() m = pq.read_table(os.path.join(galex_pulses_ROOT, "metadata.parquet")).to_pandas().set_index("eclipse") stars = core.simbad("SELECT main_id, ra, dec, pmra, pmdec, plx_value, otype, sp_type FROM basic WHERE plx_value > " "100 AND otype NOT IN ('Pl', 'Pl?') AND ra IS NOT NULL") ra0, de0 = np.radians(b["ra0"].values), np.radians(b["dec0"].values) by = {} for s in stars: a, d = math.radians(float(s["ra"])), math.radians(float(s["dec"])) cs = np.sin(d) * np.sin(de0) + np.cos(d) * np.cos(de0) * np.cos(ra0 - a) for e in set(b["eclipse"].values[np.degrees(np.arccos(np.clip(cs, -1, 1))) < 0.55]): e = int(e) if e in m.index and m.loc[e]["obstype"] != "AIS": by.setdefault(e, []).append(s) json.dump({str(e): v for e, v in by.items()}, open(cache, "w")) return by def galex_pulses_column(key, name): """One photon-list column by byte range from the public bucket; None when absent.""" import pyarrow.fs import pyarrow.parquet as pq fs = pyarrow.fs.S3FileSystem(anonymous=True, region="us-east-1", connect_timeout=60, request_timeout=600) for attempt in range(3): try: return pq.read_table("stpubdata/" + key, filesystem=fs, columns=[name]).column(name).to_numpy() except FileNotFoundError: return None except OSError as ex: if "NoSuchKey" in str(ex) or "does not exist" in str(ex): return None if attempt == 2: raise import time time.sleep(10) def galex_pulses_extract(e, stars, band="fuv"): suffix, cut, _ = galex_pulses_BANDS[band] out = os.path.join(galex_pulses_ROOT, cut, "e%05d.npz" % e) if os.path.exists(out): return 1 key = "galex/gap/galex_photon_list/e%05d/e%05d-%s.parquet" % (e, e, suffix) try: t = galex_pulses_column(key, "t") if t is None: return 0 good = (galex_pulses_column(key, "flags") == 0) & ~galex_pulses_column(key, "mask").astype(bool) de = galex_pulses_column(key, "dec") good &= np.isfinite(de) tt_all = t[good] if not len(tt_all): return 0 # keep photons in the stars' and controls' declination bands dy = galex_pulses_year_of(np.median(tt_all)) - 2000.0 band_sel = np.zeros(len(de), bool) for s in stars: sde = float(s["dec"]) + float(s.get("pmdec") or 0) / 3.6e6 * dy band_sel |= np.abs(de - sde) < (galex_pulses_R_CTRL + 2 * galex_pulses_R_SRC) / 3600 idx = np.nonzero(band_sel & good)[0] de = de[idx] ra = galex_pulses_column(key, "ra")[idx] tt = t[idx] del t, band_sel, good except Exception as ex: print("e%05d: %s" % (e, str(ex)[:120]), flush=True) return 0 res = dict(eclipse=e, band=band, t0=float(tt_all.min()), t1=float(tt_all.max()), n_all=int(len(tt_all))) # detector-wide photons per 5 ms tick, for field-wide bursts ticks, counts = np.unique(np.round(tt_all / 0.005).astype(np.int64), return_counts=True) del tt_all res["tick_hist"] = np.bincount(counts) arrays = {} for i, s in enumerate(stars): sra = float(s["ra"]) + float(s.get("pmra") or 0) / 3.6e6 * dy / math.cos(math.radians(float(s["dec"]))) sde = float(s["dec"]) + float(s.get("pmdec") or 0) / 3.6e6 * dy cd = math.cos(math.radians(sde)) for j in range(galex_pulses_N_CTRL + 1): if j == 0: cra, cde = sra, sde else: ang = 2 * math.pi * (j - 1) / galex_pulses_N_CTRL cra, cde = sra + galex_pulses_R_CTRL / 3600 * math.cos(ang) / cd, sde + galex_pulses_R_CTRL / 3600 * math.sin(ang) sel = (np.abs(de - cde) < galex_pulses_R_SRC / 3600) & (np.hypot(((ra - cra + 180) % 360 - 180) * cd, de - cde) * 3600 < galex_pulses_R_SRC) arrays["s%d_a%d" % (i, j)] = tt[sel] if j == 0: arrays["tick_s%d" % i] = counts[np.searchsorted(ticks, np.round(tt[sel] / 0.005).astype(np.int64))] os.makedirs(os.path.dirname(out), exist_ok=True) np.savez_compressed(out, meta=json.dumps(dict(res, tick_hist=res["tick_hist"].tolist(), stars=[s["main_id"] for s in stars])), **arrays) return 1 def galex_pulses_fetch(band): by = galex_pulses_targets() print("%d single-pointing eclipses cover %d objects within 10 pc" % (len(by), len({s["main_id"] for v in by.values() for s in v})), flush=True) done = [0] def one(kv): r = galex_pulses_extract(kv[0], kv[1], band) done[0] += 1 if done[0] % 10 == 0: print("%d of %d eclipses read" % (done[0], len(by)), flush=True) return r with ThreadPoolExecutor(int(os.environ.get("GALEX_THREADS", "3"))) as ex: n = sum(ex.map(one, sorted(by.items()))) print("%d eclipses with %s photon lists extracted" % (n, band)) def galex_pulses_clusters(t, span): """Per window length: largest photon-started count, its local expectation and the FAP over the exposure.""" from scipy.stats import poisson out = {} t = np.sort(t) for w in galex_pulses_WINDOWS: if len(t) < 2: out[w] = (len(t), 0.0, 1.0) continue k = np.searchsorted(t, t + w - 1e-9) - np.arange(len(t)) i = int(np.argmax(k)) n_local = np.searchsorted(t, t[i] + galex_pulses_LOCAL) - np.searchsorted(t, t[i] - galex_pulses_LOCAL) lam = max(n_local - k[i], 1) / (2 * galex_pulses_LOCAL) * w p = poisson.sf(k[i] - 1, lam) fap = -math.expm1(max(span / w, 1) * math.log1p(-min(p, 1 - 1e-16))) out[w] = (int(k[i]), float(lam), float(fap)) return out def galex_pulses_search(band): _, cut, result = galex_pulses_BANDS[band] files = sorted(os.listdir(os.path.join(galex_pulses_ROOT, cut))) rows, ctrl = [], [] for fn in files: z = np.load(os.path.join(galex_pulses_ROOT, cut, fn)) meta = json.loads(str(z["meta"])) span = meta["t1"] - meta["t0"] for i, name in enumerate(meta["stars"]): t = z["s%d_a0" % i] c = galex_pulses_clusters(t, span) best = min(c.values(), key=lambda x: x[2]) rows.append(dict(star=name, eclipse=meta["eclipse"], n=int(len(t)), span=span, clusters={str(k): v for k, v in c.items()}, fap=best[2], tick_max=int(z["tick_s%d" % i].max()) if len(t) else 0, tick_hist=meta["tick_hist"])) for j in range(1, galex_pulses_N_CTRL + 1): cc = galex_pulses_clusters(z["s%d_a%d" % (i, j)], span) ctrl.append(min(v[2] for v in cc.values())) ctrl = np.array(ctrl) print("%d star-eclipses (%d stars), median %d source photons over %.0f s; %d control apertures" % (len(rows), len({r["star"] for r in rows}), np.median([r["n"] for r in rows]), np.median([r["span"] for r in rows]), len(ctrl))) for thr in (1e-2, 1e-3, 1e-4): print(" false-alarm probability below %.0e: targets %d (%.1f expected from the controls' rate)" % (thr, sum(r["fap"] < thr for r in rows), len(rows) * np.mean(ctrl < thr))) for r in sorted(rows, key=lambda r: r["fap"])[:15]: c = r["clusters"] print(" %-24s e%05d: %d photons over %.0f s; 5 ms max %d (exp %.1e), 20 ms %d (exp %.1e), 100 ms %d (exp %.1e); " "FAP %.1e; busiest detector tick holding a source photon %d photons" % (r["star"], r["eclipse"], r["n"], r["span"], c["0.005"][0], c["0.005"][1], c["0.02"][0], c["0.02"][1], c["0.1"][0], c["0.1"][1], r["fap"], r["tick_max"])) json.dump(dict(rows=rows, ctrl=ctrl.tolist()), open(os.path.join(galex_pulses_HERE, result), "w"), indent=1) def galex_pulses_record(band): """Record one row per star with its lowest FAP over its eclipses.""" _, _, result = galex_pulses_BANDS[band] r = json.load(open(os.path.join(galex_pulses_HERE, result))) rows, ctrl = r["rows"], np.array(r["ctrl"]) best = {} for x in rows: if x["star"] not in best or x["fap"] < best[x["star"]]["fap"]: best[x["star"]] = dict(x, n_ecl=sum(1 for y in rows if y["star"] == x["star"])) run = "galex_%s_pulses" % band core.record_search( run, "GALEX time-tagged photons", archive="STScI S3 galex/gap/galex_photon_list (single-pointing eclipses)", release="GALEX %s photon lists" % band.upper(), selection="SIMBAD objects with parallax above 100 mas inside 0.55 deg of a single-pointing eclipse boresight", code="seti.py ultraviolet galex_pulses", params=dict(r_src=galex_pulses_R_SRC, n_ctrl=galex_pulses_N_CTRL, windows=list(galex_pulses_WINDOWS), local_s=galex_pulses_LOCAL), n_targets=len(best), null_summary=dict(below_1e2=sum(1 for x in rows if x["fap"] < 1e-2), below_1e3=sum(1 for x in rows if x["fap"] < 1e-3), expected_1e2=float(len(rows) * np.mean(ctrl < 1e-2)), controls=int(len(ctrl))), sensitivity="several photons within 5, 20 or 100 ms in an 8 arcsec aperture against the local rate", result="%d star-eclipses of %d stars: %d clusters at false-alarm probability below 0.01 against %.1f expected " "from %d control apertures" % (len(rows), len(best), sum(1 for x in rows if x["fap"] < 1e-2), len(rows) * np.mean(ctrl < 1e-2), len(ctrl)), ledger="GALEX time-tagged photons") core.clear_targets(run) core.record_targets(run, [ dict(target="simbad:%s" % s, p=float(x["fap"]), note="best of %d eclipses: e%05d, %d photons over %.0f s; 5 ms max %d" % (x["n_ecl"], x["eclipse"], x["n"], x["span"], x["clusters"]["0.005"][0])) for s, x in best.items()]) print("recorded") def galex_pulses_main(): p = argparse.ArgumentParser() p.add_argument("mode", choices=("fetch", "search", "record")) p.add_argument("--band", choices=tuple(galex_pulses_BANDS), default="fuv") a = p.parse_args() {"fetch": galex_pulses_fetch, "search": galex_pulses_search, "record": galex_pulses_record}[a.mode](a.band) def galex_pulses_cli(): quiet_warnings() galex_pulses_main() # -------------------------------------------------------------------------------------------------------------- galex_census_DOC = r"""Bursts across the GALEX NUV detector as the false-positive population: 100 ms clusters per 16 arcsec cell, judged off streaks and local rises, with a detector-position test. python -I seti.py ultraviolet galex_census census [--eclipses FILE] [--part k/n] python -I seti.py ultraviolet galex_census judge python -I seti.py ultraviolet galex_census detector ECLIPSE RA DEC T [--fixed] python -I seti.py ultraviolet galex_census burstpos ECLIPSE python -I seti.py ultraviolet galex_census record""" galex_census_ROOT = os.environ.get("GALEX_CENSUS_ROOT", core.data("galex_census")) galex_census_KEY = "stpubdata/galex/gap/galex_photon_list/e%05d/e%05d-nd.parquet" galex_census_W, galex_census_R, galex_census_FAP = 0.1, 16.0, 1e-3 def galex_census_photons(ecl, columns=("t", "ra", "dec")): """Good photons of one NUV eclipse from the public bucket (flags 0, unmasked, finite position).""" import pyarrow.fs import pyarrow.parquet as pq fs = pyarrow.fs.S3FileSystem(anonymous=True, region="us-east-1", connect_timeout=60, request_timeout=900) cols = list(dict.fromkeys(list(columns) + ["flags", "mask", "dec"])) tab = pq.read_table(galex_census_KEY % (ecl, ecl), filesystem=fs, columns=cols) good = ((tab.column("flags").to_numpy() == 0) & ~tab.column("mask").to_numpy().astype(bool) & np.isfinite(tab.column("dec").to_numpy())) return {c: tab.column(c).to_numpy()[good] for c in columns} def galex_census_eclipse(ecl): """Census of one eclipse into ROOT/e.json: every cluster below the FAP threshold.""" from scipy.stats import poisson W, R = galex_census_W, galex_census_R out_p = os.path.join(galex_census_ROOT, "e%05d.json" % ecl) if os.path.exists(out_p): return try: ph = galex_census_photons(ecl) except Exception as e: json.dump(dict(eclipse=ecl, error=str(e)[:200]), open(out_p, "w")) return t, ra, de = ph["t"], ph["ra"], ph["dec"] if len(t) < 1000: json.dump(dict(eclipse=ecl, error="too few photons"), open(out_p, "w")) return t0 = float(t.min()) tr = t - t0 span = float(tr.max()) ra0, de0 = float(np.median(ra)), float(np.median(de)) x = ((ra - ra0 + 180) % 360 - 180) * math.cos(math.radians(de0)) * 3600.0 y = (de - de0) * 3600.0 nwin = span / W hits, cells = [], {} for ox, oy in ((0, 0), (8, 0), (0, 8), (8, 8)): ix = np.floor((x + ox) / R).astype(np.int64) iy = np.floor((y + oy) / R).astype(np.int64) keyc = (ix - ix.min()) * 100000 + (iy - iy.min()) order = np.lexsort((tr, keyc)) kc, ts = keyc[order], tr[order] v = kc.astype(np.float64) * 4096.0 + ts cnt = np.searchsorted(v, v + W, side="left") - np.arange(len(v)) uniq, inv, ncell = np.unique(kc, return_inverse=True, return_counts=True) rate = ncell / span mu = rate * W if (ox, oy) == (0, 0): h, _ = np.histogram(rate, bins=[0, 0.1, 0.3, 1, 3, 10, 1e9]) cells = dict(zip(["<0.1", "0.1-0.3", "0.3-1", "1-3", "3-10", ">10"], h.tolist())) kthr = np.full(len(uniq), 10 ** 6, np.int64) for k in range(2, 600): ok = (kthr == 10 ** 6) & (poisson.sf(k - 1, mu) * nwin < galex_census_FAP) kthr[ok] = k if not (kthr == 10 ** 6).any(): break for h_ in np.nonzero(cnt >= kthr[inv])[0]: j = order[h_] fap = float(poisson.sf(cnt[h_] - 1, mu[inv[h_]]) * nwin) # the cell's own rate within 10 s either side, the window left out lo = np.searchsorted(v, kc[h_] * 4096.0 + max(ts[h_] - 10.0, 0.0)) hi = np.searchsorted(v, kc[h_] * 4096.0 + ts[h_] + 10.0) rate_loc = max(int(hi - lo - cnt[h_]), 1) / 20.0 fap_loc = float(poisson.sf(cnt[h_] - 1, max(rate_loc, rate[inv[h_]] * 0.5) * W) * nwin) hits.append(dict(t=float(ts[h_]), n=int(cnt[h_]), x=float(x[j]), y=float(y[j]), ra=float(ra[j]), dec=float(de[j]), rate=float(rate[inv[h_]]), rate_local=float(rate_loc), fap=fap, fap_local=fap_loc)) cl = [] for h_ in sorted(hits, key=lambda d: d["fap"]): if not any(abs(c["t"] - h_["t"]) < 0.3 and math.hypot(c["x"] - h_["x"], c["y"] - h_["y"]) < 24 for c in cl): cl.append(h_) for c in cl: c["neighbours"] = sum(1 for o in cl if o is not c and abs(o["t"] - c["t"]) < 2.0 and math.hypot(o["x"] - c["x"], o["y"] - c["y"]) < 60.0) json.dump(dict(eclipse=ecl, t0=t0, span=span, n_photons=int(len(t)), cells=cells, clusters=cl), open(out_p, "w")) def galex_census_eclipses(path=None): if path: return [int(e) for e in json.load(open(path))] rows = json.load(open(core.results("galex_pulses_nuv.json")))["rows"] return sorted({int(r["eclipse"]) for r in rows}) def galex_census_run(path=None, part=None): os.makedirs(galex_census_ROOT, exist_ok=True) ecls = galex_census_eclipses(path) k, n = (int(x) for x in part.split("/")) if part else (0, 1) for e in ecls[k::n]: t_ = time.time() galex_census_eclipse(e) print("e%05d done in %.0f s" % (e, time.time() - t_), flush=True) def galex_census_judge(): """Bursts of 7 or more photons in 100 ms on cells of 0.5-5 counts/s, off streaks and outside local rises.""" n_ecl = n_err = streaks = flares = 0 src_seconds = 0.0 bursts = [] for f in sorted(glob.glob(os.path.join(galex_census_ROOT, "e*.json"))): d = json.load(open(f)) if "error" in d: n_err += 1 continue n_ecl += 1 src_seconds += (d["cells"].get("1-3", 0) + d["cells"].get("3-10", 0)) * d["span"] for c in d["clusters"]: if c["n"] < 7 or not 0.5 <= c["rate"] <= 5.0: continue if c.get("neighbours", 0) > 0: streaks += 1 elif c.get("fap_local", 1.0) >= 1e-4: flares += 1 else: bursts.append(dict(c, eclipse=d["eclipse"])) # one burst per cell and moment across the four grids uniq = [] for b in sorted(bursts, key=lambda b: b["fap_local"]): if not any(b["eclipse"] == u["eclipse"] and abs(b["t"] - u["t"]) < 0.2 and abs(b["ra"] - u["ra"]) < 0.01 and abs(b["dec"] - u["dec"]) < 0.01 for u in uniq): uniq.append(b) print("%d eclipses judged (%d failed); %.3g source-seconds at 1-10 counts/s" % (n_ecl, n_err, src_seconds)) print("clusters of 7+ photons in 100 ms on 0.5-5 counts/s cells: %d on streaks, %d within a local rise, %d isolated" % (streaks, flares, len(uniq))) print("isolated bursts per source-second: %.2g" % (len(uniq) / max(src_seconds, 1))) for b in uniq[:30]: print(" e%05d t %.2f n %d at RA %.5f Dec %+.5f, cell %.2f c/s, local %.2f c/s, fap %.1e, local fap %.1e" % (b["eclipse"], b["t"], b["n"], b["ra"], b["dec"], b["rate"], b["rate_local"], b["fap"], b["fap_local"])) json.dump(dict(n_ecl=n_ecl, src_seconds=src_seconds, streaks=streaks, flares=flares, bursts=uniq), open(os.path.join(galex_census_ROOT, "judged.json"), "w")) def galex_census_boresight(tr, ra, de): """Boresight per photon from each second's median position (over 100 photons), 5 s running median, linear within the second; with cos(dec).""" from scipy.ndimage import median_filter sec = np.floor(tr).astype(int) ns = sec.max() + 1 order = np.argsort(sec) st = np.searchsorted(sec[order], np.arange(ns)) en = np.searchsorted(sec[order], np.arange(ns), side="right") bra, bde = np.full(ns, np.nan), np.full(ns, np.nan) for s in range(ns): idx = order[st[s]:en[s]] if len(idx) > 100: bra[s], bde[s] = np.median(ra[idx]), np.median(de[idx]) ok = np.isfinite(bra) sb_ra = np.interp(np.arange(ns), np.nonzero(ok)[0], median_filter(bra[ok], 5)) sb_de = np.interp(np.arange(ns), np.nonzero(ok)[0], median_filter(bde[ok], 5)) fr = tr - sec bore_ra = sb_ra[sec] + (np.r_[sb_ra[1:], sb_ra[-1]][sec] - sb_ra[sec]) * fr bore_de = sb_de[sec] + (np.r_[sb_de[1:], sb_de[-1]][sec] - sb_de[sec]) * fr return bore_ra, bore_de, math.cos(math.radians(np.nanmedian(bde))), (sb_ra, sb_de) def galex_census_detector(ecl, star_ra, star_de, t_burst, fixed=False): """A star's burst in detector coordinates, and that detector spot over the whole visit against the star.""" ph = galex_census_photons(ecl) t, ra, de = ph["t"], ph["ra"], ph["dec"] tr = t - t.min() bore_ra, bore_de, cd, (sb_ra, sb_de) = galex_census_boresight(tr, ra, de) dx = ((ra - bore_ra + 180) % 360 - 180) * cd * 3600.0 dy = (de - bore_de) * 3600.0 print("visit %.0f s, %d photons; boresight dither range RA %.1f arcsec, Dec %.1f arcsec" % (tr.max(), len(t), (sb_ra.max() - sb_ra.min()) * cd * 3600, (sb_de.max() - sb_de.min()) * 3600)) sp = np.hypot(np.diff(sb_ra) * cd * 3600, np.diff(sb_de) * 3600) print("boresight speed: median %.2f arcsec/s, 95th percentile %.2f" % (np.median(sp), np.percentile(sp, 95))) # star centroid: densest 2 arcsec cell within 60 arcsec, then a 6 arcsec mean x60 = ((ra - star_ra + 180) % 360 - 180) * cd * 3600 y60 = (de - star_de) * 3600 w60 = np.hypot(x60, y60) < 60 H, xe, ye = np.histogram2d(x60[w60], y60[w60], bins=60, range=[[-60, 60], [-60, 60]]) ix, iy = np.unravel_index(np.argmax(H), H.shape) cx, cy = (xe[ix] + xe[ix + 1]) / 2, (ye[iy] + ye[iy + 1]) / 2 m = w60 & (np.hypot(x60 - cx, y60 - cy) < 6) if not fixed: star_ra += np.mean(x60[m]) / 3600 / cd star_de += np.mean(y60[m]) / 3600 print("star centroid RA %.5f Dec %.5f from %d photons within 6 arcsec" % (star_ra, star_de, m.sum())) r8 = np.hypot(((ra - star_ra + 180) % 360 - 180) * cd * 3600, (de - star_de) * 3600) < 8 cand = np.nonzero(r8 & (tr > t_burst - 5) & (tr < t_burst + 5))[0] tt = tr[cand] best = None for x in tt: k = int(np.sum((tt >= x) & (tt < x + 0.1))) if best is None or k > best[0]: best = (k, x) t_best = best[1] sel = cand[(tt >= t_best) & (tt < t_best + 0.1)] px, py = dx[sel].mean(), dy[sel].mean() print("burst: %d photons from %.3f s; sky centre RA %.5f Dec %.5f; detector centre %.1f %.1f arcsec, spread %.1f arcsec" % (len(sel), t_best, ra[sel].mean(), de[sel].mean(), px, py, np.hypot(dx[sel] - px, dy[sel] - py).mean())) on_star = np.hypot(((ra - star_ra + 180) % 360 - 180) * cd * 3600, (de - star_de) * 3600) < 6 on_spot = np.hypot(dx - px, dy - py) < 6 star_det_x = ((star_ra - bore_ra + 180) % 360 - 180) * cd * 3600 star_det_y = (star_de - bore_de) * 3600 sep_t = np.hypot(star_det_x - px, star_det_y - py) print("star (sky, 6 arcsec): %d photons; detector spot (6 arcsec): %d photons; spot while the star is 20+ arcsec " "away: %d photons" % (on_star.sum(), on_spot.sum(), (on_spot & (sep_t > 20)).sum())) print("star track on the detector: x %.0f..%.0f, y %.0f..%.0f arcsec; burst spot %.1f %.1f; nearest approach of the " "track to the spot %.1f arcsec" % (star_det_x.min(), star_det_x.max(), star_det_y.min(), star_det_y.max(), px, py, sep_t.min())) def tops(mask, label): ts = np.sort(tr[mask]) idx = np.argsort(tr[mask]) cnt = np.searchsorted(ts, ts + 0.1) - np.arange(len(ts)) seen = [] for i in np.argsort(-cnt): if cnt[i] < 4: break if any(abs(ts[i] - s) < 0.2 for s in seen): continue seen.append(ts[i]) j = np.nonzero(mask)[0][idx[i]] print(" %s: %d photons in 100 ms at %.3f s; star-spot separation on the detector %.1f arcsec" % (label, cnt[i], ts[i], sep_t[j])) tops(on_spot, "detector spot") tops(on_star, "star") away = on_spot & (sep_t > 20) tops(away, "spot, star away") frac = np.mean(sep_t > 20) print("fraction of the visit the star is 20+ arcsec from the spot: %.2f; spot rate then %.3f c/s; background in " "6 arcsec %.3f c/s" % (frac, away.sum() / (frac * tr.max()), np.sum(np.hypot(dx - px - 120, dy - py) < 6) / tr.max())) def galex_census_burstpos(ecl): """Detector position and radius of every judged burst of one eclipse.""" d = json.load(open(os.path.join(galex_census_ROOT, "judged.json"))) ph = galex_census_photons(ecl, ("t", "ra", "dec", "detrad")) t, ra, de, dr = ph["t"], ph["ra"], ph["dec"], ph["detrad"] tr = t - t.min() bore_ra, bore_de, cd, _ = galex_census_boresight(tr, ra, de) for b in sorted((b for b in d["bursts"] if b["eclipse"] == ecl), key=lambda b: b["t"]): w = ((tr >= b["t"] - 0.02) & (tr < b["t"] + 0.12) & (np.hypot(((ra - b["ra"] + 180) % 360 - 180) * cd * 3600, (de - b["dec"]) * 3600) < 12)) dx = ((ra[w] - bore_ra[w] + 180) % 360 - 180) * cd * 3600 dy = (de[w] - bore_de[w]) * 3600 print("t %7.2f n %2d sky %.5f %+.5f -> detector %7.1f %7.1f arcsec (spread %.1f), detrad %.4f deg" % (b["t"], w.sum(), b["ra"], b["dec"], dx.mean(), dy.mean(), np.hypot(dx - dx.mean(), dy - dy.mean()).mean(), np.median(dr[w]))) def galex_census_record(): """Record the judged census: one row per isolated burst.""" d = json.load(open(os.path.join(galex_census_ROOT, "judged.json"))) b = d["bursts"] core.record_search( "galex_burst_census", "GALEX NUV burst census", archive="STScI S3 galex/gap/galex_photon_list (single-pointing eclipses)", release="%d NUV eclipses" % d["n_ecl"], selection="every 16 arcsec cell of the eclipses of the NUV pulse search", code=os.path.abspath(__file__), entry="galex_census", n_targets=len(b), params=dict(window_s=galex_census_W, cell_arcsec=galex_census_R, fap=galex_census_FAP, min_photons=7, cell_rate=[0.5, 5.0], local_fap=1e-4), null_summary=dict(streaks=d["streaks"], local_rises=d["flares"], source_seconds=d["src_seconds"]), sensitivity="7 or more photons in 100 ms on cells of 0.5 to 5 counts/s", result="%d eclipses: %d isolated bursts of 7 or more photons in 100 ms (%.2g per source-second), %d clusters on " "satellite streaks and %d within local rises" % (d["n_ecl"], len(b), len(b) / max(d["src_seconds"], 1), d["streaks"], d["flares"]), ledger="Ultraviolet pulses") core.clear_targets("galex_burst_census") core.record_targets("galex_burst_census", [dict( target="galex:e%05d t %.2f" % (x["eclipse"], x["t"]), stat=x["n"], p=x["fap_local"], unit="photons in 100 ms", note="RA %.5f Dec %+.5f, cell %.2f counts/s, local %.2f counts/s" % (x["ra"], x["dec"], x["rate"], x["rate_local"])) for x in b]) print("recorded %d bursts" % len(b)) def galex_census_cli(): quiet_warnings() p = argparse.ArgumentParser() p.add_argument("mode", choices=("census", "judge", "detector", "burstpos", "record")) p.add_argument("args", nargs="*") p.add_argument("--eclipses", default=None) p.add_argument("--part", default=None) p.add_argument("--fixed", action="store_true") a = p.parse_args() if a.mode == "census": galex_census_run(a.eclipses, a.part) elif a.mode == "judge": galex_census_judge() elif a.mode == "detector": e, ra_, de_, tb = a.args galex_census_detector(int(e), float(ra_), float(de_), float(tb), a.fixed) elif a.mode == "burstpos": galex_census_burstpos(int(a.args[0])) else: galex_census_record() # -------------------------------------------------------------------------------------------------------------- # ultraviolet from Earth's transit zone etzuv_DOC = r"""Ultraviolet from Earth's transit zone: Gaia Catalogue of Nearby Stars objects within 0.264 deg of the ecliptic, which see Earth cross the Sun, carried to 2007.0, against GALEX GUVcat AIS (II/335) within 5 arcsec by CDS XMatch, a detection in either band; the share with an ultraviolet source in the zone against the share 1-5 deg from the ecliptic and against the zone moved 2 arcmin in eight directions. Matches read in SIMBAD. python -I seti.py ultraviolet etzuv [record] [--radius-arcsec 5] [--epoch 2007] [--catalogue II/335/galex_ais] Another CDS catalogue of ultraviolet sources, with the epoch of its positions, is recorded under its own run.""" def etzuv_main(record=False, radius_arcsec=5.0, epoch=2007.0, catalogue="II/335/galex_ais"): from setitools.neighbourhood import etz_share, etz_vets r = etz_share(catalogue, radius_arcsec, epoch) vets = etz_vets(r["got"]) mean = float(np.mean(r["null"])) print("%d zone stars and %d stars 1-5 deg away; %s within %g arcsec: zone %d, away %d (%.2f expected in the zone, p %.3g); zone moved 2 arcmin %s (mean %.1f); " "zone matches by SIMBAD type %s" % (r["zone"], r["near"], catalogue, radius_arcsec, len(r["got"]), r["far"], r["exp"], r["p"], r["null"], mean, sorted({(v[0] or (None, "none"))[1] for v in vets.values()}))) if record: galex = catalogue == "II/335/galex_ais" core.record_run("etz_galex_stars" if galex else "etz_%s_stars" % catalogue.replace("/", "_"), "Ultraviolet from Earth's transit zone", "Gaia Catalogue of Nearby Stars (external.gaiaedr3_gcns_main_1); %s by CDS XMatch; SIMBAD" % ("GALEX GUVcat AIS (II/335)" if galex else catalogue), "stars within 100 pc and 0.264 deg of the ecliptic: %d" % r["zone"], os.path.abspath(__file__), "etzuv", r["zone"], "stars", dict(zone_deg=0.264, control_deg=[1, 5], radius_arcsec=radius_arcsec, epoch=epoch, catalogue=catalogue), "the share 1-5 deg away (%d of %d); the zone moved 2 arcmin: %s" % (r["far"], r["near"], r["null"]), "ultraviolet from a star that sees Earth transit, beyond the rate of stars beside the zone", "%d against %.2f from the share beside the zone (p %.3g) and a mean of %.1f moved" % (len(r["got"]), r["exp"], r["p"], mean), [dict(target="gaia_dr3:%s" % n[1:], vet=vets[n][1] or ("chance-rate" if r["p"] > 0.01 else None), note="SIMBAD %s; NUV %s, FUV %s at %.1f arcsec" % ( vets[n][0], m.get("NUVmag"), m.get("FUVmag"), float(m["angDist"]))) for n, m in sorted(r["got"].items())]) print("recorded") bluest_DOC = r"""Nearby light bluer than any star: Gaia Catalogue of Nearby Stars objects (probability above 0.9) whose BP-RP lies below -0.65, bluer than the Rayleigh-Jeans tail of any photosphere gives in Gaia's bands, a spectrum rising toward the blue faster than heat allows, each read for its BP/RP flux excess against the single-star relation, its white-dwarf probability and its SIMBAD type; the hottest white dwarfs and subdwarfs reach -0.6, and a blend or a faint source's BP noise makes a bluer colour without a bluer star. python -I seti.py ultraviolet bluest [record] [--colour-max -0.65] [--excess-max 0.2]""" def bluest_main(record=False, colour_max=-0.65, excess_max=0.2): from setitools import catalogues from setitools.core import num rows = catalogues.gcns("source_id, ra, dec, parallax, phot_g_mean_mag, phot_bp_mean_mag, phot_rp_mean_mag, phot_bp_rp_excess_factor, ruwe, wd_prob, ipd_frac_multi_peak", "gcns_prob > 0.9 AND phot_bp_mean_mag - phot_rp_mean_mag < %g" % colour_max, ranges=12) sim = {} ids = [str(r["source_id"]) for r in rows] for k in range(0, len(ids), 200): for s in core.simbad("SELECT i.id, b.main_id, b.otype, b.sp_type FROM ident AS i JOIN basic AS b ON i.oidref = b.oid WHERE i.id IN (%s)" % ",".join( "'Gaia DR3 %s'" % x for x in ids[k:k + 200])): sim[str(s["id"]).split()[-1]] = (" ".join(str(s["main_id"]).split()), str(s["otype"]).strip(), (s["sp_type"] or "").strip()) out = [] for r in rows: c = num(r["phot_bp_mean_mag"]) - num(r["phot_rp_mean_mag"]) # the single-star BP/RP excess (Riello et al. 2021) at the colour ex = (num(r["phot_bp_rp_excess_factor"]) or 0) - (1.154360 + 0.033772 * c + 0.032277 * c * c) o = dict(id=str(r["source_id"]), c=c, g=num(r["phot_g_mean_mag"]), plx=num(r["parallax"]), ex=ex, ruwe=num(r["ruwe"]), wd=num(r["wd_prob"]), multi=num(r["ipd_frac_multi_peak"]), simbad=sim.get(str(r["source_id"]))) typ = (o["simbad"] or (None, "", ""))[1] o["vet"] = ("artefact" if abs(ex) > excess_max or (o["multi"] or 0) > 2 or (o["ruwe"] or 0) > 1.4 else "white-dwarf" if (o["wd"] or 0) > 0.5 or typ in ("WD*", "WD?", "HS*", "HS?", "sdB", "sdO") or "D" in (o["simbad"] or (None, "", ""))[2][:2] else "natural-type" if typ not in ("", "*", "PM*", "HV*") else None) out.append(o) left = [o for o in out if o["vet"] is None] print("%d stars within 100 pc with BP-RP below %g: %d with a BP/RP flux excess beyond %g, several peaks or RUWE above 1.4, %d white dwarfs or hot subdwarfs, %d of other " "known types; unexplained %d" % (len(out), colour_max, sum(o["vet"] == "artefact" for o in out), excess_max, sum(o["vet"] == "white-dwarf" for o in out), sum(o["vet"] == "natural-type" for o in out), len(left))) for o in sorted(out, key=lambda o: o["c"]): print(" Gaia DR3 %s: BP-RP %.3f, G %.2f, parallax %.1f mas, excess %+.3f, RUWE %s, WD probability %s; SIMBAD %s; %s" % ( o["id"], o["c"], o["g"], o["plx"], o["ex"], o["ruwe"], o["wd"], o["simbad"], o["vet"] or "unexplained")) json.dump(out, open(core.results("bluest_nearby.json"), "w"), indent=1) if record: core.record_run( "gcns_bluer_than_stars", "Nearby light bluer than any star", "Gaia Catalogue of Nearby Stars (external.gaiaedr3_gcns_main_1); SIMBAD", "stars within 100 pc with BP-RP below %g: %d" % (colour_max, len(out)), os.path.abspath(__file__), "bluest", len(out), "sources", dict(colour_max=colour_max, excess_max=excess_max), "the BP/RP flux excess, the image-peak and RUWE flags, the white-dwarf probability and SIMBAD's types", "a source bluer than the Rayleigh-Jeans limit of any photosphere in Gaia's bands", "%d: %d flux-excess or blended artefacts, %d white dwarfs or hot subdwarfs, %d other known types; %d unexplained" % ( len(out), sum(o["vet"] == "artefact" for o in out), sum(o["vet"] == "white-dwarf" for o in out), sum(o["vet"] == "natural-type" for o in out), len(left)), [dict(target="gaia_dr3:%s" % o["id"], stat=o["c"], unit="BP-RP", vet=o["vet"], note="G %.2f, parallax %.1f mas, BP/RP excess %+.3f, RUWE %s, WD probability %s; SIMBAD %s" % (o["g"], o["plx"], o["ex"], o["ruwe"], o["wd"], o["simbad"])) for o in out]) print("recorded") APPROACHES = { "bluest": (lambda: core.run_cli(bluest_main), bluest_DOC), "etzuv": (lambda: core.run_cli(etzuv_main), etzuv_DOC), "starcat": (starcat_cli, starcat_DOC), "cos": (cos_cli, cos_DOC), "iue": (iue_cli, iue_DOC), "fuse": (fuse_cli, fuse_DOC), "stis": (stis_cli, stis_DOC), "limits": (limits_cli, limits_DOC), "hst_pulses": (hst_pulses_cli, hst_pulses_DOC), "galex_pulses": (galex_pulses_cli, galex_pulses_DOC), "galex_census": (galex_census_cli, galex_census_DOC), }