From 3f92ef077d944f8b18a5337977e2e3d1d9d8164a Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Mon, 4 Mar 2024 21:23:38 -0500 Subject: [PATCH] use power units under the hood --- docs/source/supported-sites.rst | 5 +- maria/atmosphere/__init__.py | 48 ++++++-------- maria/atmosphere/turbulent_layer.py | 2 +- maria/atmosphere/weather.py | 2 +- maria/{utils => }/constants.py | 2 + maria/{utils => }/errors.py | 0 maria/instrument/__init__.py | 4 +- maria/instrument/bands.py | 38 ++++++----- maria/instrument/bands/abs.yml | 2 +- maria/instrument/bands/act.yml | 12 ++-- maria/instrument/bands/alma.yml | 20 +++--- maria/instrument/bands/atlast.yml | 12 ++-- maria/instrument/bands/mustang2.yml | 2 +- maria/instrument/dets.py | 66 ++++++++++---------- maria/map/__init__.py | 10 +-- maria/map/mappers.py | 4 +- maria/noise/__init__.py | 25 ++++---- maria/sim/__init__.py | 28 ++++++--- maria/sim/base.py | 18 +++--- maria/sim/{default_params.yml => params.yml} | 25 +++----- maria/tests/test_pipeline.py | 2 - maria/tod/__init__.py | 48 +++----------- maria/utils/__init__.py | 1 - maria/utils/functions.py | 15 +++-- 24 files changed, 179 insertions(+), 212 deletions(-) rename maria/{utils => }/constants.py (89%) rename maria/{utils => }/errors.py (100%) rename maria/sim/{default_params.yml => params.yml} (73%) diff --git a/docs/source/supported-sites.rst b/docs/source/supported-sites.rst index 3d8a89a..5d8b87f 100644 --- a/docs/source/supported-sites.rst +++ b/docs/source/supported-sites.rst @@ -1,6 +1,7 @@ -############### + +=============== Supported sites -############### +=============== +---------------------+--------------------------------------------------------------+------------------+-------------------+--------------------+---------------------+------------------------------------------------------------------------+ | site | description | region | altitude (meters) | latitude (degrees) | longitude (degrees) | documentation | diff --git a/maria/atmosphere/__init__.py b/maria/atmosphere/__init__.py index 79ac2c0..bfdd1ab 100644 --- a/maria/atmosphere/__init__.py +++ b/maria/atmosphere/__init__.py @@ -7,6 +7,7 @@ from tqdm import tqdm from .. import utils +from ..constants import k_B from ..sim.base import BaseSimulation from ..site import InvalidRegionError, all_regions from .turbulent_layer import TurbulentLayer @@ -84,12 +85,6 @@ def __init__( from_cache=weather_from_cache, ) - def emission(self, nu): - return - - def transmission(self, nu): - return - class AtmosphereMixin: def _initialize_atmosphere(self): @@ -196,32 +191,18 @@ def _simulate_atmospheric_emission(self, units="K_RJ"): ) for band in bands: - band_index = self.instrument.dets.subset(band=band.name).index + band_index = self.instrument.dets.subset(band_name=band.name).index if self.verbose: bands.set_description(f"Computing atm. emission ({band.name})") - # multiply by 1e9 to go from GHz to Hz - # det_power_grid = ( - # 1e9 - # * 1.380649e-23 - # * np.trapz( - # self.atmosphere.emission_spectrum - # * band.passband(self.atmosphere.spectrum_side_nu), - # self.atmosphere.spectrum_side_nu, - # axis=-1, - # ) - # ) - - # this is NOT power - det_power_grid = np.sum( + # in watts. the 1e9 is for GHz -> Hz + det_power_grid = k_B * np.trapz( self.atmosphere.emission_spectrum * band.passband(self.atmosphere.spectrum_side_nu), + 1e9 * self.atmosphere.spectrum_side_nu, axis=-1, ) - det_power_grid /= np.sum( - band.passband(self.atmosphere.spectrum_side_nu), axis=-1 - ) band_power_interpolator = sp.interpolate.RegularGridInterpolator( ( @@ -240,6 +221,11 @@ def _simulate_atmospheric_emission(self, units="K_RJ"): ) ) + self.atmospheric_transmission = np.empty( + (self.instrument.n_dets, self.plan.n_time), dtype=np.float32 + ) + + # to make a new progress bar bands = ( tqdm(self.instrument.dets.bands) if self.verbose @@ -247,7 +233,7 @@ def _simulate_atmospheric_emission(self, units="K_RJ"): ) for band in bands: - band_index = self.instrument.dets.subset(band=band.name).index + band_index = self.instrument.dets.subset(band_name=band.name).index if self.verbose: bands.set_description(f"Computing atm. transmission ({band.name})") @@ -257,14 +243,13 @@ def _simulate_atmospheric_emission(self, units="K_RJ"): * self.atmosphere.spectrum_side_nu**2 ) - # multiply by 1e9 to go from GHz to Hz self.det_transmission_grid = np.trapz( rel_T_RJ_spectrum * self.atmosphere.transmission_spectrum, - self.atmosphere.spectrum_side_nu, + 1e9 * self.atmosphere.spectrum_side_nu, axis=-1, ) / np.trapz( rel_T_RJ_spectrum, - self.atmosphere.spectrum_side_nu, + 1e9 * self.atmosphere.spectrum_side_nu, axis=-1, ) @@ -277,7 +262,12 @@ def _simulate_atmospheric_emission(self, units="K_RJ"): self.det_transmission_grid, ) - self.atmospheric_transmission = band_transmission_interpolator( + # what's happening here? the atmosphere blocks some of the light from space. + # we want to calibrate to the stuff in space, so we make the atmosphere *hotter* + + self.atmospheric_transmission[ + band_index + ] = band_transmission_interpolator( ( self.zenith_scaled_pwv[band_index], self.atmosphere.weather.temperature[0], diff --git a/maria/atmosphere/turbulent_layer.py b/maria/atmosphere/turbulent_layer.py index d7dfff9..f195800 100644 --- a/maria/atmosphere/turbulent_layer.py +++ b/maria/atmosphere/turbulent_layer.py @@ -312,7 +312,7 @@ def sample(self): for band in self.instrument.bands: # we assume the atmosphere looks the same for every nu in the band - band_index = self.instrument.dets.subset(band=band.name).index + band_index = self.instrument.dets.subset(band_name=band.name).index band_angular_fwhm = self.instrument.angular_fwhm(z=self.depth)[ band_index ].mean() diff --git a/maria/atmosphere/weather.py b/maria/atmosphere/weather.py index 82099bf..a478af3 100644 --- a/maria/atmosphere/weather.py +++ b/maria/atmosphere/weather.py @@ -7,9 +7,9 @@ import pytz import scipy as sp +from ..constants import g from ..site import InvalidRegionError, all_regions, supported_regions_table from ..utils import get_utc_day_hour, get_utc_year_day -from ..utils.constants import g from . import utils here, this_filename = os.path.split(__file__) diff --git a/maria/utils/constants.py b/maria/constants.py similarity index 89% rename from maria/utils/constants.py rename to maria/constants.py index 6c705fe..05890f2 100644 --- a/maria/utils/constants.py +++ b/maria/constants.py @@ -8,3 +8,5 @@ k_B = 1.380649e-23 T_CMB = 2.7255 + +c = 299_792_458.0 diff --git a/maria/utils/errors.py b/maria/errors.py similarity index 100% rename from maria/utils/errors.py rename to maria/errors.py diff --git a/maria/instrument/__init__.py b/maria/instrument/__init__.py index c5a9f95..236d53c 100644 --- a/maria/instrument/__init__.py +++ b/maria/instrument/__init__.py @@ -107,7 +107,7 @@ def __repr__(self): @property def ubands(self): - return self.dets.ubands + return self.dets.bands.names @property def offset_x(self): @@ -178,7 +178,7 @@ def plot_dets(self, units="deg"): legend_handles = [] for iub, uband in enumerate(self.ubands): - band_mask = self.dets.band == uband + band_mask = self.dets.band_name == uband fwhm = np.degrees(self.fwhm[band_mask]) offsets = np.degrees(self.offsets[band_mask]) diff --git a/maria/instrument/bands.py b/maria/instrument/bands.py index 4410866..1bcda69 100644 --- a/maria/instrument/bands.py +++ b/maria/instrument/bands.py @@ -7,13 +7,14 @@ import pandas as pd from .. import utils +from ..constants import k_B from ..utils.io import flatten_config, read_yaml BAND_FIELD_TYPES = { "center": "float", "width": "float", "shape": "str", - "tau": "float", + "time_constant": "float", "white_noise": "float", "pink_noise": "float", "efficiency": "float", @@ -102,18 +103,16 @@ class Band: name: str center: float width: float - tau: float = 0.0 + time_constant: float = 0.0 white_noise: float = 0.0 pink_noise: float = 0.0 - shape: str = "gaussian" + shape: str = "top_hat" efficiency: float = 1.0 @classmethod def from_passband(cls, name, nu, pb, pb_err=None): center = np.round(np.sum(pb * nu), 3) - width = np.round( - nu[pb > pb.max() / np.e**2].ptp(), 3 - ) # width is the two-sigma interval + width = np.round(nu[pb > 1e-2 * pb.max()].ptp(), 3) band = cls(name=name, center=center, width=width, shape="custom") @@ -126,34 +125,43 @@ def from_passband(cls, name, nu, pb, pb_err=None): def nu_min(self) -> float: if self.shape == "flat": return self.center - 0.5 * self.width - if self.shape == "gaussian": - return self.center - self.width if self.shape == "custom": return self._nu[self._pb > 1e-2 * self._pb.max()].min() + return self.center - self.width + @property def nu_max(self) -> float: if self.shape == "flat": return self.center + 0.5 * self.width - if self.shape == "gaussian": - return self.center + self.width if self.shape == "custom": return self._nu[self._pb > 1e-2 * self._pb.max()].max() + return self.center + self.width + def passband(self, nu): """ - Passband response as a function of nu (in GHz). These integrate to one. + Passband response as a function of nu (in GHz). """ - _nu = np.atleast_1d(nu) - if self.shape == "gaussian": - band_sigma = self.width / 4 + if self.shape == "top_hat": + return np.exp(-np.log(2) * (2 * (_nu - self.center) / self.width) ** 8) - return np.exp(-0.5 * np.square((_nu - self.center) / band_sigma)) + if self.shape == "gaussian": + return np.exp(-np.log(2) * (2 * (_nu - self.center) / self.width) ** 2) if self.shape == "flat": return np.where((_nu > self.nu_min) & (_nu < self.nu_max), 1.0, 0.0) elif self.shape == "custom": return np.interp(_nu, self._nu, self._pb) + + @property + def abs_cal_rj(self): + """ + Absolute calibration (i.e. temperature per power) in Kelvins per Watt. + """ + nu = np.linspace(self.nu_min, self.nu_max, 256) + dP_dT = self.efficiency * np.trapz(self.passband(nu) * k_B, 1e9 * nu) + return 1 / dP_dT diff --git a/maria/instrument/bands/abs.yml b/maria/instrument/bands/abs.yml index 8362c66..ea791f4 100644 --- a/maria/instrument/bands/abs.yml +++ b/maria/instrument/bands/abs.yml @@ -1,7 +1,7 @@ f150: center: 150 width: 30 - tau: 0 + time_constant: 0 white_noise: 266.e-6 # in K_RJ s^-1/2 pink_noise: 1.e-1 # in K_RJ s^-1/2 efficiency: 0.5 diff --git a/maria/instrument/bands/act.yml b/maria/instrument/bands/act.yml index 72cba8b..275677a 100644 --- a/maria/instrument/bands/act.yml +++ b/maria/instrument/bands/act.yml @@ -3,7 +3,7 @@ pa4: center: 150 width: 30 shape: gaussian - tau: 0 + time_constant: 0 white_noise: 266.e-6 pink_noise: 1.e-1 efficiency: 0.5 @@ -11,7 +11,7 @@ pa4: center: 220 width: 30 shape: gaussian - tau: 0 + time_constant: 0 white_noise: 266.e-6 pink_noise: 1.e-1 efficiency: 0.5 @@ -20,7 +20,7 @@ pa5: center: 90 width: 30 shape: gaussian - tau: 0 + time_constant: 0 white_noise: 266.e-6 pink_noise: 1.e-1 efficiency: 0.5 @@ -28,7 +28,7 @@ pa5: center: 150 width: 30 shape: gaussian - tau: 0 + time_constant: 0 white_noise: 266.e-6 pink_noise: 1.e-1 efficiency: 0.5 @@ -37,7 +37,7 @@ pa6: center: 90 width: 30 shape: gaussian - tau: 0 + time_constant: 0 white_noise: 266.e-6 pink_noise: 1.e-1 efficiency: 0.5 @@ -45,7 +45,7 @@ pa6: center: 150 width: 30 shape: gaussian - tau: 0 + time_constant: 0 white_noise: 266.e-6 pink_noise: 1.e-1 efficiency: 0.5 diff --git a/maria/instrument/bands/alma.yml b/maria/instrument/bands/alma.yml index 831d8a3..7d0c850 100644 --- a/maria/instrument/bands/alma.yml +++ b/maria/instrument/bands/alma.yml @@ -2,7 +2,7 @@ f043: center: 43 width: 16 shape: gaussian - tau: 0 + time_constant: 0 white_noise: 0 pink_noise: 0 efficiency: 1.0 @@ -11,7 +11,7 @@ f078: center: 78 width: 22 shape: gaussian - tau: 0 + time_constant: 0 white_noise: 0 pink_noise: 0 efficiency: 1.0 @@ -20,7 +20,7 @@ f100: center: 100 width: 32 shape: gaussian - tau: 0 + time_constant: 0 white_noise: 0 pink_noise: 0 efficiency: 1.0 @@ -29,7 +29,7 @@ f144: center: 144 width: 38 shape: gaussian - tau: 0 + time_constant: 0 white_noise: 0 pink_noise: 0 efficiency: 1.0 @@ -38,7 +38,7 @@ f187: center: 187 width: 48 shape: gaussian - tau: 0 + time_constant: 0 white_noise: 0 pink_noise: 0 efficiency: 1.0 @@ -47,7 +47,7 @@ f243: center: 243 width: 64 shape: gaussian - tau: 0 + time_constant: 0 white_noise: 0 pink_noise: 0 efficiency: 1.0 @@ -56,7 +56,7 @@ f324: center: 324 width: 98 shape: gaussian - tau: 0 + time_constant: 0 white_noise: 0 pink_noise: 0 efficiency: 1.0 @@ -65,7 +65,7 @@ f447: center: 447 width: 114 shape: gaussian - tau: 0 + time_constant: 0 white_noise: 0 pink_noise: 0 efficiency: 1.0 @@ -74,7 +74,7 @@ f661: center: 661 width: 118 shape: gaussian - tau: 0 + time_constant: 0 white_noise: 0 pink_noise: 0 efficiency: 1.0 @@ -83,7 +83,7 @@ f869: center: 869 width: 163 shape: gaussian - tau: 0 + time_constant: 0 white_noise: 0 pink_noise: 0 efficiency: 1.0 diff --git a/maria/instrument/bands/atlast.yml b/maria/instrument/bands/atlast.yml index bc76f36..e0334b3 100644 --- a/maria/instrument/bands/atlast.yml +++ b/maria/instrument/bands/atlast.yml @@ -2,7 +2,7 @@ f027: center: 27 width: 5 shape: gaussian - tau: 0 + time_constant: 0 efficiency: 1.0 white_noise: 0 pink_noise: 0 @@ -11,7 +11,7 @@ f039: center: 39 width: 5 shape: gaussian - tau: 0 + time_constant: 0 efficiency: 1.0 white_noise: 0 pink_noise: 0 @@ -20,7 +20,7 @@ f093: center: 93 width: 10 shape: gaussian - tau: 0 + time_constant: 0 efficiency: 1.0 white_noise: 0 pink_noise: 0 @@ -29,7 +29,7 @@ f150: center: 150 width: 10 shape: gaussian - tau: 0 + time_constant: 0 efficiency: 1.0 white_noise: 0 pink_noise: 0 @@ -38,7 +38,7 @@ f225: center: 225 width: 30 shape: gaussian - tau: 0 + time_constant: 0 efficiency: 1.0 white_noise: 0 pink_noise: 0 @@ -47,7 +47,7 @@ f280: center: 280 width: 40 shape: gaussian - tau: 0 + time_constant: 0 efficiency: 1.0 white_noise: 0 pink_noise: 0 diff --git a/maria/instrument/bands/mustang2.yml b/maria/instrument/bands/mustang2.yml index 93bff67..2316107 100644 --- a/maria/instrument/bands/mustang2.yml +++ b/maria/instrument/bands/mustang2.yml @@ -2,7 +2,7 @@ f093: # file: data/mustang2/passbands/f093.csv center: 90 width: 30 - tau: 0 + time_constant: 0 white_noise: 266.e-6 # in K_RJ s^-1/2 pink_noise: 1.e-1 # in K_RJ s^-1/2 efficiency: 0.5 diff --git a/maria/instrument/dets.py b/maria/instrument/dets.py index 93aa868..066797e 100644 --- a/maria/instrument/dets.py +++ b/maria/instrument/dets.py @@ -15,12 +15,10 @@ for t in [*np.linspace(0.05, 0.95, 12)] ] -REQUIRED_DET_CONFIG_KEYS = ["n", "band_center", "band_width"] - DET_COLUMN_TYPES = { "tag": "str", "uid": "str", - "band": "str", + "band_name": "str", "band_center": "float", "offset_x": "float", "offset_y": "float", @@ -28,9 +26,14 @@ "baseline_y": "float", "baseline_z": "float", "pol_angle": "float", - "efficiency": "float", "primary_size": "float", "bath_temp": "float", + "abs_cal_rj": "float", + "time_constant": "float", + "white_noise": "float", + "pink_noise": "float", + "efficiency": "float", + "abs_cal_rj": "float", } SUPPORTED_GEOMETRIES = ["flower", "hex", "square", "circle"] @@ -114,12 +117,14 @@ def generate_dets( for band in bands: band_dets = pd.DataFrame( - index=np.arange(n), columns=["band", "offset_x", "offset_y"], dtype=float + index=np.arange(n), + columns=["band_name", "offset_x", "offset_y"], + dtype=float, ) - band_dets.loc[:, "band"] = band - band_dets.loc[:, "offset_x"] = array_offset[0] + detector_offsets[0] - band_dets.loc[:, "offset_y"] = array_offset[1] + detector_offsets[1] + band_dets.loc[:, "band_name"] = band + band_dets.loc[:, "offset_x"] = np.radians(array_offset[0] + detector_offsets[0]) + band_dets.loc[:, "offset_y"] = np.radians(array_offset[1] + detector_offsets[1]) band_dets.loc[:, "baseline_x"] = baseline_offset[0] + baselines[0] band_dets.loc[:, "baseline_y"] = baseline_offset[1] + baselines[1] @@ -230,8 +235,8 @@ def from_config(cls, config): bands = BandList.from_config(bands_config) for band in bands: - df.loc[df.band == band.name, "band_center"] = band.center - df.loc[df.band == band.name, "efficiency"] = band.efficiency + df.loc[df.band_name == band.name, "band_center"] = band.center + df.loc[df.band_name == band.name, "efficiency"] = band.efficiency return cls(df=df, bands=bands) @@ -245,18 +250,27 @@ def __init__(self, df: pd.DataFrame, bands: dict = {}): self.df = df self.bands = bands + for attr in [ + "time_constant", + "white_noise", + "pink_noise", + "efficiency", + "abs_cal_rj", + ]: + values = np.zeros(shape=self.n) + for band in self.bands: + values[self.band_name == band.name] = getattr(band, attr) + self.df.loc[:, attr] = values + def __getattr__(self, attr): if attr in self.df.columns: return self.df.loc[:, attr].values.astype(DET_COLUMN_TYPES[attr]) - # def __call__(self, band=None): - # if band is not None: - # return self.subset(band=band) + raise AttributeError(f"'Detectors' object has no attribute '{attr}'") - def subset(self, band=None): - bands = BandList(self.bands[band]) - mask = self.band == band - return Detectors(bands=bands, df=self.df.loc[mask]) + def subset(self, band_name=None): + bands = BandList([self.bands[band_name]]) + return Detectors(bands=bands, df=self.df.loc[self.band_name == band_name]) @property def n(self): @@ -266,22 +280,6 @@ def n(self): def offset(self): return np.c_[self.offset_x, self.offset_y] - # @property - # def band_center(self): - # centers = np.zeros(shape=self.n) - # for band in self.bands: - # centers[self.band == band.name] = band.center - - # return centers - - # @property - # def band_width(self): - # widths = np.zeros(shape=self.n) - # for band in self.bands: - # widths[self.band == band.name] = band.width - - # return widths - @property def __len__(self): return len(self.df) @@ -300,6 +298,6 @@ def passband(self, nu): PB = np.zeros((len(self.df), len(_nu))) for band in self.bands: - PB[self.band == band.name] = band.passband(_nu) + PB[self.band_name == band.name] = band.passband(_nu) return PB diff --git a/maria/map/__init__.py b/maria/map/__init__.py index c858f16..b01d8df 100644 --- a/maria/map/__init__.py +++ b/maria/map/__init__.py @@ -135,7 +135,7 @@ def plot(self, cmap="plasma", units="degrees", **kwargs): if units == "arcsec": map_extent = 3600 * np.degrees(map_extent_radians) - vmin, vmax = np.nanpercentile(self.data, q=[0.1, 99.9]) + vmin, vmax = np.nanpercentile(self.data, q=[1, 99]) map_im = ax.imshow( self.data.T, @@ -258,7 +258,9 @@ def _sample_maps(self): method="linear", )((dx[det_mask], dy[det_mask])) - if hasattr(self, "atmospheric_transmission"): - samples *= self.atmospheric_transmission - self.data["map"][det_mask] = samples + + if hasattr(self, "atmospheric_transmission"): + self.data["map"] *= self.atmospheric_transmission + + self.data["map"] /= self.instrument.dets.abs_cal_rj[:, None] diff --git a/maria/map/mappers.py b/maria/map/mappers.py index 82ddbc3..103399e 100644 --- a/maria/map/mappers.py +++ b/maria/map/mappers.py @@ -148,7 +148,7 @@ def __init__( def run(self): self.band = sorted( - [band for tod in self.tods for band in np.unique(tod.dets.band)] + [band for tod in self.tods for band in np.unique(tod.dets.band_name)] ) self.band_data = {} @@ -165,7 +165,7 @@ def run(self): # compute detector offsets WRT the map dx, dy = tod.coords.offsets(frame=self.frame, center=self.center) - band_mask = tod.dets.band == band + band_mask = tod.dets.band_name == band if self.calibrate: D = tod.data_calibrated.copy()[band_mask] diff --git a/maria/noise/__init__.py b/maria/noise/__init__.py index f141214..f378bbf 100644 --- a/maria/noise/__init__.py +++ b/maria/noise/__init__.py @@ -6,12 +6,13 @@ 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: - band_index = self.instrument.dets.subset(band=band.name).index + band_index = self.instrument.dets.subset(band_name=band.name).index if band.white_noise > 0: self.data["noise"][band_index] += ( @@ -20,14 +21,19 @@ def _simulate_noise(self): * np.random.standard_normal( size=(len(band_index), self.plan.n_time) ) + / band.abs_cal_rj ) if band.pink_noise > 0: for i in band_index: - self.data["noise"][i] += band.pink_noise * self._spectrum_noise( - spectrum_func=self._pink_spectrum, - size=int(self.plan.n_time), - dt=self.plan.dt, + self.data["noise"][i] += ( + band.pink_noise + * self._spectrum_noise( + spectrum_func=self._pink_spectrum, + size=int(self.plan.n_time), + dt=self.plan.dt, + ) + / band.abs_cal_rj ) def _spectrum_noise(self, spectrum_func, size, dt, amp=2.0): @@ -51,7 +57,7 @@ def _spectrum_noise(self, spectrum_func, size, dt, amp=2.0): return noise def _pink_spectrum(self, f, f_min=0, f_max=np.inf, amp=1.0): - s = (f / amp) ** -(self.pink_noise_slope) + s = (f / amp) ** -0.5 s[np.logical_or(f < f_min, f > f_max)] = 0 # apply band pass return s @@ -59,14 +65,7 @@ def _pink_spectrum(self, f, f_min=0, f_max=np.inf, amp=1.0): class NoiseSimulation(NoiseMixin, BaseSimulation): def __init__( self, - white_noise_level: float = 1e-2, - pink_noise_level: float = 1e-2, - pink_noise_slope: float = 0.5, *args, **kwargs, ): super(NoiseSimulation, self).__init__(*args, **kwargs) - - self.white_noise_level = white_noise_level - self.pink_noise_level = pink_noise_level - self.pink_noise_slope = pink_noise_slope diff --git a/maria/sim/__init__.py b/maria/sim/__init__.py index eb61e98..41cbbe1 100644 --- a/maria/sim/__init__.py +++ b/maria/sim/__init__.py @@ -40,6 +40,8 @@ def __init__( **self.parsed_sim_kwargs["site"], ) + self.noise = True + self.params = {} for sub_type, sub_master_params in master_params.items(): @@ -72,25 +74,35 @@ def __init__( self._initialize_atmosphere() - def _run(self, units="K_RJ"): + def _run(self): # number of bands are lost here - self._simulate_noise() + if self.noise: + self._simulate_noise() if self.atmosphere_model: self._simulate_atmospheric_emission() + # convert to source Rayleigh-Jeans + # self.data["atmosphere"] *= self.instrument.dets.abs_cal_rj[:, None] + # self.data["atmosphere"] /= self.atmospheric_transmission + if self.map_file: self._sample_maps() - if hasattr(self, "atmospheric_transmission"): - self.data["map"] *= self.atmospheric_transmission + if hasattr(self, "cmb"): + self._simulate_cmb_emission() + + self.tod_data = sum( + [self.data.get(k, 0) for k in ["atmosphere", "cmb", "map", "noise"]] + ) + + if hasattr(self, "atmospheric_transmission"): + self.tod_data /= self.atmospheric_transmission - if hasattr(self, "cmb_sim"): - self.cmb_sim._run() - self.data += self.cmb_sim.data + self.tod_data *= self.instrument.dets.abs_cal_rj[:, None] def __repr__(self): object_reprs = [ getattr(self, attr).__repr__() for attr in ["instrument", "site", "plan"] ] - return "\n\n".join(object_reprs) + return ", ".join(object_reprs) diff --git a/maria/sim/base.py b/maria/sim/base.py index 850579c..30e6579 100644 --- a/maria/sim/base.py +++ b/maria/sim/base.py @@ -19,7 +19,7 @@ def __init__(self, invalid_keys): ) -master_params = utils.io.read_yaml(f"{here}/default_params.yml") +master_params = utils.io.read_yaml(f"{here}/params.yml") def parse_sim_kwargs(kwargs, master_kwargs, strict=False): @@ -62,8 +62,6 @@ def __init__( self.verbose = verbose - self.data = {} - parsed_sim_kwargs = parse_sim_kwargs(kwargs, master_params) if isinstance(instrument, Instrument): @@ -87,6 +85,9 @@ def __init__( "The passed site must be either a Site object or a string." ) + self.data = {} + self.calibration = np.ones((self.instrument.dets.n, self.plan.n_time)) + self.boresight = Coordinates( self.plan.time, self.plan.phi, @@ -127,18 +128,15 @@ def _run(self): raise NotImplementedError() def run(self): - self._run() + self.data = {} - abscal = 1.0 - if hasattr(self, "atmospheric_transmission"): - abscal /= self.atmospheric_transmission.mean() + # Simulate all the junk + self._run() tod = TOD( - data=self.data, + data=self.tod_data, dets=self.instrument.dets.df, - # boresight=self.boresight, coords=self.det_coords, - abscal=abscal, ) return tod diff --git a/maria/sim/default_params.yml b/maria/sim/params.yml similarity index 73% rename from maria/sim/default_params.yml rename to maria/sim/params.yml index 0e713b8..f248276 100644 --- a/maria/sim/default_params.yml +++ b/maria/sim/params.yml @@ -2,25 +2,16 @@ # repeated parameters like "description" or "documentation" are not preserved in the simulation class. instrument: # defaults to a small test instrument - instrument_description: - dets: - f090: - n: 60 - bands: ["act/pa5/f090"] - f150: - n: 60 - bands: ["act/pa5/f150"] - field_of_view: 0.8 - baseline: 0 # meters - shape: hex + instrument_description: str + instrument_documentation: str primary_size: 6 + dets: dict az_bounds: [0, 360] el_bounds: [20, 90] max_az_vel: 3 max_el_vel: 2 max_az_acc: 1 max_el_acc: 0.25 - instrument_documentation: '' plan: # defaults to a 45 degree stare due north plan_description: '' @@ -31,7 +22,7 @@ plan: # defaults to a 45 degree stare due north sample_rate: 20 scan_pattern: daisy scan_center: [10, 4.5] - scan_options: {} + scan_options: mapping site: # default to the ALMA site site_description: '' @@ -52,6 +43,9 @@ atmosphere: pwv_rms_frac: 0.03 pwv: +noise: + noise: True + map: map_file: '' map_frame: ra_dec @@ -60,8 +54,3 @@ map: map_inbright: map_units: K_RJ map_freqs: [150] - -noise: - white_noise_level: 1.e-2 # in Kelvin Rayleigh-Jeans equivalent - pink_noise_level: 1.e-2 # in Kelvin Rayleigh-Jeans equivalent amplitude in fourier domain - pink_noise_slope: 0.5 diff --git a/maria/tests/test_pipeline.py b/maria/tests/test_pipeline.py index 949fe2e..8e4c9e1 100644 --- a/maria/tests/test_pipeline.py +++ b/maria/tests/test_pipeline.py @@ -40,8 +40,6 @@ def test_mustang2(): site="green_bank", # Site plan="daisy", # Scanning strategy atmosphere_model=atm_model, # atmospheric model - white_noise_level=white_noise_level, # white noise level - pink_noise_level=pink_noise_level, # pink noise level # True sky input # --------------------- map_file=inputfile, # Input files must be a fits file. diff --git a/maria/tod/__init__.py b/maria/tod/__init__.py index e29beb2..f4ee3b9 100644 --- a/maria/tod/__init__.py +++ b/maria/tod/__init__.py @@ -1,5 +1,4 @@ import json -from functools import lru_cache import h5py import numpy as np @@ -14,32 +13,9 @@ class TOD: """ - Time-ordered data. + Time-ordered data. This has per-detector pointing and data. """ - @property - @lru_cache(maxsize=None) - def data(self): - """ - Combine all the fields into the total. - """ - return sum(self._data.values()) # * self.abscal - - @property - @lru_cache(maxsize=None) - def data_calibrated(self): - """ - Combine all the fields into the total. - """ - T_ant = np.zeros(self.data.shape) - if "atmosphere" in self._data.keys(): - T_ant += 0.8 * self._data["atmosphere"] - if "map" in self._data.keys(): - T_ant += 0.3 * self._data["map"] - if "noise" in self._data.keys(): - T_ant += 1.0 * self._data["noise"] - return T_ant / 0.3 * self.abscal - @staticmethod def from_fits(fname: str, format: str, **kwargs): if format.lower() == "mustang-2": @@ -108,7 +84,7 @@ def _from_mustang2(cls, fname: str, hdu: int = 1): n_dets = len(det_uids) n_samp = det_counts.max() - data = {"data": raw["FNU"].astype("float32").reshape((n_dets, n_samp))} + data = raw["FNU"].astype("float32").reshape((n_dets, n_samp)) ra = raw["dx"].astype(float).reshape((n_dets, n_samp)) dec = raw["dy"].astype(float).reshape((n_dets, n_samp)) @@ -132,7 +108,7 @@ def _from_mustang2(cls, fname: str, hdu: int = 1): def __init__( self, coords: Coordinates, - data: dict = {}, + data: np.array, units: dict = {}, dets: pd.DataFrame = None, boresight: Coordinates = None, @@ -140,7 +116,7 @@ def __init__( ): self.coords = coords self.dets = dets - self._data = data + self.data = data self.header = fits.header.Header() self.abscal = abscal @@ -179,7 +155,7 @@ def nt(self): def subset(self, det_mask=None, time_mask=None, band: str = None): if band is not None: - det_mask = self.dets.band == band + det_mask = self.dets.band_name == band if not det_mask.sum() > 0: raise ValueError(f"There are no detectors for band '{band}'.") return self.subset(det_mask=det_mask) @@ -196,12 +172,8 @@ def subset(self, det_mask=None, time_mask=None, band: str = None): frame="az_el", ) - subset_data = {} - for k, v in self._data.items(): - subset_data[k] = v[:, time_mask] - return TOD( - data=subset_data, + data=self.data[:, time_mask], coords=subset_coords, dets=self.dets, units=self.units, @@ -221,12 +193,8 @@ def subset(self, det_mask=None, time_mask=None, band: str = None): frame="az_el", ) - subset_data = {} - for k, v in self._data.items(): - subset_data[k] = v[det_mask] - return TOD( - data=subset_data, + data=self.data[det_mask], coords=subset_coords, dets=subset_dets, units=self.units, @@ -332,7 +300,7 @@ def to_fits(self, fname, format="MUSTANG-2"): col03 = fits.Column( name="FNU ", format="E", - array=self._data["data"].flatten(), + array=self.data.flatten(), unit=self.units["data"], ) col04 = fits.Column(name="UFNU ", format="E") diff --git a/maria/utils/__init__.py b/maria/utils/__init__.py index 2a23d98..d3d1ce9 100644 --- a/maria/utils/__init__.py +++ b/maria/utils/__init__.py @@ -5,7 +5,6 @@ import numpy as np import pytz -from . import constants # noqa F401 from . import functions # noqa F401 from . import io # noqa F401 from . import linalg # noqa F401 diff --git a/maria/utils/functions.py b/maria/utils/functions.py index 994e77a..8b0042f 100644 --- a/maria/utils/functions.py +++ b/maria/utils/functions.py @@ -29,16 +29,19 @@ def normalized_matern(r, nu): ) -def approximate_normalized_matern(r, nu, n_test_points=1024): +def approximate_normalized_matern(r, nu, n_test_points=256): """ Computing BesselK[nu,z] for arbitrary nu is expensive. This is good for casting over huge matrices. """ r_safe = np.atleast_1d(np.abs(r)) - r_min = np.maximum(r_safe[r_safe > 0].min(), 1e-6) - r_max = np.minimum(r_safe.max(), 1e3) + r_min = np.maximum(1e-6, 0.5 * r_safe[r_safe > 0].min()) + r_max = np.minimum(1e3, 2.0 * r_safe.max()) - test_values = np.r_[0, np.geomspace(r_min, r_max, n_test_points - 1)] - data_values = normalized_matern(test_values, nu) * np.exp(test_values) + r_samples = np.geomspace(r_min, r_max, n_test_points) - return np.interp(r_safe, test_values, data_values) * np.exp(-r_safe) + sf_samples = 1 - normalized_matern(r_samples, nu=nu) + + sf = np.exp(np.interp(np.log(r), np.log(r_samples), np.log(sf_samples))) + + return 1 - sf