diff --git a/docs/references/Behroozi+2013_ApJ_770_57.pdf b/docs/references/Behroozi+2013_ApJ_770_57.pdf new file mode 100644 index 00000000..71d73d91 Binary files /dev/null and b/docs/references/Behroozi+2013_ApJ_770_57.pdf differ diff --git a/docs/references/Finn+Thorne-2000.pdf b/docs/references/Finn+Thorne-2000.pdf new file mode 100644 index 00000000..618ba813 Binary files /dev/null and b/docs/references/Finn+Thorne-2000.pdf differ diff --git a/docs/references/Guo+2010_mnras0404-1111.pdf b/docs/references/Guo+2010_mnras0404-1111.pdf new file mode 100644 index 00000000..2bdbdc7a Binary files /dev/null and b/docs/references/Guo+2010_mnras0404-1111.pdf differ diff --git a/docs/references/Leja+2020_ApJ_893_111.pdf b/docs/references/Leja+2020_ApJ_893_111.pdf new file mode 100644 index 00000000..cb3ee7eb Binary files /dev/null and b/docs/references/Leja+2020_ApJ_893_111.pdf differ diff --git a/docs/references/Leja_2020_ApJ_893_111.pdf b/docs/references/Leja_2020_ApJ_893_111.pdf new file mode 100644 index 00000000..cb3ee7eb Binary files /dev/null and b/docs/references/Leja_2020_ApJ_893_111.pdf differ diff --git a/docs/references/MacLeod+2010_1004.0276.pdf b/docs/references/MacLeod+2010_1004.0276.pdf new file mode 100644 index 00000000..206371ba Binary files /dev/null and b/docs/references/MacLeod+2010_1004.0276.pdf differ diff --git a/docs/references/Moore+Cole+Berry-2015_1408.0740.pdf b/docs/references/Moore+Cole+Berry-2015_1408.0740.pdf new file mode 100644 index 00000000..e5afe86a Binary files /dev/null and b/docs/references/Moore+Cole+Berry-2015_1408.0740.pdf differ diff --git a/docs/references/Navarro+1997_ApJ_490_493.pdf b/docs/references/Navarro+1997_ApJ_490_493.pdf new file mode 100644 index 00000000..38f90a46 Binary files /dev/null and b/docs/references/Navarro+1997_ApJ_490_493.pdf differ diff --git a/docs/references/Rodriguez-Gomez+2015_2015MNRAS.449...49R.pdf b/docs/references/Rodriguez-Gomez+2015_2015MNRAS.449...49R.pdf new file mode 100644 index 00000000..20cdc5f0 Binary files /dev/null and b/docs/references/Rodriguez-Gomez+2015_2015MNRAS.449...49R.pdf differ diff --git a/docs/references/Romano+Cornish-2017_1608.06889.pdf b/docs/references/Romano+Cornish-2017_1608.06889.pdf new file mode 100644 index 00000000..4facfc91 Binary files /dev/null and b/docs/references/Romano+Cornish-2017_1608.06889.pdf differ diff --git a/docs/references/Runnoe+2012_1201.5155.pdf b/docs/references/Runnoe+2012_1201.5155.pdf new file mode 100644 index 00000000..37c6884f Binary files /dev/null and b/docs/references/Runnoe+2012_1201.5155.pdf differ diff --git a/docs/references/Sesana+Vecchio+Colacino-2008.pdf b/docs/references/Sesana+Vecchio+Colacino-2008.pdf new file mode 100644 index 00000000..0205e6b9 Binary files /dev/null and b/docs/references/Sesana+Vecchio+Colacino-2008.pdf differ diff --git a/docs/references/Tomczak_2014_ApJ_783_85.pdf b/docs/references/Tomczak_2014_ApJ_783_85.pdf new file mode 100644 index 00000000..bbd595d8 Binary files /dev/null and b/docs/references/Tomczak_2014_ApJ_783_85.pdf differ diff --git a/holodeck/ems/__init__.py b/holodeck/ems/__init__.py new file mode 100644 index 00000000..62ef9099 --- /dev/null +++ b/holodeck/ems/__init__.py @@ -0,0 +1,13 @@ +"""Holodeck - Electromagnetic Signals (EMS) +""" + +from . import basics # noqa +from . import drw # noqa +from . runnoe2012 import Runnoe2012 # noqa + + +bands_sdss = basics.Bands_SDSS() +bands_lsst = basics.Bands_LSST() +MacLeod2010 = drw.MacLeod2010 +runnoe2012 = Runnoe2012() + diff --git a/holodeck/ems/basics.py b/holodeck/ems/basics.py new file mode 100644 index 00000000..0afa9b23 --- /dev/null +++ b/holodeck/ems/basics.py @@ -0,0 +1,398 @@ +"""Astronomical observations calculations. + +http://svo2.cab.inta-csic.es/theory/fps/index.php?mode=browse&gname=SLOAN&asttype= +See also: https://ui.adsabs.harvard.edu/abs/2018ApJS..236...47W/abstract + +- lots of reference fluxes: https://coolwiki.ipac.caltech.edu/index.php/Central_wavelengths_and_zero_points +- VEGA/Johnson/Bessell: http://web.ipac.caltech.edu/staff/fmasci/home/astro_refs/magsystems.pdf +- SDSS/AB/Fukugita: http://www.astronomy.ohio-state.edu/~martini/usefuldata.html + + + "u": { + "wlen": * ap.units.angstrom, # wavelength effective [Angstrom] + "bandwidth_wlen": * ap.units.angstrom, # bandwidth effective [Angstrom] + "AB": { + "flux_ref_wlen": * _units_erg_s_cm2_angstrom, + "flux_ref_freq": * ap.units.jansky, + }, + "vega": { + "flux_ref_wlen": * _units_erg_s_cm2_angstrom, + "flux_ref_freq": * ap.units.jansky, + }, + }, + +""" + +import numpy as np +import astropy as ap + +UNITS_WLEN = 'angstrom' +UNITS_FREQ = 'Hz' + +UNITS_FLUX_WLEN = 'erg/(s cm2 angstrom)' +UNITS_FLUX_FREQ = 'erg/(s cm2 Hz)' + +_units_erg_s_cm2_angstrom = ap.units.erg / ap.units.second / ap.units.cm**2 / ap.units.angstrom + +# SDSS AB Magnitudes +# from http://svo2.cab.inta-csic.es/theory/fps/index.php?mode=browse&gname=SLOAN&asttype= +# Using "full transmission" filters, obtained on 2023-09-06 +BANDS_SDSS = { + "u": { + "wlen": 3608.04 * ap.units.angstrom, + "bandwidth_wlen": 540.97 * ap.units.angstrom, + "vega": { + "flux_ref_freq": 1582.54 * ap.units.jansky, + "flux_ref_wlen": 3.75079e-9 * _units_erg_s_cm2_angstrom, + }, + "AB": { + "flux_ref_freq": 3631.00 * ap.units.jansky, + "flux_ref_wlen": 8.60588e-9 * _units_erg_s_cm2_angstrom, + }, + }, + "g": { + "wlen": 4671.78 * ap.units.angstrom, + "bandwidth_wlen": 1064.68 * ap.units.angstrom, + "vega": { + "flux_ref_freq": 4023.57 * ap.units.jansky, + "flux_ref_wlen": 5.45476e-9 * _units_erg_s_cm2_angstrom, + }, + "AB": { + "flux_ref_freq": 3631.00 * ap.units.jansky, + "flux_ref_wlen": 4.92255e-9 * _units_erg_s_cm2_angstrom, + }, + }, + "r": { + "wlen": 6141.12 * ap.units.angstrom, + "bandwidth_wlen": 1055.51 * ap.units.angstrom, + "vega": { + "flux_ref_freq": 3177.38 * ap.units.jansky, + "flux_ref_wlen": 2.49767e-9 * _units_erg_s_cm2_angstrom, + }, + "AB": { + "flux_ref_freq": 3631.00 * ap.units.jansky, + "flux_ref_wlen": 2.85425e-9 * _units_erg_s_cm2_angstrom, + }, + }, + "i": { + "wlen": 7457.89 * ap.units.angstrom, + "bandwidth_wlen": 1102.57 * ap.units.angstrom, + "vega": { + "flux_ref_freq": 2593.40 * ap.units.jansky, + "flux_ref_wlen": 1.38589e-9 * _units_erg_s_cm2_angstrom, + }, + "AB": { + "flux_ref_freq": 3631.00 * ap.units.jansky, + "flux_ref_wlen": 1.94038e-9 * _units_erg_s_cm2_angstrom, + }, + }, + "z": { + "wlen": 8922.78 * ap.units.angstrom, + "bandwidth_wlen": 1164.01 * ap.units.angstrom, + "vega": { + "flux_ref_freq": 2238.99 * ap.units.jansky, + "flux_ref_wlen": 8.38585e-10 * _units_erg_s_cm2_angstrom, + }, + "AB": { + "flux_ref_freq": 3631.00 * ap.units.jansky, + "flux_ref_wlen": 1.35994e-9 * _units_erg_s_cm2_angstrom, + }, + }, +} + +BANDS_LSST = { + "u": { + "wlen": 3751.20 * ap.units.angstrom, + "bandwidth_wlen": 473.19 * ap.units.angstrom, + "AB": { + "flux_ref_wlen": 8.03787e-9 * _units_erg_s_cm2_angstrom, + "flux_ref_freq": 3631.00 * ap.units.jansky, + }, + }, + "g": { + "wlen": 4740.66 * ap.units.angstrom, + "bandwidth_wlen": 1253.26 * ap.units.angstrom, + "AB": { + "flux_ref_wlen": 4.7597e-9 * _units_erg_s_cm2_angstrom, + "flux_ref_freq": 3631.00 * ap.units.jansky, + }, + }, + "r": { + "wlen": 6172.34 * ap.units.angstrom, + "bandwidth_wlen": 1206.92 * ap.units.angstrom, + "AB": { + "flux_ref_wlen": 2.8156e-9 * _units_erg_s_cm2_angstrom, + "flux_ref_freq": 3631.00 * ap.units.jansky, + }, + }, + "i": { + "wlen": 7500.97 * ap.units.angstrom, + "bandwidth_wlen": 1174.77 * ap.units.angstrom, + "AB": { + "flux_ref_wlen": 1.91864e-9 * _units_erg_s_cm2_angstrom, + "flux_ref_freq": 3631.00 * ap.units.jansky, + }, + }, + "z": { + "wlen": 8678.90 * ap.units.angstrom, + "bandwidth_wlen": 997.51 * ap.units.angstrom, + "AB": { + "flux_ref_wlen": 1.44312e-9 * _units_erg_s_cm2_angstrom, + "flux_ref_freq": 3631.00 * ap.units.jansky, + }, + }, + "y": { + "wlen": 9711.82 * ap.units.angstrom, + "bandwidth_wlen": 871.83 * ap.units.angstrom, + "AB": { + "flux_ref_wlen": 1.14978e-9 * _units_erg_s_cm2_angstrom, + "flux_ref_freq": 3631.00 * ap.units.jansky, + }, + }, +} + + +def _get_wlen_freq(wlen, freq, error_if_neither): + if error_if_neither and (wlen is None) and (freq is None): + raise ValueError("neither 'wlen' or 'freq' is given!") + + if (wlen is None) and (freq is not None): + freq = ap.units.Quantity(freq, UNITS_FREQ) + wlen = Band.freq_to_wlen(freq) + if (freq is None) and (wlen is not None): + wlen = ap.units.Quantity(wlen, UNITS_WLEN) + freq = Band.wlen_to_freq(wlen) + + return wlen, freq + + +def _get_flux_wlen_flux_freq(flux_wlen, flux_freq, freq, wlen): + if (flux_wlen is None) and (flux_freq is None): + raise ValueError("neither 'flux_wlen' or 'flux_freq' has been specified!") + + if (flux_wlen is None): + flux_freq = ap.units.Quantity(flux_freq, UNITS_FLUX_FREQ) + flux_wlen = Band.spectral_wlen(flux_freq, freq=freq) + if (flux_freq is None): + flux_wlen = ap.units.Quantity(flux_wlen, UNITS_FLUX_WLEN) + flux_freq = Band.spectral_freq(flux_wlen, wlen=wlen) + + return flux_wlen, flux_freq + + +class Band: + + def __init__(self, name, wlen, freq, ref_flux_wlen, flux_ref_freq, + bandwidth_wlen=None, bandwidth_freq=None): + wlen, freq = _get_wlen_freq(wlen, freq, error_if_neither=True) + ref_flux_wlen, flux_ref_freq = _get_flux_wlen_flux_freq(ref_flux_wlen, flux_ref_freq, freq, wlen) + bandwidth_wlen, bandwidth_freq = _get_wlen_freq(bandwidth_wlen, bandwidth_freq, error_if_neither=False) + + self.name = name + self.wlen = wlen + self.freq = freq + self.flux_ref_wlen = ap.units.Quantity(ref_flux_wlen, UNITS_FLUX_WLEN) + self.flux_ref_freq = ap.units.Quantity(flux_ref_freq, UNITS_FLUX_FREQ) + self.bandwidth_freq = bandwidth_freq + self.bandwidth_freq = bandwidth_wlen + return + + def __str__(self): + rv = ( + f"{self.name}: w={self.wlen:.4e}, f={self.freq:.4e} | " + f"F_w={self.flux_ref_wlen:.4e}, F_f={self.flux_ref_freq:.4e}" + ) + return rv + + @classmethod + def wlen_to_freq(cls, wlen): + return ap.units.Quantity(ap.constants.c / wlen, UNITS_FREQ) + + @classmethod + def freq_to_wlen(cls, freq): + return ap.units.Quantity(ap.constants.c / freq, UNITS_WLEN) + + @classmethod + def spectral_freq(cls, spec_wlen, wlen): + spec_freq = ap.units.Quantity(spec_wlen * wlen**2 / ap.constants.c, UNITS_FLUX_FREQ) + return spec_freq + + @classmethod + def spectral_wlen(cls, spec_freq, freq): + spec_wlen = ap.units.Quantity(spec_freq * freq**2 / ap.constants.c, UNITS_FLUX_WLEN) + return spec_wlen + + def mag_to_flux(self, mag, type): + """Convert from broad-band filter magnitude to spectral flux. + + Returns + ------- + flux : () scalar + Flux in either [erg/s/cm^2/Hz] or [erg/s/cm^2/Angstrom] depending on `type`. + + """ + ref_flux = self._ref_flux_for_type(type) + flux = ref_flux * np.power(10.0, mag/-2.5) + return flux + + def flux_to_mag(self, flux, type, units=None): + """Convert from spectral flux to broad-band filter magnitude. + + Arguments + --------- + flux : () scalar + Flux in either [erg/s/cm^2/Hz] or [erg/s/cm^2/Angstrom] depending on `type`. + type + + Returns + ------- + mag + + """ + if units is not None: + flux = ap.units.Quantity(flux, units) + + ref_flux = self._ref_flux_for_type(type) + mag = flux / ref_flux + + try: + mag = mag.to('') + except ap.units.UnitConversionError as err: + msg = ( + "Could not convert 'flux' to a spectral flux. " + f"Try using the `units` argument, e.g. ``units='erg/s/cm2/Hz'`` or ``units='erg/s/cm2/Angstrom'``. " + f"({err})" + ) + raise ValueError(msg) + + mag = -2.5 * np.log10(mag) + return mag + + def abs_mag_to_lum(self, abs_mag, type): + """Convert from broad-band filter absolute-magnitude to spectral luminosity. + """ + ref_flux = self._ref_flux_for_type(type) + lum = 4.0 * np.pi * ref_flux * (10.0 * ap.units.parsec)**2 * np.power(10.0, -abs_mag/2.5) + return lum + + def lum_to_abs_mag(self, lum, type, units=None): + if units is not None: + lum = ap.units.Quantity(lum, units) + + ref_flux = self._ref_flux_for_type(type) + mag = lum / (ref_flux * 4.0 * np.pi * (10.0 * ap.units.parsec)**2) + + try: + mag = mag.to('') + except ap.units.UnitConversionError as err: + msg = ( + "Could not convert 'lum' to a spectral luminosity. " + f"Try using the `units` argument, e.g. ``units='erg/s/Hz'`` or ``units='erg/s/Angstrom'``. " + f"({err})" + ) + raise ValueError(msg) + + mag = -2.5 * np.log10(mag) + return mag + + def _ref_flux_for_type(self, type): + """Get the appropriate reference flux for the given 'type'. + + If `type` is '[f]requency' : return `flux_ref_freq` + If `type` is '[w]avelength' : return `flux_ref_wlen` + + Arguments + --------- + type : string, + Specification of wavelength or frequency. + + Returns + ------- + ref_flux : astropy.units.Quantity, + Reference flux for this band, either F_nu or F_lambda based on the `type` argument. + If `type` == '[f]requency', then F_nu is returned (e.g. in units of erg/s/cm^2/Hz) + If `type` == '[w]avelength', then F_lambda is returned (e.g. in units of erg/s/cm^2/Angstrom) + + """ + if type.startswith('f'): + ref_flux = self.flux_ref_freq + elif type.startswith('w'): + ref_flux = self.flux_ref_wlen + else: + raise ValueError(f"`type` ({type}) should be '[f]requency' or '[w]avelength'!") + return ref_flux + + +class BANDS: + + def __init__(self, bands_dict, mag_type): + bands = {} + for name, values in bands_dict.items(): + + wlen = values.get('wlen', None) + freq = values.get('freq', None) + if (wlen is None) and (freq is None): + raise ValueError(f"Band {name} has neither 'wlen' or 'freq' specification!") + if wlen is None: + wlen = Band.freq_to_wlen(freq) + if freq is None: + freq = Band.wlen_to_freq(wlen) + + bandwidth_wlen = values.get('bandwidth_wlen', None) + bandwidth_freq = values.get('bandwidth_freq', None) + if (bandwidth_wlen is None) and (bandwidth_freq is not None): + bandwidth_wlen = Band.freq_to_wlen(bandwidth_freq) + if (bandwidth_freq is None) and (bandwidth_wlen is not None): + bandwidth_freq = Band.wlen_to_freq(bandwidth_wlen) + + zero_points = values.get(mag_type, None) + if zero_points is None: + raise ValueError(f"Band '{name}' does not have specification for mag_type '{mag_type}'!") + flux_ref_wlen = zero_points.get('flux_ref_wlen', None) + flux_ref_freq = zero_points.get('flux_ref_freq', None) + if (flux_ref_wlen is None) and (flux_ref_freq is None): + raise ValueError(f"Band '{name}' '{mag_type}' has neither 'flux_ref_wlen' nor 'flux_ref_freq'!") + if (flux_ref_wlen is None): + flux_ref_wlen = Band.spectral_wlen(flux_ref_freq, freq=freq) + if (flux_ref_freq is None): + flux_ref_freq = Band.spectral_freq(flux_ref_wlen, wlen=wlen) + + band = Band( + name, wlen, freq, flux_ref_wlen, flux_ref_freq, + bandwidth_wlen=bandwidth_wlen, bandwidth_freq=bandwidth_freq + ) + bands[name] = band + + if len(bands) == 0: + raise ValueError("No bands provided!") + + self._bands = bands + return + + def __getitem__(self, name): + return self._bands[name] + + def __call__(self, name): + return self._bands[name] + + @property + def names(self): + return self._bands.keys() + + +class Bands_SDSS(BANDS): + """SDSS Generally uses AB magnitudes. + """ + + def __init__(self): + super().__init__(BANDS_SDSS, "AB") + return + + +class Bands_LSST(BANDS): + """LSST Generally uses AB magnitudes. + """ + + def __init__(self): + super().__init__(BANDS_LSST, "AB") + return \ No newline at end of file diff --git a/holodeck/ems/drw.py b/holodeck/ems/drw.py new file mode 100644 index 00000000..a0673454 --- /dev/null +++ b/holodeck/ems/drw.py @@ -0,0 +1,140 @@ +"""Damped Random Walks +""" + +import numpy as np + +import holodeck as holo +from holodeck.constants import MSOL +# from holodeck.ems import runnoe2012, bands_sdss + + +class MacLeod2010: + """Damped Random Walk models from MacLeod+2010. + + MacLeod et al. 2010 - Modeling the Time Variability of SDSS Stripe 82 Quasars as a Damped Random Walk + https://arxiv.org/abs/1004.0276 + https://ui.adsabs.harvard.edu/abs/2010ApJ...721.1014M/abstract + """ + + @classmethod + def _fit_func(cls, pars, errs, lambda_rf, imag, mbh, randomize=False): + if (randomize is not None) and (randomize is not False): + shape = np.shape(mbh) + if int(randomize) > 1: + shape = shape + (int(randomize),) + if not np.isscalar(imag): + imag = imag[..., np.newaxis] + if not np.isscalar(mbh): + mbh = mbh[..., np.newaxis] + + shape = shape + (len(pars),) + pars = np.random.normal(pars, errs, size=shape) + pars = np.moveaxis(pars, -1, 0) + + aa, bb, cc, dd = pars + + # lf = aa + bb*np.log(lambda_rf/4000e-8) + cc*(imag + 23) + dd*np.log(mbh/(1e9*MSOL)) + # ff = np.exp(lf) + rv = aa + bb*np.log10(lambda_rf/4000e-8) + cc*(imag + 23) + dd*np.log10(mbh/(1e9*MSOL)) + rv = 10**rv + return rv + + @classmethod + def sfinf(cls, imag, mbh, randomize=False): + """`mbh` should be in grams (NOTE: I THINK!) + """ + lambda_iband = 7690e-8 + pars = [-0.51, -0.479, 0.131, 0.18] + errs = [0.02, 0.005, 0.008, 0.03] + return cls._fit_func(pars, errs, lambda_iband, imag, mbh, randomize=randomize) + + @classmethod + def tau(cls, imag, mbh, randomize=False): + """`mbh` should be in grams (NOTE: I THINK!) + `tau` is returned in units of days (NOTE: I THINK!) + """ + lambda_iband = 7690e-8 + pars = [2.4, 0.17, 0.03, 0.21] + errs = [0.2, 0.02, 0.04, 0.07] + return cls._fit_func(pars, errs, lambda_iband, imag, mbh, randomize=randomize) + + +def drw_lightcurve(times, tau, mean_mag, sfinf, size=None): + """Construct an AGN DRW lightcurve based on MacLeod+2010 model. + + Arguments + --------- + times : (N,) array_like of scalar + Times at which to sample lightcurve + tau : scalar, + correlation timescale + mean_mag : scalar, + Mean absolute magnitude of AGN. + sfinf : scalar, + Structure-Function at Infinity + (i.e. measure of variance of DRW, `sigma = sfinf / sqrt(2)`). + + Returns + ------- + mags : (N,) ndarray of scalar + Flux of continuum source in magnitudes. + lums : (N,) ndarray of scalar + Flux of continuum source in luminosity. + + """ + # sfinf = np.sqrt(2) * sigma_drw # [MacLeod+2010] Eq. 4 + num = 1 if (size is None) else size + + mags = mean_mag * np.ones((num, times.size)) + dt = np.diff(times) + + if np.isscalar(tau): + tau = np.ones(num) * tau + if np.isscalar(mean_mag): + mean_mag = np.ones(num) * mean_mag + if np.isscalar(sfinf): + sfinf = np.ones(num) * sfinf + + # [MacLeod+2010] Eq. 5 + exp = np.exp(-dt[np.newaxis, :] / tau[:, np.newaxis]) + rand = np.random.normal(size=exp.shape) + var = 0.5 * np.square(sfinf[:, np.newaxis]) * (1 - exp**2) + + for i0, (ee, vv, rr) in enumerate(zip(exp.T, var.T, rand.T)): + i1 = i0 + 1 + l0 = mags[:, i0] + mean = ee * l0 + mean_mag * (1 - ee) + + # Transform the normal random variable to mean `mean` and variance `vv` + temp = rr * np.sqrt(vv) + mean + mags[:, i1] = temp + + lums = np.power(10.0, -0.4*mags) + if size is None: + mags = np.squeeze(mags) + lums = np.squeeze(lums) + else: + mags = mags.T + lums = lums.T + + return mags, lums + + +def drw_params(mass, fedd, eps=0.1, scatter=False): + """DRW Parameters + + Returns + ------- + imag + Mean i-band absolute magnitude + tau + Correlation time of DRW + sfi + Structure-Function at Infinity + + """ + imag = holo.ems.runnoe2012.iband_from_mass_fedd(mass, fedd, eps=eps, magnitude=True).value + taus = MacLeod2010.tau(imag, mass, randomize=scatter) + sfis = MacLeod2010.sfinf(imag, mass, randomize=scatter) + return imag, taus, sfis + diff --git a/holodeck/ems/runnoe2012.py b/holodeck/ems/runnoe2012.py new file mode 100644 index 00000000..ba30f9ae --- /dev/null +++ b/holodeck/ems/runnoe2012.py @@ -0,0 +1,205 @@ +"""Functions and relations from Runnoe et al. 2012. + +Runnoe+2012 [1201.5155] - Updating quasar bolometric luminosity corrections +https://ui.adsabs.harvard.edu/abs/2012MNRAS.422..478R/abstract + +""" + +import numpy as np +import astropy as ap +import astropy.units + +import holodeck as holo +from holodeck import utils + + +__all__ = ["Runnoe2012"] + +FRAC_ISO = 0.75 + + +class Runnoe2012: + + _FITS = { + '5100': { + 'alpha': (4.89, 1.66), + 'beta': (0.91, 0.04), + 'wlen': 5100.0 * ap.units.angstrom, + }, + '3000': { + 'alpha': (1.85, 1.27), + 'beta': (0.98, 0.03), + 'wlen': 3000.0 * ap.units.angstrom, + }, + '1450': { + 'alpha': (4.74, 1.00), + 'beta': (0.91, 0.02), + 'wlen': 1450.0 * ap.units.angstrom, + }, + '2to10': { + 'alpha': (25.14, 1.93), + 'beta': (0.47, 0.043), + 'wlen': None, + }, + '2to10rl': { + 'alpha': (23.04, 3.60), + 'beta': (0.52, 0.080), + 'wlen': None, + }, + '2to10rq': { + 'alpha': (33.06, 3.17), + 'beta': (0.29, 0.072), + 'wlen': None, + }, + } + + def __init__(self): + self._names = self._FITS.keys() + return + + @property + def names(self): + return self._names + + def _fit_params_for_band(self, band): + if band not in self.names: + raise KeyError(f"`band` ({band}) must be one of {self._options}!") + + vals = self._FITS[band] + alpha = vals['alpha'] + beta = vals['beta'] + wlen = ap.units.Quantity(vals['wlen'], 'angstrom') + return alpha, beta, wlen + + def lband_from_lbol(self, band, lbol, scatter=False, fiso=FRAC_ISO): + """Convert from bolometric luminosity to luminosity in photometric band. + + Arguments + --------- + band : str + Specification of which photometric band. One of `Runnoe2012.names`: + {'5100', '3000', '1450', '2to10', '2to10rl', '2to10rq'} + lbol : array_like, units of [erg/s] + Bolometric luminosity. + + Returns + ------- + lband : array_like, + Luminosity in photometric band. + * If `band` is one of the optical bands, return spectral luminosity, + with units of [erg/s/Angstrom]. + * If `band` is one of the x-ray bands, return luminosity across the band, + with units of [erg/s]. + + """ + alpha, beta, wlen = self._fit_params_for_band(band) + if not scatter: + alpha = alpha[0] + beta = beta[0] + + lbol = ap.units.Quantity(lbol, 'erg/s') + lband = _lbol_to_lband__pow_law(lbol, alpha, beta, fiso=fiso) + # if this is one of the x-ray bands, then wlen is None, and `lband` is the luminosity across + # the band, in units of erg/s + if wlen is None: + # units = 'erg / s' + pass + # if this is one of the (near-)optical bands, then `wlen` should be wavelength in Angstroms + # and `lband` starts out as lambda * F_lambda in units of [erg/s] + # convert that to just F_lambda in units of [erg/s/Angstrom] + else: + lband = lband / wlen + # units = 'erg / (s angstrom)' + + return lband + + def lbol_from_lband(self, band_name, lband, scatter=False, fiso=FRAC_ISO): + """Convert from luminosity in photometric band to bolometric luminosity. + + Arguments + --------- + band_name : str + Specification of which photometric band. One of `Runnoe2012.names`: + {'5100', '3000', '1450', '2to10', '2to10rl', '2to10rq'} + lband : array_like, + Luminosity in photometric band. + * If `band_name` is one of the optical bands, `lband` must be spectral luminosity, + with units of [erg/s/Angstrom]. + * If `band_name` is one of the x-ray bands, `lband` must be luminosity across the band, + with units of [erg/s]. + + Returns + ------- + lbol : array_like, units of [erg/s] + Bolometric luminosity. + + """ + + alpha, beta, wlen = self._fit_params_for_band(band_name) + if not scatter: + alpha = alpha[0] + beta = beta[0] + + # For the x-ray bands, `wlen` is `None`, and `lband` must have units of erg/s + if wlen is None: + pass + # For the (near-)optical bands, `wlen` is the wavelength in Angstroms, and the relation + # requires lambda * L_lambda, with units of erg/s. + # Convert from `L_lambda` (units of [erg/s/Angstrom]) to lambda * L_lambda + else: + lband = lband * wlen + + units = 'erg / s' + lband = ap.units.Quantity(lband, units) + + lbol = _lband_to_lbol__pow_law(lband, alpha, beta, fiso=0.75) + lbol = ap.units.Quantity(lbol, units) + return lbol + + def iband_from_mass_fedd(self, mass, fedd, eps=0.1, magnitude=True): + lbol = utils.eddington_luminosity(mass, eps) * fedd + lbol = ap.units.Quantity(lbol, 'erg/s') + iband = self.lband_from_lbol('5100', lbol) + if magnitude: + iband = holo.ems.bands_sdss['i'].lum_to_abs_mag(iband, type='w') + return iband + + +def _dist_pars(arg, num): + """If `arg` is tuple (2,), draw `num` points from normal distribution with those parameters. + """ + if isinstance(arg, tuple): + if len(arg) != 2: + raise ValueError("`arg` must be a tuple of (mean, std)!") + arg = np.random.normal(*arg, size=num) + + return arg + + +def _lband_to_lbol__pow_law(lam_lum_lam, alpha, beta, fiso=1.0): + """ + log(L_iso) = alpha + beta * log10(lambda * L_lambda) + L_bol = fiso * L_iso + """ + liso = alpha + beta*np.log10(lam_lum_lam.to('erg/s').value) + lbol = fiso * (10**liso) + lbol = ap.units.Quantity(lbol, 'erg/s') + return lbol + + +def _lbol_to_lband__pow_law(lbol, alpha, beta, fiso=1.0): + """Returns lambda*L_lambda + + log(L_iso) = alpha + beta * log10(lambda * L_lambda) + L_iso = L_bol / fiso + """ + liso_log = np.log10(lbol.to('erg/s').value/fiso) + num = np.size(lbol) + + alpha = _dist_pars(alpha, num) + beta = _dist_pars(beta, num) + + lband = np.power(10, (liso_log - alpha)/beta) + lband = ap.units.Quantity(lband, 'erg/s') + return lband + diff --git a/holodeck/ems/tests/basics_sdss_bands_2023-09-05.json b/holodeck/ems/tests/basics_sdss_bands_2023-09-05.json new file mode 100644 index 00000000..09408c3a --- /dev/null +++ b/holodeck/ems/tests/basics_sdss_bands_2023-09-05.json @@ -0,0 +1,97 @@ +{ + "g": { + "flux": [ + 1.6607997049926968e-31, + 1.0560362777516895e-32, + 3.7361459583907544e-33 + ], + "mag_angstrom": [ + 56.122266537265055, + 59.113862590665974, + 60.24200009911879 + ], + "mag_hz": [ + 28.34927247472027, + 31.34086852812119, + 32.469006036574 + ], + "units_angstrom": "erg/(s cm2 Angstrom)", + "units_hz": "erg/(s cm2 Hz)" + }, + "i": { + "flux": [ + 1.446549480694367e-36, + 7.605178897041373e-35, + 1.9644890737814164e-35 + ], + "mag_angstrom": [ + 67.76826922155983, + 63.46632886975941, + 64.93597841169418 + ], + "mag_hz": [ + 40.99923238797727, + 36.697292036176854, + 38.16694157811162 + ], + "units_angstrom": "erg/(s cm2 Angstrom)", + "units_hz": "erg/(s cm2 Hz)" + }, + "r": { + "flux": [ + 2.8803660328108007e-32, + 5.6190278092283624e-33, + 8.583359922817238e-36 + ], + "mag_angstrom": [ + 57.461492787819935, + 59.235959036108696, + 66.27596867975863 + ], + "mag_hz": [ + 30.251446420306976, + 32.02591266859573, + 39.06592231224567 + ], + "units_angstrom": "erg/(s cm2 Angstrom)", + "units_hz": "erg/(s cm2 Hz)" + }, + "u": { + "flux": [ + 1.9022465038219896e-35, + 9.130625455831697e-31, + 2.4280002928166076e-36 + ], + "mag_angstrom": [ + 66.63744771612129, + 54.934363382294585, + 68.87249286560993 + ], + "mag_hz": [ + 38.20189863584585, + 26.49881430201914, + 40.436943785334485 + ], + "units_angstrom": "erg/(s cm2 Angstrom)", + "units_hz": "erg/(s cm2 Hz)" + }, + "z": { + "flux": [ + 1.9884770712558305e-32, + 1.3810180549144868e-35, + 5.204097063040092e-35 + ], + "mag_angstrom": [ + 57.05101291310126, + 64.94681599101527, + 63.506450911788306 + ], + "mag_hz": [ + 30.65376415331905, + 38.549567231233056, + 37.1092021520061 + ], + "units_angstrom": "erg/(s cm2 Angstrom)", + "units_hz": "erg/(s cm2 Hz)" + } +} \ No newline at end of file diff --git a/holodeck/ems/tests/test_basics.py b/holodeck/ems/tests/test_basics.py new file mode 100644 index 00000000..967599a2 --- /dev/null +++ b/holodeck/ems/tests/test_basics.py @@ -0,0 +1,139 @@ +""" +""" + +import numpy as np +import pytest + +# import holodeck as holo +from holodeck.ems import basics + + +def test_band_init_wlen(): + wlen = 123.0 # Angstrom + freq = 2.437337e16 # Hz + flux_wlen = "2.0 erg/(s cm^2 angstrom)" + + b1a = basics.Band('z', f"{wlen} angstrom", None, flux_wlen, None, "21.3 angstrom") + b1b = basics.Band('z', f"{wlen} angstrom", None, flux_wlen, None, None) + b2a = basics.Band('z', wlen, None, flux_wlen, None, 21.3) + b2b = basics.Band('z', wlen, None, flux_wlen, None, None) + + for band in [b1a, b1b, b2a, b2b]: + assert np.isclose(b1a.freq.to('Hz').value, freq, atol=0.0, rtol=1e-3) + + return + + +def test_band_init_freq(): + wlen = 123.0 # Angstrom + freq = 2.437337e16 # Hz + flux_freq = "2.0 erg/(s cm^2 Hz)" + b1a = basics.Band('z', None, f"{freq} Hz", None, flux_freq, None, "21.3e15 Hz") + b1b = basics.Band('z', None, f"{freq} Hz", None, flux_freq, None, None) + b2a = basics.Band('z', None, freq, None, flux_freq, 21.3e15) + b2b = basics.Band('z', None, freq, None, flux_freq, None) + + for band in [b1a, b1b, b2a, b2b]: + assert np.isclose(b1a.wlen.to('Angstrom').value, wlen, atol=0.0, rtol=1e-3) + + return + + +def test_band_init_fail(): + wlen = 123.0 # Angstrom + freq = 2.437337e16 # Hz + flux_freq = "2.0 erg/(s cm^2 Hz)" + bw_freq = "21.3e15 Hz" + + basics.Band('z', None, wlen, None, flux_freq, None, None) + basics.Band('z', None, wlen, None, flux_freq, None, bw_freq) + with pytest.raises(ValueError): + basics.Band('z', None, None, None, flux_freq, None, bw_freq) + + basics.Band('z', freq, None, None, flux_freq, None, None) + basics.Band('z', freq, None, None, flux_freq, None, bw_freq) + with pytest.raises(ValueError): + basics.Band('z', None, None, None, flux_freq, None, bw_freq) + + basics.Band('z', freq, None, None, flux_freq, None, bw_freq) + with pytest.raises(ValueError): + basics.Band('z', freq, None, None, None, None, bw_freq) + + return + + +def test_sdss_bands(): + bands = basics.SDSS_Bands() + + names = "ugriz" + for name in names: + assert name in bands.names + b1 = bands(name) + b2 = bands[name] + assert b1 == b2 + + # this is the reference flux for all bands, so by definition magnitude should be zero + flux = "3631.0 Jansky" + rmag = b1.flux_to_mag(flux, 'f') + assert np.isclose(rmag, 0.0) + check = b1.mag_to_flux(rmag, 'f') + assert np.isclose(flux, check) + + rmag = b1.flux_to_mag("3631.0", 'f', units='jansky') + assert np.isclose(rmag, 0.0) + check = b1.mag_to_flux(rmag, 'f') + assert np.isclose(flux, check) + + rmag = b1.flux_to_mag(3631.0, 'f', units='jansky') + assert np.isclose(rmag, 0.0) + check = b1.mag_to_flux(rmag, 'f') + assert np.isclose(flux, check) + + # ---- Absolute Magnitude and Luminosity + + lum = b1.abs_mag_to_lum(0, 'f').to('erg/s/Hz') + + rmag = b1.lum_to_abs_mag(lum, 'f', units='erg/(s Hz)') + assert np.isclose(rmag, 0.0, atol=1e-5) + + rmag = b1.lum_to_abs_mag(lum, 'f') + assert np.isclose(rmag, 0.0, atol=1e-5) + + rmag = b1.lum_to_abs_mag(4.34447e20, 'f', units='erg/(s Hz)') + assert np.isclose(rmag, 0.0, atol=1e-5) + + rmag = b1.lum_to_abs_mag("4.34447e20 erg/(s Hz)", 'f') + assert np.isclose(rmag, 0.0, atol=1e-5) + + return + + + +def test_sdss_bands__regression_2023_09_05(): + """Compare `flux_to_mag` calculation using current code, and previously computed data to check for regression. + """ + + bands = basics.SDSS_Bands() + + import json + from pathlib import Path + path = Path(__file__).resolve().parent + FNAME = "./basics_sdss_bands_2023-09-05.json" + data = json.load(open(path.joinpath(FNAME), 'r')) + + names = data.keys() + for band in names: + flux = data[band]['flux'] + mag_hz = data[band]['mag_hz'] + mag_angstrom = data[band]['mag_angstrom'] + units_hz = data[band]['units_hz'] + units_angstrom = data[band]['units_angstrom'] + + b1 = bands[band].flux_to_mag(flux, 'f', units=units_hz).value + b2 = bands[band].flux_to_mag(flux, 'w', units=units_angstrom).value + + assert np.allclose(b1, mag_hz, atol=0.0, rtol=1e-4) + assert np.allclose(b2, mag_angstrom, atol=0.0, rtol=1e-4) + + return + diff --git a/holodeck/utils.py b/holodeck/utils.py index 9853f78f..37b541f4 100644 --- a/holodeck/utils.py +++ b/holodeck/utils.py @@ -18,7 +18,7 @@ import subprocess import warnings from pathlib import Path -from typing import Optional, Tuple, Union, List, Callable, TypeVar, Any # , TypeAlias # , Sequence, +from typing import Optional, Tuple, Union, List, # Callable, TypeVar, Any # , TypeAlias # , Sequence, # try: # from typing import ParamSpec @@ -1851,17 +1851,32 @@ def angs_from_sepa(sepa, dcom, redz): return angs -def eddington_luminosity(mass, eps=0.1): - ledd = EDDT * mass / eps - return ledd +def eddington_accretion(mass, eps=0.1): + """Eddington Accretion rate, $\\dot{M}_{Edd} = L_{Edd}/\\epsilon c^2$. + Arguments + --------- + mass : array_like of scalar + BH Mass. + eps : array_like of scalar + Efficiency parameter. -def eddington_accretion(mass, eps=0.1): + Returns + ------- + mdot : array_like of scalar + Eddington accretion rate. + + """ edd_lum = eddington_luminosity(mass, eps=eps) - mdot = edd_lum / np.square(SPLC) + # NOTE: no `epsilon` (efficiency) in this equation, because included in `eddington_luminosity` + mdot = edd_lum/np.square(SPLC) return mdot +def eddington_luminosity(mass, eps=0.1): + ledd = EDDT * mass / eps + return ledd + # ================================================================================================= # ==== Gravitational Waves ==== diff --git a/notebooks/devs/ems/basics.ipynb b/notebooks/devs/ems/basics.ipynb new file mode 100644 index 00000000..f423af0a --- /dev/null +++ b/notebooks/devs/ems/basics.ipynb @@ -0,0 +1,112 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import astropy as ap\n", + "import astropy.units\n", + "import matplotlib.pyplot as plt\n", + "\n", + "import holodeck as holo\n", + "import holodeck.ems" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bands = holo.ems.basics.SDSS_Bands()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "lum = (ap.units.Quantity(\"3631.0 Jansky\") * (4.0 * np.pi * (10.0 * ap.units.parsec)**2)).to('erg / s / Hz')\n", + "lum" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bands('u').flux_to_mag(3631.0, type='f', units='jansky')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "bands('u').lum_to_abs_mag(lum, type='f')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Generate Regression-Test data for basics.SDSS_Bands" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# from pathlib import Path\n", + "# import zcode.inout as zio\n", + "\n", + "# results = {}\n", + "\n", + "# units_hz = 'erg/(s cm2 Hz)'\n", + "# units_angstrom = 'erg/(s cm2 Angstrom)'\n", + "\n", + "# for band in \"ugriz\":\n", + "# flux = 10.0 ** np.random.uniform(-36, -30, 3)\n", + "# b1 = bands[band].flux_to_mag(flux, 'f', units=units_hz).value\n", + "# b2 = bands[band].flux_to_mag(flux, 'w', units=units_angstrom).value\n", + "# results[band] = dict(flux=flux, mag_hz=b1, mag_angstrom=b2, units_hz=units_hz, units_angstrom=units_angstrom)\n", + " \n", + "# fname = 'basics_sdss_bands_2023-09-05.json'\n", + "# path = Path(holo.ems.basics.__file__).parent.resolve()\n", + "# print(path)\n", + "# print(results)\n", + "# # rv = zio.str_format_dict(results, file=open(path.joinpath(\"tests\", fname), 'w'))\n", + "# print(zio.str_format_dict(results))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "py311", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/notebooks/devs/ems/drw.ipynb b/notebooks/devs/ems/drw.ipynb new file mode 100644 index 00000000..e5400bb3 --- /dev/null +++ b/notebooks/devs/ems/drw.ipynb @@ -0,0 +1,115 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "import holodeck as holo\n", + "import holodeck.ems\n", + "import holodeck.plot\n", + "from holodeck.constants import MSOL" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "num = 3000\n", + "mass = (10.0**np.random.uniform(6, 10, num)) * MSOL\n", + "fedd = np.random.uniform(0.0, 1.0, num)\n", + "imag, taus, sfis = holo.ems.drw.drw_params(mass, fedd, eps=1.0, samples=False)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "num = 3\n", + "mass = (10.0**np.random.uniform(6, 10, num)) * MSOL\n", + "# fedd = np.random.uniform(0.0, 1.0, num)\n", + "fedd = np.ones_like(mass)\n", + "imag, taus, sfis = holo.ems.drw.drw_params(mass, fedd, eps=1.0, samples=False)\n", + "print(mass/MSOL)\n", + "print(imag)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, axes = holo.plot.figax(figsize=[12, 4], ncols=3, xlabel='Mass [Msol]', wspace=0.35, yscale='lin')\n", + "\n", + "values = [imag, np.log10(taus), np.log10(sfis)]\n", + "smap = holo.plot.smap(fedd, log=False)\n", + "colors = smap.to_rgba(fedd)\n", + "\n", + "labels = ['$M_i$', r'$\\tau$ [days]', 'SF$_\\infty$ [mag]']\n", + "for ii, (val, ax) in enumerate(zip(values, axes)):\n", + " ax.scatter(mass/MSOL, val, color=colors)\n", + " ax.set(ylabel=labels[ii])\n", + " plt.colorbar(smap, ax=ax, orientation='horizontal', label='Fedd')\n", + "\n", + "axes[0].set(yscale='linear')\n", + "axes[0].invert_yaxis()\n", + "axes[0].set(ylim=[-25, -33])\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig, axes = holo.plot.figax(figsize=[12, 4], ncols=2, wspace=0.35, \n", + " xlabel='Mass [Msol]', yscale='lin', ylim=[-23, -30], xlim=[1e8, 1e10])\n", + "\n", + "values = [np.log10(taus), np.log10(sfis)]\n", + "\n", + "labels = [r'$\\tau$ [days]', 'SF$_\\infty$ [mag]']\n", + "for ii, (val, ax) in enumerate(zip(values, axes)):\n", + " print(holo.utils.stats(val))\n", + " smap = holo.plot.smap(val, log=False, cmap='Spectral_r')\n", + " colors = smap.to_rgba(val)\n", + " ax.scatter(mass/MSOL, imag, color=colors)\n", + " ax.set(ylabel=\"$M_i$\")\n", + " plt.colorbar(smap, ax=ax, orientation='horizontal', label=labels[ii])\n", + "\n", + "axes[0].invert_yaxis()\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "py311", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/notebooks/devs/ems/lsst_gwb/lsst_gwb_simple.ipynb b/notebooks/devs/ems/lsst_gwb/lsst_gwb_simple.ipynb new file mode 100644 index 00000000..24f7a2b0 --- /dev/null +++ b/notebooks/devs/ems/lsst_gwb/lsst_gwb_simple.ipynb @@ -0,0 +1,544 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "\n", + "import numpy as np\n", + "import scipy as sp\n", + "import matplotlib.pyplot as plt\n", + "import h5py\n", + "\n", + "import kalepy as kale\n", + "\n", + "import holodeck as holo\n", + "from holodeck.constants import MSOL, PC, YR, GYR, SPLC, EDDT" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 15yr Population Posteriors" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Load chains from 15yr Binary Astrophysics analysis to get population parameter posteriors" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "path_data = Path(\"./data/astroprior_hdall\").resolve()\n", + "print(path_data)\n", + "assert path_data.is_dir()\n", + "fname_pars = path_data.joinpath(\"pars.txt\")\n", + "fname_chains = path_data.joinpath(\"chain_1.0.txt\")\n", + "print(fname_pars)\n", + "print(fname_chains)\n", + "assert fname_chains.is_file() and fname_pars.is_file()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "chain_pars = np.loadtxt(fname_pars, dtype=str)\n", + "chains = np.loadtxt(fname_chains)\n", + "npars = len(chain_pars)\n", + "\n", + "# Get maximum likelihood parameters (estimate using KDE)\n", + "mlpars = {}\n", + "fig, axes = plt.subplots(figsize=[10, 1.5*npars], nrows=npars)\n", + "plt.subplots_adjust(hspace=0.75)\n", + "for ii, ax in enumerate(axes):\n", + " ax.set(xlabel=chain_pars[ii])\n", + " vals = chains[:, ii]\n", + " extr = holo.utils.minmax(vals)\n", + " xx, yy = kale.density(vals, reflect=extr)\n", + " kale.dist1d(chains[:, ii], ax=ax, density=True, carpet=1000)\n", + " idx = np.argmax(yy)\n", + " xmax = xx[idx]\n", + " ax.axvline(xmax, color='firebrick')\n", + " mlpars[chain_pars[ii]] = xmax\n", + " \n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(f\"Maximum Likelihood binary population parameters:\")\n", + "for kk, vv in mlpars.items():\n", + " print(f\"\\t{kk:>20s}: {vv:+.2e}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Generate Population with ML Parameters" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Construct model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose the appropriate Parameter Space (from 15yr astro analysis)\n", + "pspace = holo.param_spaces.PS_Uniform_09B\n", + "# Load SAM and hardening model for desired parameters\n", + "sam, hard = pspace.model_for_params(mlpars)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# compare a couple of the parameters to make sure things look right\n", + "print(hard._target_time/GYR, mlpars['hard_time'])\n", + "assert np.isclose(hard._target_time/GYR, mlpars['hard_time'])\n", + "print(sam._gsmf._phi0, mlpars['gsmf_phi0'])\n", + "assert np.isclose(sam._gsmf._phi0, mlpars['gsmf_phi0'])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Calculate number of binaries in a target frequency (period) range.\n", + "Takes about 1 minute" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose range of orbital periods of interest\n", + "tvals = [10.0, 0.1] # yr\n", + "NFBINS = 10\n", + "print(f\"Considering orbital periods between {tvals} yrs, {NFBINS} bins\")\n", + "# convert to frequencies\n", + "fobs_orb_edges = 1 / np.array(tvals) # 1/yr\n", + "# construct bins \n", + "fobs_orb_edges = np.logspace(*np.log10(fobs_orb_edges/YR), NFBINS+1) # 1/sec\n", + "fobs_orb_cents = holo.utils.midpoints(fobs_orb_edges)\n", + "fobs = fobs_orb_cents\n", + "\n", + "# calculate (differential) number of binaries\n", + "redz_final, diff_num = holo.sams.cyutils.dynamic_binary_number_at_fobs(\n", + " fobs_orb_cents, sam, hard, holo.cosmo\n", + ")\n", + "# integrate to find total number of binaries in each bin\n", + "edges = [sam.mtot, sam.mrat, sam.redz, fobs_orb_edges]\n", + "number = holo.sams.cyutils.integrate_differential_number_3dx1d(edges, diff_num)\n", + "print(f\"Loaded {number.sum():.1e} binaries across frequency range\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "temp = number.sum(axis=(0, 1, 2))\n", + "fig, ax = plt.subplots()\n", + "\n", + "xx = fobs_orb_cents*YR\n", + "yy = temp/np.diff(fobs_orb_edges*YR)\n", + "ax.plot(xx, yy)\n", + "\n", + "ax.set(xscale='log', yscale='log', ylabel='Number of Binaries ($dN/df$)', xlabel='Orbial Frequency [1/yr]')\n", + "tw = holo.plot._twin_hz(ax)\n", + "tw.set_xlabel('orbital frequency [nHz]')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Observability\n", + "\n", + "Choose which bins of SAM population are 'observable'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fedd = 0.1 # eddington fraction\n", + "bcorr = 0.1 # bolometric correction, bolometric ==> optical\n", + "\n", + "# LSST V-band sensitivity [erg/s/cm^2/Hz]\n", + "# see: https://ui.adsabs.harvard.edu/abs/2019MNRAS.485.1579K/abstract\n", + "flux_sens_lsst = 3.0e-30\n", + "vband_wlen = 551.0e-7 # [cm]\n", + "\n", + "# get V-band frequency\n", + "vband_freq = SPLC/(vband_wlen) # [Hz]\n", + "\n", + "# get bin-center values for population\n", + "mtot = holo.utils.midpoints(sam.mtot)\n", + "mrat = holo.utils.midpoints(sam.mrat)\n", + "redz = holo.utils.midpoints(sam.redz)\n", + "# convert redshift to luminosity-distance\n", + "dlum = holo.cosmo.z_to_dlum(redz)\n", + "\n", + "# calculate luminosity of binaries based on Eddington fraction and bolometric correction\n", + "lum = EDDT * mtot * fedd * bcorr # [erg/s]\n", + "\n", + "# calculate flux at observer\n", + "# TODO: should really divide by the width of the V-band\n", + "flux_tot = lum[:, np.newaxis] / (4*np.pi*dlum[np.newaxis, :]**2) / vband_freq\n", + "# get the flux of the secondary, assume that it is what's needed\n", + "flux_sec = flux_tot[:, np.newaxis, :] * (mrat / (1.0 + mrat))[np.newaxis, :, np.newaxis]\n", + "\n", + "# select \"observable\" systems\n", + "obs_flag = (flux_sec > flux_sens_lsst)\n", + "num_obs = np.sum(obs_flag[..., np.newaxis]*number)\n", + "num_all = np.sum(number)\n", + "frac_obs = num_obs / num_all\n", + "print(f\"observable: {num_obs:.2e}/{num_all:.2e} = {frac_obs:.2e}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Detectability (Test) Data from Caitlin " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fname_data = Path(\"./data/export_for_gwb_test.txt\")\n", + "fname_data = fname_data.absolute().resolve()\n", + "print(fname_data, fname_data.exists())\n", + "with open(fname_data, 'r') as input:\n", + " det_header = None\n", + " for ii, line in enumerate(input.readlines()):\n", + " line = line.strip()\n", + " print(line)\n", + " if det_header is None:\n", + " if not line.startswith('#'):\n", + " raise ValueError(\n", + " \"First line of file should have stared with a comment including header \"\n", + " \"information about the columns!\"\n", + " )\n", + " det_header = line.strip(' #').split(\" \")\n", + " det_header = [head.strip() for head in det_header]\n", + " print(len(det_header), det_header)\n", + " \n", + " if ii > 3:\n", + " break\n", + "\n", + "det_data = np.loadtxt(fname_data)\n", + "print(f\"{det_data.shape=}\")\n", + "\n", + "injected = (det_data[:, 11] > 0)\n", + "print(\"injected: \", holo.utils.frac_str(injected))\n", + "detected = (det_data[:, 12] > 0)\n", + "print(\"detected: \", holo.utils.frac_str(detected))\n", + "print(\" both : \", holo.utils.frac_str(detected & injected))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Look at the data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "indices = [0, 2, 4, 6, 8]\n", + "\n", + "num = len(indices)\n", + "fig, axes = plt.subplots(figsize=[20, num*5], ncols=4, nrows=num, sharex='row')\n", + "\n", + "for ii, axrow in enumerate(axes):\n", + " idx = indices[ii]\n", + "\n", + " xx = det_data[:, idx]\n", + " yy = det_data[:, idx+1]\n", + "\n", + " for ax in axrow: \n", + " ax.set(xlabel=det_header[idx], ylabel=det_header[idx+1])\n", + "\n", + " ax = axrow[0]\n", + " if ii == 0:\n", + " ax.set(title='all')\n", + " ax.scatter(xx, yy, alpha=0.2, s=5)\n", + " \n", + " ax = axrow[1]\n", + " if ii == 0:\n", + " ax.set(title='injected')\n", + " ax.scatter(xx, yy, alpha=0.2, s=14)\n", + " ax.scatter(xx[injected], yy[injected], alpha=0.75, marker='x', s=8, lw=0.5)\n", + "\n", + " ax = axrow[2]\n", + " if ii == 0:\n", + " ax.set(title='detected')\n", + " ax.scatter(xx, yy, alpha=0.2, s=14)\n", + " ax.scatter(xx[detected], yy[detected], alpha=0.75, marker='x', s=8, lw=0.5)\n", + " \n", + " ax = axrow[3]\n", + " if ii == 0:\n", + " ax.set(title='both')\n", + " ax.scatter(xx, yy, alpha=0.2, s=14)\n", + " ax.scatter(xx[detected & injected], yy[detected & injected], alpha=0.75, marker='x', s=8, lw=0.5)\n", + " \n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "indices = [0, 2, 4, 6, 8]\n", + "\n", + "num = len(indices)\n", + "fig, axes = plt.subplots(figsize=[20, num*5], ncols=2, nrows=num, sharex='row')\n", + "kwargs = dict(alpha=0.5, lw=2.0, density=True, histtype='step')\n", + "\n", + "for ii, axrow in enumerate(axes):\n", + " idx = indices[ii]\n", + " bins = 20\n", + "\n", + " for jj, ax in enumerate(axrow):\n", + " ax.set(yscale='log', xlabel=det_header[idx+jj], ylabel='Number')\n", + " vals = det_data[:, idx+jj]\n", + " hist, bins, patches = ax.hist(vals, bins=bins, label='all', **kwargs)\n", + " hist, bins, patches = ax.hist(vals[injected], bins=bins, label='injected', **kwargs)\n", + " hist, bins, patches = ax.hist(vals[detected], bins=bins, label='detected', **kwargs)\n", + " hist, bins, patches = ax.hist(vals[injected & detected], bins=bins, label='both', **kwargs)\n", + "\n", + " ax.legend()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Calculate a 'detectability' metric" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "NAMPS = 7\n", + "NPERS = 9\n", + "\n", + "def get_idx(key, header):\n", + " for ii, hh in enumerate(header):\n", + " if key.lower() in hh.lower():\n", + " return ii\n", + " else:\n", + " return None\n", + "\n", + "# grab amplitude and period data\n", + "amp_idx = get_idx('amp_in', det_header)\n", + "period_idx = get_idx('period_in', det_header)\n", + "det_amps = det_data[:, amp_idx]\n", + "# convert periods from [day] to [yr]\n", + "det_pers = det_data[:, period_idx]*24*60*60/YR\n", + "\n", + "# Choose a 2D grid of bin-edges based on the detected amplitudes and periods\n", + "sel_flag = injected & detected\n", + "amp_edges = np.linspace(*holo.utils.minmax(det_amps[sel_flag]), NAMPS)\n", + "per_edges = np.linspace(*holo.utils.minmax(det_pers[sel_flag]), NPERS)\n", + "print(f\"{amp_edges=}\")\n", + "print(f\"{per_edges=}\")\n", + "bins = [amp_edges, per_edges]\n", + "\n", + "# find the number of points in each bin\n", + "num_all, *_ = sp.stats.binned_statistic_2d(\n", + " det_amps, det_pers, np.ones_like(det_amps), statistic='sum', bins=bins\n", + ")\n", + "# find the number of injected & detected points in each bin\n", + "num_det, *_ = sp.stats.binned_statistic_2d(\n", + " det_amps, det_pers, sel_flag*np.ones_like(det_amps), statistic='sum', bins=bins\n", + ")\n", + "# The detection fraction is the number of injected & detected points divided by all points\n", + "# TODO: should denominator just be the number of injected points??? How to handle false-positives???\n", + "det_frac = num_det / num_all\n", + "\n", + "fig, ax = plt.subplots(figsize=[10, 7])\n", + "ax.set(xlabel='amplitude [frac]', ylabel='period [yr]')\n", + "pcm = ax.pcolormesh(*bins, det_frac.T, shading='auto')\n", + "plt.colorbar(pcm, ax=ax, label='detection fraction')\n", + "\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Calculate LSST Detections" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have a number of binaries `number` over a grid of total-mass, mass-ratio, redshift, and orbital frequency. The shape is (M, Q, Z, F). The bin-center values for each dimension are given in:\n", + " * `mtot` (M,) [gram]\n", + " * `mrat` (Q,) [-]\n", + " * `redz` (Z,) [-]\n", + " * `fobs` (F,) [1/sec]\n", + "\n", + "There is a boolean grid of which bins are observable given in `obs_flag`, with the same shape as `number`.\n", + "\n", + "We want to determine what fraction of binaries in each bin are detectable in LSST variability surveys." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have a grid of detectability fractions for periodic variable AGN in `det_frac` with shape (A, P). This is over a grid of variability amplitudes given by the array `amp_edges` shaped (A+1,), and orbital periods `per_edges` shaped (P+1,). We need to map the simulated binaries to this parameter space." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# assume that the variability amplitude exactly equals the mass-ratio\n", + "bin_amp = mrat[np.newaxis, :, np.newaxis, np.newaxis] * np.ones_like(number)\n", + "\n", + "# assume that the variability period is exactly the orbital period\n", + "_per = (1/fobs/YR)\n", + "bin_per = _per[np.newaxis, np.newaxis, np.newaxis, :] * np.ones_like(number)\n", + "\n", + "# convert to 1D arrays, and select out the 'observable' binaries\n", + "bin_amp = bin_amp[obs_flag].flatten()\n", + "bin_per = bin_per[obs_flag].flatten()\n", + "bin_num = number[obs_flag].flatten()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have the binaries in terms of the variable-detectability parameter space. So we just need to find the detectability fraction (`det_frac`) for each binary grid-point now." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# returned indices `idx` give the bin number\n", + "amp_idx = np.digitize(bin_amp, amp_edges) - 1\n", + "per_idx = np.digitize(bin_per, per_edges) - 1\n", + "# idx values of -1 mean the value is below the lowest bin, values of B+1 (for B bins) are above the highest bin\n", + "# put amplitudes above the highest bin into the highest bin\n", + "nbins = amp_edges.size - 1\n", + "amp_idx[amp_idx >= nbins] = nbins - 1\n", + "# set amplitudes below smallest bin to be invalid, i.e. select only values above the lowest bin\n", + "sel_amp = (amp_idx >= 0)\n", + "# put periods below the lowest bin, into the lowest bin\n", + "per_idx[per_idx < 0] = 0\n", + "# set periods above highest bin to be invalid, i.e. select only values below the highest bin\n", + "nbins = per_edges.size - 1\n", + "print(per_edges.size, holo.utils.stats(per_idx))\n", + "sel_per = (per_idx < nbins)\n", + "\n", + "# select valid entires\n", + "sel = sel_amp & sel_per\n", + "amp_idx = amp_idx[sel]\n", + "per_idx = per_idx[sel]\n", + "# grab the corresponding numbers of binaries in each of these 'selected' bins\n", + "sel_bin_num = bin_num[sel]\n", + "print(f\"{holo.utils.stats(sel_bin_num)=}\")\n", + "\n", + "# convert from indices in each dimension, to an index for the flattened array\n", + "sel_dfracs = det_frac[(amp_idx, per_idx)]\n", + "print(f\"{holo.utils.stats(sel_dfracs)=}\")\n", + "\n", + "# find the total number of detectable binaries\n", + "# multiply the number of binaries in each bin, by the detection fraction in that bin\n", + "num_det_bins = sel_bin_num * sel_dfracs\n", + "print(f\"{num_det_bins.sum()=:.2e}\")\n", + "num_all_bins = number.sum()\n", + "frac_det_bins = num_det_bins.sum() / num_all_bins\n", + "print(f\"Total detection fraction: {frac_det_bins:.2e}\")\n", + "\n", + "# remind us the fraction of binaries that were 'observable'\n", + "print(f\"Total 'observability' fraction: {frac_obs:.2e}\")\n", + "\n", + "frac_obs_det = num_det_bins.sum() / num_obs\n", + "print(f\"Det frac of observable: {frac_obs_det:.2e}\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "py311", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/notebooks/devs/ems/lsst_gwb/pop_drw_params.ipynb b/notebooks/devs/ems/lsst_gwb/pop_drw_params.ipynb new file mode 100644 index 00000000..135b4e00 --- /dev/null +++ b/notebooks/devs/ems/lsst_gwb/pop_drw_params.ipynb @@ -0,0 +1,222 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "\n", + "import numpy as np\n", + "import scipy as sp\n", + "import matplotlib.pyplot as plt\n", + "import h5py\n", + "\n", + "import kalepy as kale\n", + "\n", + "import holodeck as holo\n", + "import holodeck.ems\n", + "from holodeck.constants import MSOL, PC, YR, GYR, SPLC, EDDT" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Get population parameters from 15yr data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "path_data = \"/Users/lzkelley/programs/nanograv/15yr_astro_data/phenom/ceffyl_chains/astroprior_hdall\"\n", + "path_data = Path(path_data).resolve()\n", + "print(path_data)\n", + "assert path_data.is_dir()\n", + "fname_pars = path_data.joinpath(\"pars.txt\")\n", + "fname_chains = path_data.joinpath(\"chain_1.0.txt\")\n", + "print(fname_pars)\n", + "print(fname_chains)\n", + "assert fname_chains.is_file() and fname_pars.is_file()\n", + "\n", + "chain_pars = np.loadtxt(fname_pars, dtype=str)\n", + "chains = np.loadtxt(fname_chains)\n", + "npars = len(chain_pars)\n", + "\n", + "# Get maximum likelihood parameters (estimate using KDE)\n", + "mlpars = {}\n", + "fig, axes = plt.subplots(figsize=[10, 1.5*npars], nrows=npars)\n", + "plt.subplots_adjust(hspace=0.75)\n", + "for ii, ax in enumerate(axes):\n", + " ax.set(xlabel=chain_pars[ii])\n", + " vals = chains[:, ii]\n", + " extr = holo.utils.minmax(vals)\n", + " xx, yy = kale.density(vals, reflect=extr)\n", + " kale.dist1d(chains[:, ii], ax=ax, density=True, carpet=1000)\n", + " idx = np.argmax(yy)\n", + " xmax = xx[idx]\n", + " ax.axvline(xmax, color='firebrick')\n", + " mlpars[chain_pars[ii]] = xmax\n", + " \n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(f\"Maximum Likelihood binary population parameters:\")\n", + "for kk, vv in mlpars.items():\n", + " print(f\"\\t{kk:>20s}: {vv:+.4e}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Construct population" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose the appropriate Parameter Space (from 15yr astro analysis)\n", + "pspace = holo.param_spaces.PS_Uniform_09B\n", + "# Load SAM and hardening model for desired parameters\n", + "sam, hard = pspace.model_for_params(mlpars)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Choose range of orbital periods of interest\n", + "tvals = [10.0, 0.1] # yr\n", + "NFBINS = 10\n", + "print(f\"Considering orbital periods between {tvals} yrs, {NFBINS} bins\")\n", + "# convert to frequencies\n", + "fobs_orb_edges = 1 / np.array(tvals) # 1/yr\n", + "# construct bins \n", + "fobs_orb_edges = np.logspace(*np.log10(fobs_orb_edges/YR), NFBINS+1) # 1/sec\n", + "fobs_orb_cents = holo.utils.midpoints(fobs_orb_edges)\n", + "fobs = fobs_orb_cents\n", + "\n", + "# calculate (differential) number of binaries\n", + "redz_final, diff_num = holo.sams.cyutils.dynamic_binary_number_at_fobs(\n", + " fobs_orb_cents, sam, hard, holo.cosmo\n", + ")\n", + "# integrate to find total number of binaries in each bin\n", + "edges = [sam.mtot, sam.mrat, sam.redz, fobs_orb_edges]\n", + "number = holo.sams.cyutils.integrate_differential_number_3dx1d(edges, diff_num)\n", + "print(f\"Loaded {number.sum():.1e} binaries across frequency range\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Calculate DRW parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# breaker()\n", + "fedd_num = 10\n", + "# we dont care about orbital frequency for this, so ignore\n", + "cents = [holo.utils.midpoints(ee, log=True) for ee in edges[:-1]]\n", + "mesh = [mm.flatten() for mm in np.meshgrid(*cents, indexing='ij')]\n", + "size = mesh[0].size\n", + "shape = (size, fedd_num, 2)\n", + "fedd = holo.utils.log_normal_base_10(0.1, 0.5, shape)\n", + "fedd[fedd > 1.0] = 1.0/fedd[fedd > 1.0]\n", + "fedd = fedd.reshape(-1, 2)\n", + "mesh = [mm[:, np.newaxis] * np.ones(shape[:-1]) for mm in mesh]\n", + "mt, mr, rz = [mm.flatten() for mm in mesh]\n", + "m1, m2 = holo.utils.m1m2_from_mtmr(mt, mr)\n", + "\n", + "num = number.sum(axis=-1).flatten()\n", + "num = num[:, np.newaxis] * np.ones(shape[:-1])\n", + "num = num.flatten() / fedd_num\n", + "\n", + "scatter = False\n", + "imag_1, taus_1, sfis_1 = holo.ems.drw.drw_params(m1, fedd[:, 0], scatter=scatter)\n", + "imag_2, taus_2, sfis_2 = holo.ems.drw.drw_params(m2, fedd[:, 1], scatter=scatter)\n", + "\n", + "scatter = True\n", + "imag_1_scatter, taus_1_scatter, sfis_1_scatter = holo.ems.drw.drw_params(m1, fedd[:, 0], scatter=scatter)\n", + "imag_2_scatter, taus_2_scatter, sfis_2_scatter = holo.ems.drw.drw_params(m2, fedd[:, 1], scatter=scatter)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fname = \"./drw_params.hdf5\"\n", + "fname = Path(fname).resolve()\n", + "with h5py.File(fname, 'w') as out:\n", + " out.create_dataset(\"m1\", data=m1)\n", + " out.create_dataset(\"m2\", data=m2)\n", + " out.create_dataset(\"num\", data=num)\n", + " out.create_dataset(\"redz\", data=rz)\n", + " out.create_dataset(\"fedd1\", data=fedd[:, 0])\n", + " out.create_dataset(\"fedd2\", data=fedd[:, 1])\n", + "\n", + " group = out.create_group(\"mean\")\n", + " group.create_dataset(\"imag1\", data=imag_1)\n", + " group.create_dataset(\"imag2\", data=imag_2)\n", + " group.create_dataset(\"taus1\", data=taus_1)\n", + " group.create_dataset(\"taus2\", data=taus_2)\n", + " group.create_dataset(\"sfis1\", data=sfis_1)\n", + " group.create_dataset(\"sfis2\", data=sfis_2)\n", + "\n", + " group = out.create_group(\"scatter\")\n", + " group.create_dataset(\"imag1\", data=imag_1_scatter)\n", + " group.create_dataset(\"imag2\", data=imag_2_scatter)\n", + " group.create_dataset(\"taus1\", data=taus_1_scatter)\n", + " group.create_dataset(\"taus2\", data=taus_2_scatter)\n", + " group.create_dataset(\"sfis1\", data=sfis_1_scatter)\n", + " group.create_dataset(\"sfis2\", data=sfis_2_scatter)\n", + " \n", + "print(f\"Saved to {fname}, size {holo.utils.get_file_size(fname)}\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "py311", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.4" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}