From d96e8c2afdc2a1eccbc94fb0a4d59b5cb2ac6eec Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Wed, 22 May 2024 01:49:46 +0200 Subject: [PATCH] bug fixes --- maria/atmosphere/sim.py | 19 +++-- maria/instrument/__init__.py | 8 +- maria/instrument/bands/__init__.py | 95 ++++++++++++--------- maria/instrument/bands/configs/abs.yml | 2 +- maria/instrument/bands/configs/act.yml | 24 +++--- maria/instrument/bands/configs/alma.yml | 40 ++++----- maria/instrument/bands/configs/atlast.yml | 24 +++--- maria/instrument/bands/configs/mustang2.yml | 4 +- maria/instrument/bands/format.csv | 7 ++ maria/instrument/configs.yml | 1 + maria/map/__init__.py | 95 ++++++++++++++------- maria/noise/__init__.py | 10 ++- maria/sim/__init__.py | 12 +-- 13 files changed, 207 insertions(+), 134 deletions(-) create mode 100644 maria/instrument/bands/format.csv diff --git a/maria/atmosphere/sim.py b/maria/atmosphere/sim.py index ca7b906..c0eaed2 100644 --- a/maria/atmosphere/sim.py +++ b/maria/atmosphere/sim.py @@ -35,11 +35,13 @@ def _initialize_2d_atmosphere( ) self.atmosphere.layers = [] - depths = enumerate(self.turbulent_layer_depths) - if self.verbose: - depths = tqdm(depths, desc="Initializing atmospheric layers") + depths = tqdm( + self.turbulent_layer_depths, + desc="Initializing atmospheric layers", + disable=not self.verbose, + ) - for _, layer_depth in depths: + for layer_depth in depths: layer_res = ( self.instrument.physical_fwhm(z=layer_depth).min() / min_atmosphere_beam_res @@ -72,7 +74,12 @@ def _simulate_2d_turbulence(self): (len(self.atmosphere.layers), self.instrument.n_dets, self.plan.n_time) ) - pbar = tqdm(enumerate(self.atmosphere.layers), disable=not self.verbose) + pbar = tqdm( + enumerate(self.atmosphere.layers), + desc="Generating atmosphere", + disable=not self.verbose, + ) + for layer_index, layer in pbar: layer.generate() layer_data[layer_index] = sp.interpolate.interp1d( @@ -83,7 +90,7 @@ def _simulate_2d_turbulence(self): bounds_error=False, fill_value="extrapolate", )(self.boresight.time) - pbar.set_description(f"Generating atmosphere (z={layer.depth:.00f}m)") + # pbar.set_description(f"Generating atmosphere (z={layer.depth:.00f}m)") return layer_data diff --git a/maria/instrument/__init__.py b/maria/instrument/__init__.py index 8a93355..e42239f 100644 --- a/maria/instrument/__init__.py +++ b/maria/instrument/__init__.py @@ -29,6 +29,10 @@ INSTRUMENT_CONFIGS = read_yaml(f"{here}/configs.yml") +for name, config in INSTRUMENT_CONFIGS.items(): + config["aliases"] = config.get("aliases", []) + config["aliases"].append(name.lower()) + INSTRUMENT_DISPLAY_COLUMNS = [ "description", # "field_of_view", @@ -50,7 +54,7 @@ def get_instrument(instrument_name="default", **kwargs): """ if instrument_name: for key, config in INSTRUMENT_CONFIGS.items(): - if instrument_name in config.get("aliases", []): + if instrument_name.lower() in config.get("aliases", []): instrument_name = key if instrument_name not in INSTRUMENT_CONFIGS.keys(): raise InvalidInstrumentError(instrument_name) @@ -149,7 +153,7 @@ def get_subarrays(instrument_config): subarray["bands"][band_name] = all_bands[band_name] for param in subarray_params_to_inherit: - if param in config: + if (param in config) and (param not in subarray): subarray[param] = config[param] if "bands" not in subarray: diff --git a/maria/instrument/bands/__init__.py b/maria/instrument/bands/__init__.py index 606cae6..f2af403 100644 --- a/maria/instrument/bands/__init__.py +++ b/maria/instrument/bands/__init__.py @@ -8,21 +8,14 @@ import pandas as pd from ...atmosphere.spectrum import Spectrum -from ...constants import T_CMB, c -from ...functions import planck_spectrum, rayleigh_jeans_spectrum +from ...constants import T_CMB, c, k_B +from ...functions import planck_spectrum from ...io import flatten_config, read_yaml -BAND_FIELDS = { - "center": {"units": "GHz", "dtype": "float"}, - "width": {"units": "GHz", "dtype": "float"}, - "shape": {"units": None, "dtype": "str"}, - "efficiency": {"units": None, "dtype": "float"}, - "sensitivity": {"units": "K sqrt(s)", "dtype": "float"}, - "NEP": {"units": "pW sqrt(s)", "dtype": "float"}, -} - here, this_filename = os.path.split(__file__) +FIELD_FORMATS = pd.read_csv(f"{here}/format.csv", index_col=0) + all_bands = {} for path in glob.glob(f"{here}/configs/*.yml"): tag = os.path.split(path)[1].split(".")[0] @@ -115,16 +108,35 @@ def __init__( sensitivity, kind=sensitivity_kind ) # this sets the NEP automatically + @property + def summary(self): + summary = pd.Series(index=FIELD_FORMATS.index, dtype=str) + + for field, entry in FIELD_FORMATS.iterrows(): + value = getattr(self, field) + + if entry["dtype"] == "float": + if entry["format"] == "e": + s = f"{value:.02e}" + else: + s = f"{value}" + + if entry.units != "none": + s = f"{s} {entry.units}" + + elif entry["dtype"] == "str": + s = f"{value}" + + summary[field] = s + + return summary + def __repr__(self): + summary = self.summary parts = [] - for field, d in BAND_FIELDS.items(): - value = getattr(self, field) - if d["dtype"] == "str": - s = f"{field}='{value}'" - elif d["units"] is not None: - s = f"{field}={value} {d['units']}" - else: - s = f"{field}={value}" + for field, entry in FIELD_FORMATS.iterrows(): + value = summary[field] + s = f"{field}='{value}'" if entry["dtype"] == str else f"{field}={value}" parts.append(s) return f"Band({', '.join(parts)})" @@ -219,36 +231,40 @@ def passband(self, nu): @property def dP_dTRJ(self) -> float: """ - In watts per kelvin Rayleigh-Jeans, assuming perfect transmission. + In picowatts per kelvin Rayleigh-Jeans, assuming perfect transmission. """ nu = np.linspace(self.nu_min, self.nu_max, 256) - dI_dTRJ = rayleigh_jeans_spectrum( - nu=1e9 * nu, T=1 - ) # this is the same as the derivative - dP_dTRJ = np.trapz(dI_dTRJ * self.passband(nu), 1e9 * nu) + # dI_dTRJ = rayleigh_jeans_spectrum(nu=1e9 * nu, T=1) # this is the same as the derivative + # dP_dTRJ = np.trapz(dI_dTRJ * self.passband(nu), 1e9 * nu) + dP_dTRJ = k_B * np.trapz(self.passband(nu), 1e9 * nu) - return self.efficiency * dP_dTRJ + return 1e12 * self.efficiency * dP_dTRJ @property def dP_dTCMB(self) -> float: """ - In watts per kelvin CMB, assuming perfect transmission. + In picowatts per kelvin CMB, assuming perfect transmission. """ - eps = 1e-3 + eps = 1e-2 + delta_T = np.array([-eps / 2, eps / 2]) nu = np.linspace(self.nu_min, self.nu_max, 256) + TRJ = ( + planck_spectrum(nu=1e9 * nu, T=T_CMB + delta_T[:, None]) + * c**2 + / (2 * k_B * (1e9 * nu) ** 2) + ) - delta_T = np.array([-eps / 2, eps / 2]) - dI_dTCMB = ( - np.diff(planck_spectrum(nu=1e9 * nu, T=T_CMB + delta_T[:, None]), axis=0)[0] + return ( + 1e12 + * self.efficiency + * k_B + * np.diff(np.trapz(TRJ * self.passband(nu), 1e9 * nu))[0] / eps ) - dP_dTCMB = self.efficiency * np.trapz(dI_dTCMB * self.passband(nu), 1e9 * nu) - - return self.efficiency * dP_dTCMB @property def wavelength(self): @@ -321,12 +337,11 @@ def names(self): @property def summary(self) -> pd.DataFrame: - table = pd.DataFrame(columns=list(BAND_FIELDS.keys()), index=self.names) + summary = pd.DataFrame(index=self.names) - for attr, d in BAND_FIELDS.items(): - dtype = d["dtype"] - for band in self.bands: - table.at[band.name, attr] = getattr(band, attr) - table[attr] = table[attr].astype(dtype) + for band in self.bands: + band_summary = band.summary + for field, entry in FIELD_FORMATS.iterrows(): + summary.loc[band.name, field] = band_summary[field] - return table + return summary diff --git a/maria/instrument/bands/configs/abs.yml b/maria/instrument/bands/configs/abs.yml index f1b7361..7c4d535 100644 --- a/maria/instrument/bands/configs/abs.yml +++ b/maria/instrument/bands/configs/abs.yml @@ -2,6 +2,6 @@ f150: center: 150 width: 30 time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz efficiency: 0.5 diff --git a/maria/instrument/bands/configs/act.yml b/maria/instrument/bands/configs/act.yml index e99a8fc..502bf5e 100644 --- a/maria/instrument/bands/configs/act.yml +++ b/maria/instrument/bands/configs/act.yml @@ -1,51 +1,51 @@ pa4: f150: center: 150 - width: 50 + width: 30 shape: gaussian time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz efficiency: 0.5 f220: center: 220 - width: 60 + width: 40 shape: gaussian time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz efficiency: 0.5 pa5: f090: center: 90 - width: 30 + width: 20 shape: gaussian time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz efficiency: 0.5 f150: center: 150 - width: 50 + width: 30 shape: gaussian time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz efficiency: 0.5 pa6: f090: center: 90 - width: 30 + width: 20 shape: gaussian time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz efficiency: 0.5 f150: center: 150 - width: 50 + width: 30 shape: gaussian time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz efficiency: 0.5 diff --git a/maria/instrument/bands/configs/alma.yml b/maria/instrument/bands/configs/alma.yml index e226057..25bf499 100644 --- a/maria/instrument/bands/configs/alma.yml +++ b/maria/instrument/bands/configs/alma.yml @@ -3,87 +3,87 @@ f043: width: 16 shape: gaussian time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz - efficiency: 1.0 + efficiency: 0.5 f078: center: 78 width: 22 shape: gaussian time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz - efficiency: 1.0 + efficiency: 0.5 f100: center: 100 width: 32 shape: gaussian time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz - efficiency: 1.0 + efficiency: 0.5 f144: center: 144 width: 38 shape: gaussian time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz - efficiency: 1.0 + efficiency: 0.5 f187: center: 187 width: 48 shape: gaussian time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz - efficiency: 1.0 + efficiency: 0.5 f243: center: 243 width: 64 shape: gaussian time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz - efficiency: 1.0 + efficiency: 0.5 f324: center: 324 width: 98 shape: gaussian time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz - efficiency: 1.0 + efficiency: 0.5 f447: center: 447 width: 114 shape: gaussian time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz - efficiency: 1.0 + efficiency: 0.5 f661: center: 661 width: 118 shape: gaussian time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz - efficiency: 1.0 + efficiency: 0.5 f869: center: 869 width: 163 shape: gaussian time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz - efficiency: 1.0 + efficiency: 0.5 diff --git a/maria/instrument/bands/configs/atlast.yml b/maria/instrument/bands/configs/atlast.yml index 0a2de4e..b3990aa 100644 --- a/maria/instrument/bands/configs/atlast.yml +++ b/maria/instrument/bands/configs/atlast.yml @@ -3,8 +3,8 @@ f027: width: 5 shape: gaussian time_constant: 0 - efficiency: 1.0 - NEP: 1.0 # pW / sqrt(s) + efficiency: 0.5 + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz f039: @@ -12,8 +12,8 @@ f039: width: 5 shape: gaussian time_constant: 0 - efficiency: 1.0 - NEP: 1.0 # pW / sqrt(s) + efficiency: 0.5 + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz f093: @@ -21,8 +21,8 @@ f093: width: 10 shape: gaussian time_constant: 0 - efficiency: 1.0 - NEP: 1.0 # pW / sqrt(s) + efficiency: 0.5 + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz f150: @@ -30,8 +30,8 @@ f150: width: 10 shape: gaussian time_constant: 0 - efficiency: 1.0 - NEP: 1.0 # pW / sqrt(s) + efficiency: 0.5 + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz f225: @@ -39,8 +39,8 @@ f225: width: 30 shape: gaussian time_constant: 0 - efficiency: 1.0 - NEP: 1.0 # pW / sqrt(s) + efficiency: 0.5 + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz f280: @@ -48,6 +48,6 @@ f280: width: 40 shape: gaussian time_constant: 0 - efficiency: 1.0 - NEP: 1.0 # pW / sqrt(s) + efficiency: 0.5 + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz diff --git a/maria/instrument/bands/configs/mustang2.yml b/maria/instrument/bands/configs/mustang2.yml index d21c19d..fbc74a2 100644 --- a/maria/instrument/bands/configs/mustang2.yml +++ b/maria/instrument/bands/configs/mustang2.yml @@ -3,7 +3,7 @@ f093: center: 90 width: 30 time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz efficiency: 0.5 @@ -11,6 +11,6 @@ f093: f093-custom: passband: data/passbands/mustang2/f093.csv time_constant: 0 - NEP: 1.0 # pW / sqrt(s) + NEP: 3.e-5 # pW / sqrt(s) knee: 1.0 # Hz efficiency: 0.2 diff --git a/maria/instrument/bands/format.csv b/maria/instrument/bands/format.csv new file mode 100644 index 0000000..b90f361 --- /dev/null +++ b/maria/instrument/bands/format.csv @@ -0,0 +1,7 @@ +,units,dtype,format +center,GHz,float,f +width,GHz,float,f +shape,none,str,f +efficiency,none,float,f +sensitivity,K√s,float,e +NEP,pW√s,float,e diff --git a/maria/instrument/configs.yml b/maria/instrument/configs.yml index de0b0a0..c8f2abe 100644 --- a/maria/instrument/configs.yml +++ b/maria/instrument/configs.yml @@ -8,6 +8,7 @@ default: f150: center: 150 width: 30 + primary_size: 5 diff --git a/maria/map/__init__.py b/maria/map/__init__.py index d3d522a..7682d7c 100644 --- a/maria/map/__init__.py +++ b/maria/map/__init__.py @@ -5,6 +5,7 @@ import scipy as sp from tqdm import tqdm +from ..constants import k_B from ..instrument import beams from .map import Map @@ -24,46 +25,82 @@ def _run(self, **kwargs): def _sample_maps(self): dx, dy = self.coords.offsets(frame=self.map.frame, center=self.map.center) - self.data["map"] = np.zeros((dx.shape)) + self.data["map"] = 1e-16 * np.random.standard_normal(size=dx.shape) - pbar = tqdm(enumerate(self.map.frequency), disable=not self.verbose) - for i, nu in pbar: - pbar.set_description(f"Sampling input map ({nu} GHz)") + pbar = tqdm(self.instrument.bands, disable=not self.verbose) + + for band in pbar: + pbar.set_description(f"Sampling map ({band.name})") + + band_mask = self.instrument.dets.band_name == band.name + + nu = np.linspace(band.nu_min, band.nu_max, 64) + + TRJ = sp.interpolate.interp1d( + self.map.frequency, + self.map.data, + axis=0, + kind="nearest", + bounds_error=False, + fill_value="extrapolate", + )(nu) + + power_map = ( + 1e12 + * k_B + * np.trapz(band.passband(nu)[:, None, None] * TRJ, axis=0, x=1e9 * nu) + ) # nu is in GHz, f is in Hz nu_fwhm = beams.compute_angular_fwhm( - fwhm_0=self.instrument.dets.primary_size.mean(), z=np.inf, f=1e9 * nu + fwhm_0=self.instrument.dets.primary_size.mean(), + z=np.inf, + f=1e9 * band.center, ) nu_map_filter = beams.construct_beam_filter(fwhm=nu_fwhm, res=self.map.res) - filtered_nu_map_data = beams.separably_filter( - self.map.data[i], nu_map_filter - ) - - # band_res_radians = 1.22 * (299792458 / (1e9 * nu)) / self.instrument.primary_size - # band_res_pixels = band_res_radians / self.map.res - # FWHM_TO_SIGMA = 2.355 - # band_beam_sigma_pixels = band_res_pixels / FWHM_TO_SIGMA - - # # band_beam_filter = self.instrument. - - # # filtered_map_data = sp.ndimage.convolve() - det_freq_response = self.instrument.dets.passband(nu=np.array([nu]))[:, 0] - det_mask = det_freq_response > -np.inf # -1e-3 + filtered_power_map = beams.separably_filter(power_map, nu_map_filter) - samples = sp.interpolate.RegularGridInterpolator( + self.data["map"][band_mask] += sp.interpolate.RegularGridInterpolator( (self.map.x_side, self.map.y_side), - filtered_nu_map_data, + filtered_power_map, bounds_error=False, fill_value=0, method="linear", - )((dx[det_mask], dy[det_mask])) - - self.data["map"][det_mask] = samples - - if hasattr(self, "atmospheric_transmission"): - self.data["map"] *= self.atmospheric_transmission - - self.data["map"] *= 1e12 * self.instrument.dets.dP_dTRJ[:, None] + )((dx[band_mask], dy[band_mask])) + + # pbar = tqdm(enumerate(self.map.frequency), disable=not self.verbose) + # for i, nu in pbar: + # pbar.set_description(f"Sampling input map ({nu} GHz)") + + # # nu is in GHz, f is in Hz + # nu_fwhm = beams.compute_angular_fwhm( + # fwhm_0=self.instrument.dets.primary_size.mean(), z=np.inf, f=1e9 * nu + # ) + # nu_map_filter = beams.construct_beam_filter(fwhm=nu_fwhm, res=self.map.res) + # filtered_nu_map_data = beams.separably_filter( + # self.map.data[i], nu_map_filter + # ) + + # # band_res_radians = 1.22 * (299792458 / (1e9 * nu)) / self.instrument.primary_size + # # band_res_pixels = band_res_radians / self.map.res + # # FWHM_TO_SIGMA = 2.355 + # # band_beam_sigma_pixels = band_res_pixels / FWHM_TO_SIGMA + + # # # band_beam_filter = self.instrument. + + # # # filtered_map_data = sp.ndimage.convolve() + # det_freq_response = self.instrument.dets.passband(nu=np.array([nu]))[:, 0] + # det_mask = det_freq_response > -np.inf # -1e-3 + + # samples = sp.interpolate.RegularGridInterpolator( + # (self.map.x_side, self.map.y_side), + # filtered_nu_map_data, + # bounds_error=False, + # fill_value=0, + # method="linear", + # )((dx[det_mask], dy[det_mask])) + + # self.data["map"][det_mask] = self.instrument.dets.dP_dTRJ[:, None] * samples def from_fits( diff --git a/maria/noise/__init__.py b/maria/noise/__init__.py index e10c03a..aee28ab 100644 --- a/maria/noise/__init__.py +++ b/maria/noise/__init__.py @@ -1,5 +1,6 @@ import numpy as np from todder.sim.noise import generate_noise_with_knee +from tqdm import tqdm from ..sim.base import BaseSimulation @@ -7,12 +8,17 @@ class NoiseMixin: def _run(self): self._simulate_noise() - self.tod_data = self.data["noise"] def _simulate_noise(self): self.data["noise"] = np.zeros((self.instrument.n_dets, self.plan.n_time)) - for band in self.instrument.dets.bands: + bands = tqdm( + self.instrument.dets.bands, + desc="Generating noise", + disable=not self.verbose, + ) + + for band in bands: band_mask = self.instrument.dets.band_name == band.name self.data["noise"][band_mask] = generate_noise_with_knee( diff --git a/maria/sim/__init__.py b/maria/sim/__init__.py index ef82368..620665f 100644 --- a/maria/sim/__init__.py +++ b/maria/sim/__init__.py @@ -107,18 +107,14 @@ def _run(self): if hasattr(self, "cmb"): self._simulate_cmb_emission() - # convert to source Rayleigh-Jeans - # self.data["atmosphere"] *= self.instrument.dets.pW_to_KRJ[:, None] - # self.data["atmosphere"] /= self.atmospheric_transmission - if hasattr(self, "map"): self._sample_maps() - # calibrate so that there is unit efficiency to celestial sources + # scale the noise so that there is if hasattr(self, "atmospheric_transmission"): - for k in self.data: - self.data[k] /= self.atmospheric_transmission - # self.data[k] *= self.instrument.dets.pW_to_KRJ[:, None] + for k in ["cmb", "map"]: + if k in self.data: + self.data[k] *= self.atmospheric_transmission def __repr__(self): object_reprs = [