#!/usr/bin/env python3 """Signal model: WISE synthetic photometry of (i) an old stellar population and (ii) a blackbody that re-radiates a fraction of the stellar bolometric luminosity, both at redshift z. Outputs (data/): model_colour_corrections.csv my synthetic WISE colour corrections f_c vs the published values (validation) model_templates.csv per Brown+14 elliptical: L_bol/(nu L_nu)(W1), synthetic W1-W2, W1-W3, W1-W4 at z=0 model_grid.npz R_b(T, z) for b = W1..W4: catalogue-equivalent flux density of a blackbody of luminosity L in band b, divided by the catalogue W1 flux density of the stellar template with bolometric luminosity L, both at redshift z; and the stellar colours C_b(z) = F_b/F_1 of the template (b = 2, 3, 4) Conventions. WISE catalogue magnitudes are converted to "catalogue flux densities" with the zero-magnitude flux densities for a nu^-2 reference spectrum (Jarrett+ 2011): F_cat = F0 * 10^(-0.4 m). For a source with spectrum S_nu the same quantity is F_cat = INT S_nu R dnu/nu / INT (nu_iso/nu)^2 R dnu/nu with R the relative response per photon (Wright+ 2010). This is exactly f_c * S_nu(nu_iso) in the notation of the WISE Explanatory Supplement, so no separate colour correction is needed anywhere in the analysis. """ import glob import os import numpy as np import pandas as pd import mirq_common as M H_K = 4.799243e-11 # h/k in K s C_UM = 2.99792458e14 # micron / s SIGMA_SB = 5.670374e-5 # erg cm^-2 s^-1 K^-4 H_CGS = 6.62607015e-27 C_CGS = 2.99792458e10 # Ellipticals of the Brown+14 atlas without AGN / star-formation classification (table1 Class '---', T <= -3), # excluding NGC 4486 (M87: synchrotron-dominated mid-IR) and NGC 0750 (interacting pair). PASSIVE_E = ["NGC_0584", "NGC_3379", "NGC_4125", "NGC_4168", "NGC_4365", "NGC_4387", "NGC_4458", "NGC_4473", "NGC_4551", "NGC_4552", "NGC_4621", "NGC_4660", "NGC_4860", "NGC_4889", "IC_4051", "NGC_4926"] # Published colour corrections f_c (Wright et al. 2010, Table 1 / WISE Explanatory Supplement IV.4.h Table 2). # The two power-law rows are quoted in the text of the Explanatory Supplement; the blackbody rows are from the # published table (entered by hand). FC_PUBLISHED = { "nu^0": (0.9907, 0.9935, 0.9169, 0.9905), "nu^-1": (0.9921, 0.9943, 0.9373, 0.9926), "B(100K)": (17.2062, 3.9096, 2.6588, 1.0032), "B(141K)": (4.0882, 1.9739, 1.4002, 0.9852), "B(200K)": (2.0577, 1.3448, 1.0006, 0.9833), "B(283K)": (1.3917, 1.1124, 0.8791, 0.9865), "B(400K)": (1.1316, 1.0229, 0.8622, 0.9903), "B(566K)": (1.0263, 0.9919, 0.8833, 0.9935), "B(800K)": (0.9884, 0.9853, 0.9125, 0.9958), } _RSR = None def load_rsr(): global _RSR if _RSR is None: _RSR = {} for b in (1, 2, 3, 4): a = np.genfromtxt(os.path.join(M.RAW, "wise_rsr", "RSR-W%d.txt" % b), comments="#", usecols=(0, 1)) _RSR[b] = (a[:, 0], a[:, 1]) # full tabulated curve, wings included return _RSR def planck_nu(lam_um, T): """B_nu (erg s^-1 cm^-2 Hz^-1 sr^-1) at wavelength lam_um (micron).""" nu = C_UM / np.asarray(lam_um, float) x = H_K * nu / T with np.errstate(over="ignore"): return 2 * H_CGS * nu ** 3 / C_CGS ** 2 / np.expm1(x) def synth(fnu_of_lam, band, z=0.0): """Catalogue-equivalent flux density in a WISE band of a source whose *rest-frame* S_nu(lambda_rest) is given by the callable fnu_of_lam, observed at redshift z (the common (1+z)/(4 pi D_L^2) factor is omitted).""" lam, r = load_rsr()[band] nu = C_UM / lam s = fnu_of_lam(lam / (1.0 + z)) ref = (C_UM / M.WISE_LAM_ISO[band] / nu) ** 2 # integrate over ln(nu): dnu/nu. lam grid is uniform, so use trapz in ln(lam) (|dln nu| = dln lam) x = np.log(lam) return np.trapezoid(s * r, x) / np.trapezoid(ref * r, x) def bb_cat_flux(T, band, z=0.0): """Catalogue W-band flux density per unit bolometric luminosity of a blackbody: L_nu/L = pi B_nu / (sigma T^4).""" return synth(lambda lr: np.pi * planck_nu(lr, T) / (SIGMA_SB * T ** 4), band, z) def validate_colour_corrections(): rows = [] for name, pub in FC_PUBLISHED.items(): for b in (1, 2, 3, 4): if name.startswith("nu"): p = float(name.split("^")[1]) f = lambda lr, p=p: (C_UM / lr) ** p # noqa: E731 else: T = float(name[2:-2]) f = lambda lr, T=T: planck_nu(lr, T) # noqa: E731 fc = synth(f, b) / f(M.WISE_LAM_ISO[b]) rows.append(dict(spectrum=name, band="W%d" % b, fc_this_work=fc, fc_published=pub[b - 1], ratio=fc / pub[b - 1])) out = pd.DataFrame(rows) out.to_csv(os.path.join(M.DATA, "model_colour_corrections.csv"), index=False, float_format="%.5g") w34 = out[out.band.isin(["W3", "W4"])] print("colour-correction check: W3/W4 max |ratio-1| = %.3f ; all bands, T>=200K: %.3f" % (np.abs(w34.ratio - 1).max(), np.abs(out[~out.spectrum.isin(["B(100K)", "B(141K)"])].ratio - 1).max())) return out def load_template(name): a = np.loadtxt(os.path.join(M.RAW, "brown14", "sp", name + "_spec.dat"), comments="#", usecols=(0, 1)) lam = a[:, 0] * 1e-4 # micron (rest frame) flam = a[:, 1] # erg/s/cm^2/A o = np.argsort(lam) lam, flam = lam[o], flam[o] u, iu = np.unique(lam, return_index=True) return u, flam[iu] def template_props(name): lam, flam = load_template(name) fnu = flam * (lam * 1e4) ** 2 / 2.99792458e18 # erg/s/cm^2/Hz # stellar bolometric flux: integrate F_lambda from 0.0912 to 5 micron on the data, plus a Rayleigh-Jeans # tail beyond 5 micron anchored on the 3-5 micron continuum (dust emission at >5 micron is NOT starlight) m = (lam >= 0.0912) & (lam <= 5.0) fbol_data = np.trapezoid(flam[m], lam[m] * 1e4) k = (lam > 3.0) & (lam < 5.0) a5 = np.median(flam[k] * lam[k] ** 4) # F_lambda ~ a5 lam^-4 fbol_tail = a5 / 3.0 / 5.0 ** 3 * 1e4 fbol = fbol_data + fbol_tail f_of = lambda lr: np.interp(lr, lam, fnu, left=0.0, right=np.nan) # noqa: E731 # extend the template beyond its red end with nu^2 (only matters for W4 at z~0 if the spectrum stops early) lmax = lam.max() fmax = np.median(fnu[lam > 0.95 * lmax]) def f_ext(lr): lr = np.asarray(lr, float) out = np.interp(lr, lam, fnu, left=0.0) return np.where(lr > lmax, fmax * (lmax / lr) ** 2, out) cat = {b: synth(f_ext, b) for b in (1, 2, 3, 4)} nu1 = C_UM / M.WISE_LAM_ISO[1] row = dict(name=name, lam_max_um=lmax, frac_bol_beyond_5um=fbol_tail / fbol, Lbol_over_nuLnu_W1=fbol / (nu1 * cat[1])) for b in (2, 3, 4): row["W1-W%d" % b] = 2.5 * np.log10(cat[b] / cat[1] * M.WISE_F0[1] / M.WISE_F0[b]) return row, f_ext, fbol def build(zmax=0.14, nz=57): validate_colour_corrections() zs = np.linspace(0.0, zmax, nz) Ts = np.unique(np.concatenate([M.T_GRID, [125.0, 250.0, 500.0, 800.0, 1000.0]])) rows, star = [], [] for name in PASSIVE_E: if not os.path.exists(os.path.join(M.RAW, "brown14", "sp", name + "_spec.dat")): continue row, f_ext, fbol = template_props(name) rows.append(row) # catalogue fluxes of the stellar template per unit stellar bolometric flux, vs z star.append(np.array([[synth(f_ext, b, z) / fbol for z in zs] for b in (1, 2, 3, 4)])) tp = pd.DataFrame(rows) tp.to_csv(os.path.join(M.DATA, "model_templates.csv"), index=False, float_format="%.5g") star = np.array(star) # (ntemplate, 4, nz) print(tp.to_string(index=False, float_format=lambda v: "%.3f" % v)) kb = tp.Lbol_over_nuLnu_W1 print("L_bol/(nu L_nu)(W1): median %.2f, min %.2f, max %.2f, rms %.2f" % (kb.median(), kb.min(), kb.max(), kb.std())) star_med = np.median(star, axis=0) # (4, nz) per unit L_bol bb = np.array([[[bb_cat_flux(T, b, z) for z in zs] for b in (1, 2, 3, 4)] for T in Ts]) # (nT, 4, nz) R = bb / star_med[0][None, None, :] # blackbody band flux / stellar W1 flux at equal L R_lo = bb / star.max(axis=0)[0][None, None, :] # template with the largest W1 per L_bol -> smallest R R_hi = bb / star.min(axis=0)[0][None, None, :] C = star_med[1:] / star_med[0][None, :] # stellar colours F_b/F_1, b=2,3,4 C_all = star[:, 1:, :] / star[:, :1, :] np.savez(os.path.join(M.DATA, "model_grid.npz"), z=zs, T=Ts, R=R, R_lo=R_lo, R_hi=R_hi, C=C, C_min=C_all.min(axis=0), C_max=C_all.max(axis=0), kbol=kb.values) return tp class Model: """Interpolator for R_b(T, z) and stellar colours C_b(z).""" def __init__(self, which="R"): g = np.load(os.path.join(M.DATA, "model_grid.npz")) self.z, self.T, self.R, self.C = g["z"], g["T"], g[which], g["C"] self.C_min = g["C_min"] def R_of(self, T, band, z): i = int(np.where(np.isclose(self.T, T))[0][0]) return np.interp(z, self.z, self.R[i, band - 1]) def C_of(self, band, z, lower=False): return np.interp(z, self.z, (self.C_min if lower else self.C)[band - 2]) def main(): build() if __name__ == "__main__": main()