From 9c7e0e1c68c0ffcf35a197ca7a8bdfab3a70a694 Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Mon, 20 May 2024 20:38:46 +0200 Subject: [PATCH 1/2] a bunch of stuff good lord --- .cookiecutter-maria.yaml | 2 +- AUTHORS.rst | 2 +- docs/source/instruments.rst | 9 +- maria/__init__.py | 1 - maria/atmosphere/__init__.py | 293 +---------- maria/atmosphere/sim.py | 221 +++++++++ maria/atmosphere/spectrum.py | 78 +++ maria/atmosphere/turbulent_layer.py | 6 +- maria/atmosphere/weather.py | 75 ++- maria/cmb/__init__.py | 64 +-- maria/coords/__init__.py | 463 ------------------ maria/errors.py | 4 + maria/{utils => }/functions.py | 6 +- maria/instrument/__init__.py | 48 +- maria/instrument/bands.py | 226 --------- maria/instrument/bands/__init__.py | 328 +++++++++++++ maria/instrument/bands/configs/abs.yml | 7 + .../{data/bands => bands/configs}/act.yml | 36 +- .../{data/bands => bands/configs}/alma.yml | 60 +-- .../{data/bands => bands/configs}/atlast.yml | 36 +- .../bands => bands/configs}/mustang2.yml | 8 +- .../passbands/mustang2/f093.csv | 0 .../{instruments.yml => configs.yml} | 5 +- maria/instrument/data/bands/abs.yml | 7 - .../{detectors.py => detectors/__init__.py} | 3 +- maria/instrument/{ => detectors}/arrays.py | 35 +- .../{data => detectors}/arrays/alma/aca.cfg | 0 .../arrays/alma/alma.cycle1.csv | 0 .../arrays/alma/alma.cycle1.f043.csv | 0 .../arrays/alma/alma.cycle1.f078.csv | 0 .../arrays/alma/alma.cycle1.f100.csv | 0 .../arrays/alma/alma.cycle1.f144.csv | 0 .../arrays/alma/alma.cycle1.f187.csv | 0 .../arrays/alma/alma.cycle1.f243.csv | 0 .../arrays/alma/alma.cycle1.f324.csv | 0 .../arrays/alma/alma.cycle1.f447.csv | 0 .../arrays/alma/alma.cycle1.f661.csv | 0 .../arrays/alma/alma.cycle1.f869.csv | 0 .../arrays/alma/alma.cycle1.total.csv | 0 .../arrays/alma/alma.cycle10.2.cfg | 0 .../arrays/alma/alma.cycle10.3.cfg | 0 .../arrays/alma/alma.cycle10.4.cfg | 0 .../arrays/alma/alma.cycle10.5.cfg | 0 .../arrays/alma/alma.cycle10.6.cfg | 0 .../arrays/alma/alma.cycle10.7.cfg | 0 .../arrays/alma/alma.cycle10.8.cfg | 0 maria/{utils => }/io.py | 46 +- maria/map/__init__.py | 287 +++-------- maria/map/map.py | 191 ++++++++ maria/map/mappers.py | 18 +- maria/noise/__init__.py | 91 ++-- maria/plan/__init__.py | 23 +- maria/plan/plans.yml | 8 +- maria/sim/__init__.py | 112 +++-- maria/sim/base.py | 22 +- maria/sim/params.yml | 2 +- maria/site/__init__.py | 14 +- maria/site/sites.yml | 32 +- maria/tests/test_atmosphere.py | 15 +- maria/tests/test_cmb.py | 10 + maria/tests/test_coordinates.py | 3 +- maria/tests/test_noise.py | 33 ++ maria/tests/test_noise_sims.py | 14 - maria/tests/test_pipeline.py | 54 +- maria/tests/test_tod.py | 20 +- maria/tests/test_weather.py | 15 - maria/tod/__init__.py | 268 ---------- maria/tod/atlast.py | 48 -- maria/tod/mustang2.py | 40 -- maria/{utils => }/units.py | 0 maria/utils/__init__.py | 162 ------ maria/utils/weather.py | 34 -- pyproject.toml | 5 +- pytest.ini | 5 +- 74 files changed, 1403 insertions(+), 2192 deletions(-) create mode 100644 maria/atmosphere/sim.py create mode 100644 maria/atmosphere/spectrum.py delete mode 100644 maria/coords/__init__.py rename maria/{utils => }/functions.py (91%) delete mode 100644 maria/instrument/bands.py create mode 100644 maria/instrument/bands/__init__.py create mode 100644 maria/instrument/bands/configs/abs.yml rename maria/instrument/{data/bands => bands/configs}/act.yml (54%) rename maria/instrument/{data/bands => bands/configs}/alma.yml (53%) rename maria/instrument/{data/bands => bands/configs}/atlast.yml (53%) rename maria/instrument/{data/bands => bands/configs}/mustang2.yml (57%) rename maria/instrument/{data => bands}/passbands/mustang2/f093.csv (100%) rename maria/instrument/{instruments.yml => configs.yml} (95%) delete mode 100644 maria/instrument/data/bands/abs.yml rename maria/instrument/{detectors.py => detectors/__init__.py} (97%) rename maria/instrument/{ => detectors}/arrays.py (87%) rename maria/instrument/{data => detectors}/arrays/alma/aca.cfg (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle1.csv (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle1.f043.csv (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle1.f078.csv (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle1.f100.csv (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle1.f144.csv (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle1.f187.csv (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle1.f243.csv (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle1.f324.csv (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle1.f447.csv (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle1.f661.csv (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle1.f869.csv (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle1.total.csv (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle10.2.cfg (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle10.3.cfg (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle10.4.cfg (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle10.5.cfg (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle10.6.cfg (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle10.7.cfg (100%) rename maria/instrument/{data => detectors}/arrays/alma/alma.cycle10.8.cfg (100%) rename maria/{utils => }/io.py (66%) create mode 100644 maria/map/map.py create mode 100644 maria/tests/test_cmb.py create mode 100644 maria/tests/test_noise.py delete mode 100644 maria/tests/test_noise_sims.py delete mode 100644 maria/tests/test_weather.py delete mode 100644 maria/tod/__init__.py delete mode 100644 maria/tod/atlast.py delete mode 100644 maria/tod/mustang2.py rename maria/{utils => }/units.py (100%) delete mode 100644 maria/utils/weather.py diff --git a/.cookiecutter-maria.yaml b/.cookiecutter-maria.yaml index d21c8fa..9f89568 100644 --- a/.cookiecutter-maria.yaml +++ b/.cookiecutter-maria.yaml @@ -1,6 +1,6 @@ default_context: full_name: "Thomas Morris" - email: "thomasmorris@princeton.edu" + email: "thomas.w.morris@yale.edu" github_username: "thomaswmorris" project_name: "maria" package_dist_name: "maria" diff --git a/AUTHORS.rst b/AUTHORS.rst index 5f7d223..47589b5 100644 --- a/AUTHORS.rst +++ b/AUTHORS.rst @@ -5,7 +5,7 @@ Credits Maintainer ---------- -* Thomas Morris +* Thomas Morris Contributors ------------ diff --git a/docs/source/instruments.rst b/docs/source/instruments.rst index 16180d8..be71dd0 100644 --- a/docs/source/instruments.rst +++ b/docs/source/instruments.rst @@ -31,15 +31,12 @@ Custom arrays of detectors are a bit more complicated. For example: .. code-block:: python - f090 = {"center": 90, "width": 30} - f150 = {"center": 150, "width": 30} - - dets = {"n": 500, + array = {"n": 500, "field_of_view": 2 "array_shape": "hex", - "bands": [f090, f150]} + "bands": {"center": 150, "width": 30}} - my_custom_array = maria.get_instrument(dets=dets) + my_custom_array = maria.get_instrument(array=array) Actually, there are several valid ways to define an array of detectors. These are all valid: diff --git a/maria/__init__.py b/maria/__init__.py index e1edca3..f0e4d80 100644 --- a/maria/__init__.py +++ b/maria/__init__.py @@ -9,4 +9,3 @@ from .plan import all_plans, get_plan # noqa F401 from .sim import Simulation # noqa F401 from .site import all_regions, all_sites, get_site # noqa F401 -from .tod import TOD # noqa F401 diff --git a/maria/atmosphere/__init__.py b/maria/atmosphere/__init__.py index a4dcab5..59850b6 100644 --- a/maria/atmosphere/__init__.py +++ b/maria/atmosphere/__init__.py @@ -1,23 +1,13 @@ import os from datetime import datetime -import h5py -import numpy as np -import scipy as sp -from tqdm import tqdm - -from .. import utils -from ..constants import k_B -from ..plan import validate_pointing -from ..sim.base import BaseSimulation -from ..site import InvalidRegionError, all_regions -from .turbulent_layer import TurbulentLayer +from .spectrum import Spectrum from .weather import Weather # noqa F401 here, this_filename = os.path.split(__file__) SPECTRA_DATA_DIRECTORY = f"{here}/data" -SPECTRA_DATA_CACHE_DIRECTORY = "/tmp/maria/atmosphere/spectra" +SPECTRA_DATA_CACHE_DIRECTORY = "/tmp/maria-data/atmosphere/spectra" SPECTRA_DATA_URL_BASE = "https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/spectra" # noqa F401 CACHE_MAX_AGE_SECONDS = 30 * 86400 @@ -29,290 +19,25 @@ def __init__( region: str = "princeton", altitude: float = None, weather_quantiles: dict = {}, - weather_override: dict = {}, + weather_kwargs: dict = {}, weather_source: str = "era5", - weather_from_cache: bool = None, spectrum_source: str = "am", - spectrum_from_cache: bool = None, + pwv_rms_frac: float = 0.03, ): - if region not in all_regions: - raise InvalidRegionError(region) - t = t or datetime.utcnow().timestamp() - self.region = region - self.spectrum_source = spectrum_source - self.spectrum_source_path = ( - f"{SPECTRA_DATA_DIRECTORY}/{self.spectrum_source}/{self.region}.h5" - ) - - # if the data isn't in the module, default to use the cache - self.spectrum_from_cache = ( - spectrum_from_cache - if spectrum_from_cache is not None - else not os.path.exists(self.spectrum_source_path) + self.spectrum = Spectrum( + region=region, + source=spectrum_source, ) - if self.spectrum_from_cache: - self.spectrum_source_path = f"{SPECTRA_DATA_CACHE_DIRECTORY}/{self.spectrum_source}/{self.region}.h5" - utils.io.fetch_cache( - source_url=f"{SPECTRA_DATA_URL_BASE}/{self.spectrum_source}/{self.region}.h5", - cache_path=self.spectrum_source_path, - CACHE_MAX_AGE=CACHE_MAX_AGE_SECONDS, - ) - self.spectrum_from_cache = True - - with h5py.File(self.spectrum_source_path, "r") as f: - self.spectrum_side_nu = f["side_nu_GHz"][:] - self.spectrum_side_elevation = f["side_elevation_deg"][:] - self.spectrum_side_zenith_pwv = f["side_zenith_pwv_mm"][:] - self.spectrum_side_base_temperature = f["side_base_temperature_K"][:] - - self.emission_spectrum = f["emission_temperature_rayleigh_jeans_K"][:] - self.transmission_spectrum = np.exp(-f["opacity_nepers"][:]) - self.excess_path_spectrum = 1e6 * ( - f["excess_path"][:] + f["offset_excess_path_m"][:] - ) - self.weather = Weather( t=t, region=region, altitude=altitude, quantiles=weather_quantiles, - override=weather_override, + override=weather_kwargs, source=weather_source, - from_cache=weather_from_cache, - ) - - @property - def emission_interpolator(self): - return sp.interpolate.RegularGridInterpolator( - points=( - self.spectrum_side_zenith_pwv, - self.spectrum_side_base_temperature, - self.spectrum_side_elevation, - self.spectrum_side_nu, - ), - values=self.emission_spectrum, - ) - - def emission(self, nu, zenith_pwv=None, base_temperature=None, elevation=45): - if zenith_pwv is None: - zenith_pwv = np.median(self.spectrum_side_zenith_pwv) - if base_temperature is None: - base_temperature = np.median(self.spectrum_side_base_temperature) - - return self.emission_interpolator((zenith_pwv, base_temperature, elevation, nu)) - - -class AtmosphereMixin: - def _initialize_atmosphere(self): - """ - This assume that BaseSimulation.__init__() has been called. - """ - - validate_pointing(self.coords.az, self.coords.el) - - if self.atmosphere_model == "2d": - self.turbulent_layer_depths = np.linspace( - self.min_atmosphere_height, - self.max_atmosphere_height, - self.n_atmosphere_layers, - ) - self.atmosphere.layers = [] - - depths = enumerate(self.turbulent_layer_depths) - if self.verbose: - depths = tqdm(depths, desc="Initializing atmospheric layers") - - for _, layer_depth in depths: - layer_res = ( - self.instrument.physical_fwhm(z=layer_depth).min() - / self.min_atmosphere_beam_res - ) # in meters - - layer = TurbulentLayer( - instrument=self.instrument, - boresight=self.boresight, - weather=self.atmosphere.weather, - depth=layer_depth, - res=layer_res, - turbulent_outer_scale=self.turbulent_outer_scale, - ) - - self.atmosphere.layers.append(layer) - - if self.atmosphere_model == "3d": - self.initialize_3d_atmosphere() - - def _simulate_atmospheric_fluctuations(self): - if self.atmosphere_model == "2d": - self._simulate_2d_atmospheric_fluctuations() - - if self.atmosphere_model == "3d": - self._simulate_3d_atmospheric_fluctuations() - - def _simulate_2d_turbulence(self): - """ - Simulate layers of two-dimensional turbulence. - """ - - layer_data = np.zeros( - (self.n_atmosphere_layers, self.instrument.n_dets, self.plan.n_time) - ) - - layers = ( - tqdm(self.atmosphere.layers) if self.verbose else self.atmosphere.layers ) - for layer_index, layer in enumerate(layers): - layer.generate() - layer_data[layer_index] = sp.interpolate.interp1d( - layer.sim_time, - layer.sample(), - axis=-1, - kind="cubic", - bounds_error=False, - fill_value="extrapolate", - )(self.boresight.time) - - if self.verbose: - layers.set_description(f"Generating atmosphere (z={layer.depth:.00f}m)") - - return layer_data - - def _simulate_2d_atmospheric_fluctuations(self): - """ - Simulate layers of two-dimensional turbulence. - """ - - turbulence = self._simulate_2d_turbulence() - - rel_layer_scaling = sp.interpolate.interp1d( - self.atmosphere.weather.altitude_levels, - self.atmosphere.weather.absolute_humidity, - kind="linear", - )( - self.site.altitude - + self.turbulent_layer_depths[:, None, None] * np.sin(self.coords.el) - ) - rel_layer_scaling /= np.sqrt(np.square(rel_layer_scaling).sum(axis=0)[None]) - - self.layer_scaling = ( - self.pwv_rms_frac * self.atmosphere.weather.pwv * rel_layer_scaling - ) - - self.zenith_scaled_pwv = self.atmosphere.weather.pwv + ( - self.layer_scaling * turbulence - ).sum(axis=0) - - def _simulate_atmospheric_emission(self, units="K_RJ"): - if units == "K_RJ": # Kelvin Rayleigh-Jeans - self._simulate_atmospheric_fluctuations() - self.data["atmosphere"] = np.empty( - (self.instrument.n_dets, self.plan.n_time), dtype=np.float32 - ) - - bands = ( - tqdm(self.instrument.dets.bands) - if self.verbose - else self.instrument.dets.bands - ) - - for band in bands: - band_index = self.instrument.dets.subset(band_name=band.name).index - - if self.verbose: - bands.set_description(f"Computing atm. emission ({band.name})") - - # in picowatts. the 1e9 is for GHz -> Hz - det_power_grid = ( - 1e12 - * k_B - * np.trapz( - self.atmosphere.emission_spectrum - * band.passband(self.atmosphere.spectrum_side_nu), - 1e9 * self.atmosphere.spectrum_side_nu, - axis=-1, - ) - ) - - band_power_interpolator = sp.interpolate.RegularGridInterpolator( - ( - self.atmosphere.spectrum_side_zenith_pwv, - self.atmosphere.spectrum_side_base_temperature, - self.atmosphere.spectrum_side_elevation, - ), - det_power_grid, - ) - - self.data["atmosphere"][band_index] = band_power_interpolator( - ( - self.zenith_scaled_pwv[band_index], - self.atmosphere.weather.temperature[0], - np.degrees(self.coords.el[band_index]), - ) - ) - - 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 - else self.instrument.dets.bands - ) - - for band in bands: - band_index = self.instrument.dets.subset(band_name=band.name).index - - if self.verbose: - bands.set_description(f"Computing atm. transmission ({band.name})") - - rel_T_RJ_spectrum = ( - band.passband(self.atmosphere.spectrum_side_nu) - * self.atmosphere.spectrum_side_nu**2 - ) - - self.det_transmission_grid = np.trapz( - rel_T_RJ_spectrum * self.atmosphere.transmission_spectrum, - 1e9 * self.atmosphere.spectrum_side_nu, - axis=-1, - ) / np.trapz( - rel_T_RJ_spectrum, - 1e9 * self.atmosphere.spectrum_side_nu, - axis=-1, - ) - - band_transmission_interpolator = sp.interpolate.RegularGridInterpolator( - ( - self.atmosphere.spectrum_side_zenith_pwv, - self.atmosphere.spectrum_side_base_temperature, - self.atmosphere.spectrum_side_elevation, - ), - self.det_transmission_grid, - ) - - # 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], - np.degrees(self.coords.el[band_index]), - ) - ) - - if units == "F_RJ": # Fahrenheit Rayleigh-Jeans 🇺🇸 - self._simulate_atmospheric_emission(self, units="K_RJ") - self.data["atmosphere"] = 1.8 * (self.data["atmosphere"] - 273.15) + 32 - -class AtmosphereSimulation(BaseSimulation, AtmosphereMixin): - def __init__(self, *args, **kwargs): - super(AtmosphereSimulation, self).__init__(*args, **kwargs) - self._initialize_atmosphere() + self.pwv_rms_frac = pwv_rms_frac diff --git a/maria/atmosphere/sim.py b/maria/atmosphere/sim.py new file mode 100644 index 0000000..ca7b906 --- /dev/null +++ b/maria/atmosphere/sim.py @@ -0,0 +1,221 @@ +import os + +import numpy as np +import scipy as sp +from tqdm import tqdm + +from ..constants import k_B +from ..plan import validate_pointing +from .turbulent_layer import TurbulentLayer + +here, this_filename = os.path.split(__file__) + + +class AtmosphereMixin: + def _initialize_2d_atmosphere( + self, + min_atmosphere_height=500, + max_atmosphere_height=5000, + n_atmosphere_layers=4, + min_atmosphere_beam_res=4, + turbulent_outer_scale=500, + ): + """ + This assume that BaseSimulation.__init__() has been called. + """ + + validate_pointing(self.coords.az, self.coords.el) + + self.atmosphere_model = "2d" + + self.turbulent_layer_depths = np.linspace( + min_atmosphere_height, + max_atmosphere_height, + n_atmosphere_layers, + ) + self.atmosphere.layers = [] + + depths = enumerate(self.turbulent_layer_depths) + if self.verbose: + depths = tqdm(depths, desc="Initializing atmospheric layers") + + for _, layer_depth in depths: + layer_res = ( + self.instrument.physical_fwhm(z=layer_depth).min() + / min_atmosphere_beam_res + ) # in meters + + layer = TurbulentLayer( + instrument=self.instrument, + boresight=self.boresight, + weather=self.atmosphere.weather, + depth=layer_depth, + res=layer_res, + turbulent_outer_scale=turbulent_outer_scale, + ) + + self.atmosphere.layers.append(layer) + + def _simulate_atmospheric_fluctuations(self): + if self.atmosphere_model == "2d": + self._simulate_2d_atmospheric_fluctuations() + + if self.atmosphere_model == "3d": + self._simulate_3d_atmospheric_fluctuations() + + def _simulate_2d_turbulence(self): + """ + Simulate layers of two-dimensional turbulence. + """ + + layer_data = np.zeros( + (len(self.atmosphere.layers), self.instrument.n_dets, self.plan.n_time) + ) + + pbar = tqdm(enumerate(self.atmosphere.layers), disable=not self.verbose) + for layer_index, layer in pbar: + layer.generate() + layer_data[layer_index] = sp.interpolate.interp1d( + layer.sim_time, + layer.sample(), + axis=-1, + kind="cubic", + bounds_error=False, + fill_value="extrapolate", + )(self.boresight.time) + pbar.set_description(f"Generating atmosphere (z={layer.depth:.00f}m)") + + return layer_data + + def _simulate_2d_atmospheric_fluctuations(self): + """ + Simulate layers of two-dimensional turbulence. + """ + + turbulence = self._simulate_2d_turbulence() + + rel_layer_scaling = sp.interpolate.interp1d( + self.atmosphere.weather.altitude_levels, + self.atmosphere.weather.absolute_humidity, + kind="linear", + )( + self.site.altitude + + self.turbulent_layer_depths[:, None, None] * np.sin(self.coords.el) + ) + rel_layer_scaling /= np.sqrt(np.square(rel_layer_scaling).sum(axis=0)[None]) + + self.layer_scaling = ( + self.atmosphere.pwv_rms_frac + * self.atmosphere.weather.pwv + * rel_layer_scaling + ) + + self.zenith_scaled_pwv = self.atmosphere.weather.pwv + ( + self.layer_scaling * turbulence + ).sum(axis=0) + + def _simulate_atmospheric_emission(self, units="K_RJ"): + if units == "K_RJ": # Kelvin Rayleigh-Jeans + self._simulate_atmospheric_fluctuations() + self.data["atmosphere"] = np.empty( + (self.instrument.n_dets, self.plan.n_time), dtype=np.float32 + ) + + bands = ( + tqdm(self.instrument.dets.bands) + if self.verbose + else self.instrument.dets.bands + ) + + for band in bands: + band_index = self.instrument.dets.subset(band_name=band.name).index + + if self.verbose: + bands.set_description(f"Computing atm. emission ({band.name})") + + # in picowatts. the 1e9 is for GHz -> Hz + det_power_grid = ( + 1e12 + * k_B + * np.trapz( + self.atmosphere.spectrum._emission + * band.passband(self.atmosphere.spectrum._side_nu), + 1e9 * self.atmosphere.spectrum._side_nu, + axis=-1, + ) + ) + + band_power_interpolator = sp.interpolate.RegularGridInterpolator( + ( + self.atmosphere.spectrum._side_zenith_pwv, + self.atmosphere.spectrum._side_base_temperature, + self.atmosphere.spectrum._side_elevation, + ), + det_power_grid, + ) + + self.data["atmosphere"][band_index] = band_power_interpolator( + ( + self.zenith_scaled_pwv[band_index], + self.atmosphere.weather.temperature[0], + np.degrees(self.coords.el[band_index]), + ) + ) + + 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 + else self.instrument.dets.bands + ) + + for band in bands: + band_index = self.instrument.dets.subset(band_name=band.name).index + + if self.verbose: + bands.set_description(f"Computing atm. transmission ({band.name})") + + rel_T_RJ_spectrum = ( + band.passband(self.atmosphere.spectrum._side_nu) + * self.atmosphere.spectrum._side_nu**2 + ) + + self.det_transmission_grid = np.trapz( + rel_T_RJ_spectrum * self.atmosphere.spectrum._transmission, + 1e9 * self.atmosphere.spectrum._side_nu, + axis=-1, + ) / np.trapz( + rel_T_RJ_spectrum, + 1e9 * self.atmosphere.spectrum._side_nu, + axis=-1, + ) + + band_transmission_interpolator = sp.interpolate.RegularGridInterpolator( + ( + self.atmosphere.spectrum._side_zenith_pwv, + self.atmosphere.spectrum._side_base_temperature, + self.atmosphere.spectrum._side_elevation, + ), + self.det_transmission_grid, + ) + + # 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], + np.degrees(self.coords.el[band_index]), + ) + ) + + if units == "F_RJ": # Fahrenheit Rayleigh-Jeans 🇺🇸 + self._simulate_atmospheric_emission(self, units="K_RJ") + self.data["atmosphere"] = 1.8 * (self.data["atmosphere"] - 273.15) + 32 diff --git a/maria/atmosphere/spectrum.py b/maria/atmosphere/spectrum.py new file mode 100644 index 0000000..a5d8ae0 --- /dev/null +++ b/maria/atmosphere/spectrum.py @@ -0,0 +1,78 @@ +import os + +import h5py +import numpy as np +import scipy as sp + +from ..io import fetch_cache +from ..site import InvalidRegionError, all_regions + +here, this_filename = os.path.split(__file__) + +SPECTRA_CACHE_BASE = "/tmp/maria-data/atmosphere/spectra" +SPECTRA_SOURCE_BASE = ( + "https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/spectra" # noqa +) + + +class Spectrum: + def __init__(self, region, source="am", refresh_cache=False): + if region not in all_regions: + raise InvalidRegionError(region) + + self.region = region + self.source = source + + self.source_url = f"{SPECTRA_SOURCE_BASE}/{source}/{self.region}.h5" # noqa + self.cache_path = f"{SPECTRA_CACHE_BASE}/{source}/{self.region}.h5" # noqa + + fetch_cache( + source_url=self.source_url, + cache_path=self.cache_path, + max_age=30 * 86400, + refresh=refresh_cache, + ) + + with h5py.File(self.cache_path, "r") as f: + self._side_nu = f["side_nu_GHz"][:] + self._side_elevation = f["side_elevation_deg"][:] + self._side_zenith_pwv = f["side_zenith_pwv_mm"][:] + self._side_base_temperature = f["side_base_temperature_K"][:] + + self._emission = f["emission_temperature_rayleigh_jeans_K"][:] + self._transmission = np.exp(-f["opacity_nepers"][:]) + self._excess_path = 1e6 * ( + f["excess_path"][:] + f["offset_excess_path_m"][:] + ) + + def emission(self, nu, zenith_pwv=None, base_temperature=None, elevation=45): + if zenith_pwv is None: + zenith_pwv = np.median(self._side_zenith_pwv) + if base_temperature is None: + base_temperature = np.median(self._side_base_temperature) + + return sp.interpolate.RegularGridInterpolator( + points=( + self._side_zenith_pwv, + self._side_base_temperature, + self._side_elevation, + self._side_nu, + ), + values=self._emission, + )((zenith_pwv, base_temperature, elevation, nu)) + + def transmission(self, nu, zenith_pwv=None, base_temperature=None, elevation=45): + if zenith_pwv is None: + zenith_pwv = np.median(self._side_zenith_pwv) + if base_temperature is None: + base_temperature = np.median(self._side_base_temperature) + + return sp.interpolate.RegularGridInterpolator( + points=( + self._side_zenith_pwv, + self._side_base_temperature, + self._side_elevation, + self._side_nu, + ), + values=self._transmission, + )((zenith_pwv, base_temperature, elevation, nu)) diff --git a/maria/atmosphere/turbulent_layer.py b/maria/atmosphere/turbulent_layer.py index e1d77c1..102668a 100644 --- a/maria/atmosphere/turbulent_layer.py +++ b/maria/atmosphere/turbulent_layer.py @@ -1,9 +1,10 @@ import dask.array as da import numpy as np import scipy as sp +from todder.coords import Coordinates, get_center_phi_theta from .. import utils -from ..coords import Coordinates, get_center_phi_theta +from ..functions import approximate_normalized_matern from ..instrument import Instrument from ..instrument.beams import construct_beam_filter, separably_filter from .weather import Weather @@ -12,8 +13,7 @@ RIBBON_SAMPLE_DECAY = 2 JITTER_LEVEL = 1e-4 -matern_callback = utils.functions.approximate_normalized_matern -# matern_callback = utils.functions.normalized_matern +matern_callback = approximate_normalized_matern class TurbulentLayer: diff --git a/maria/atmosphere/weather.py b/maria/atmosphere/weather.py index 074126e..5b3804f 100644 --- a/maria/atmosphere/weather.py +++ b/maria/atmosphere/weather.py @@ -1,5 +1,4 @@ import os -from dataclasses import dataclass, field from datetime import datetime import h5py @@ -8,16 +7,14 @@ import scipy as sp from ..constants import g +from ..io import fetch_cache from ..site import InvalidRegionError, all_regions, supported_regions_table from ..utils import get_utc_day_hour, get_utc_year_day -from . import utils here, this_filename = os.path.split(__file__) -WEATHER_DATA_DIRECTORY = f"{here}/data" -WEATHER_DATA_CACHE_DIRECTORY = "/tmp/maria/weather" -WEATHER_DATA_URL_BASE = "https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/weather" # noqa F401 -CACHE_MAX_AGE_SECONDS = 30 * 86400 +WEATHER_CACHE_BASE = "/tmp/maria-data/weather" +WEATHER_SOURCE_BASE = "https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/weather" # noqa F401 def get_vapor_pressure(air_temp, rel_hum): # units are (°K, %) @@ -53,43 +50,43 @@ def absolute_to_relative_humidity(air_temp, abs_hum): return 1e2 * 461.5 * air_temp * abs_hum / get_saturation_pressure(air_temp) -@dataclass class Weather: - region: str = "chajnantor" - t: float = 0 - altitude: float = None - utc_time: str = "" - local_time: str = "" - time_zone: str = "" - quantiles: dict = field(default_factory=dict) - override: dict = field(default_factory=dict) - source: str = "era5" - from_cache: bool = None - - def __post_init__(self): - if self.region not in all_regions: + def __init__( + self, + region: str = "chajnantor", + t: float = 0, + altitude: float = None, + utc_time: str = "", + local_time: str = "", + time_zone: str = "", + quantiles: dict = {}, + override: dict = {}, + source: str = "era5", + refresh_cache: bool = False, + ): + if region not in all_regions: raise InvalidRegionError(self.region) - self.source_path = f"{WEATHER_DATA_DIRECTORY}/{self.source}/{self.region}.h5" - - # if the data isn't in the module, default to use the cache - self.from_cache = ( - self.from_cache - if self.from_cache is not None - else not os.path.exists(self.source_path) + self.region = region + self.t = t + self.altitude = altitude + self.utc_time = utc_time + self.local_time = local_time + self.time_zone = time_zone + self.quantiles = quantiles + self.override = override + self.source = source + + self.cache_path = f"{WEATHER_CACHE_BASE}/{source}/{region}.h5" + self.source_url = f"{WEATHER_SOURCE_BASE}/{source}/{region}.h5" + + fetch_cache( + source_url=self.source_url, + cache_path=self.cache_path, + max_age=30 * 86400, + refresh=refresh_cache, ) - # if we don't have the data, download and cache it - if self.from_cache: - self.source_path = ( - f"{WEATHER_DATA_CACHE_DIRECTORY}/{self.source}/{self.region}.h5" - ) - utils.io.fetch_cache( - source_url=f"{WEATHER_DATA_URL_BASE}/{self.source}/{self.region}.h5", - cache_path=self.source_path, - ) - self.from_cache = True - if self.altitude is None: self.altitude = supported_regions_table.loc[self.region, "altitude"] @@ -102,7 +99,7 @@ def __post_init__(self): pytz.timezone(self.time_zone) ).ctime() - with h5py.File(self.source_path, "r") as f: + with h5py.File(self.cache_path, "r") as f: self.utc_day_hour = get_utc_day_hour(self.t) self.utc_year_day = get_utc_year_day(self.t) diff --git a/maria/cmb/__init__.py b/maria/cmb/__init__.py index 38735f4..1c7f7f6 100644 --- a/maria/cmb/__init__.py +++ b/maria/cmb/__init__.py @@ -1,24 +1,25 @@ import healpy as hp import numpy as np import pandas as pd +import scipy as sp +from tqdm import tqdm -from maria.constants import T_CMB -from maria.utils.functions import planck_spectrum - -from ..utils.io import fetch_cache +from ..constants import T_CMB +from ..functions import planck_spectrum +from ..io import fetch_cache CMB_SPECTRUM_SOURCE_URL = ( "https://github.com/thomaswmorris/maria-data/raw/master/cmb/spectra/" "COM_PowerSpect_CMB-base-plikHM-TTTEEE-lowl-lowE-lensing-minimum-theory_R3.01.txt" ) -CMB_SPECTRUM_CACHE_PATH = "/tmp/maria/cmb/spectrum.txt" +CMB_SPECTRUM_CACHE_PATH = "/tmp/maria-data/cmb/spectrum.txt" CMB_SPECTRUM_CACHE_MAX_AGE = 30 * 86400 # one month CMB_MAP_SOURCE_URL = ( "https://irsa.ipac.caltech.edu/data/Planck/release_3/all-sky-maps/maps/component-maps/cmb/" "COM_CMB_IQU-143-fgsub-sevem_2048_R3.00_full.fits" ) -CMB_MAP_CACHE_PATH = "/tmp/maria/cmb/planck/map.fits" +CMB_MAP_CACHE_PATH = "/tmp/maria-data/cmb/planck/map.fits" CMB_MAP_CACHE_MAX_AGE = 30 * 86400 # one month @@ -28,7 +29,7 @@ def __init__(self, data, fields): raise ValueError("Data and labels must have the same shape!") self.maps = {} - for i, (field, M) in enumerate(zip(fields, data)): + for field, M in zip(fields, data): self.maps[field] = M self.nside = int(np.sqrt(len(M) / 12)) @@ -51,7 +52,7 @@ def plot(self, field=None, units="uK"): ) -def generate_cmb(nside=2048, seed=123): +def generate_cmb(nside=1024, seed=123456, **kwargs): """ Taken from https://www.zonca.dev/posts/2020-09-30-planck-spectra-healpy.html """ @@ -61,7 +62,7 @@ def generate_cmb(nside=2048, seed=123): fetch_cache( source_url=CMB_SPECTRUM_SOURCE_URL, cache_path=CMB_SPECTRUM_CACHE_PATH, - CACHE_MAX_AGE=CMB_SPECTRUM_CACHE_MAX_AGE, + max_age=CMB_SPECTRUM_CACHE_MAX_AGE, ) cl = pd.read_csv(CMB_SPECTRUM_CACHE_PATH, delim_whitespace=True, index_col=0) @@ -81,11 +82,11 @@ def generate_cmb(nside=2048, seed=123): return cmb -def get_cmb(): +def get_cmb(**kwargs): fetch_cache( source_url=CMB_MAP_SOURCE_URL, cache_path=CMB_MAP_CACHE_PATH, - CACHE_MAX_AGE=CMB_MAP_CACHE_MAX_AGE, + max_age=CMB_MAP_CACHE_MAX_AGE, ) field_dtypes = { @@ -109,41 +110,40 @@ def get_cmb(): class CMBMixin: - def _initialize_cmb(self, source): - if source == "generate": - if self.verbose: - print("Generating CMB realization...") - self.cmb = generate_cmb() - - elif source == "map": - if self.verbose: - print("Loading CMB...") - self.cmb = get_cmb() - def _simulate_cmb_emission(self): pixel_index = hp.ang2pix( nside=self.cmb.nside, phi=self.coords.l, theta=np.pi / 2 - self.coords.b - ).compute() + ).compute() # noqa cmb_temperatures = self.cmb.T[pixel_index] - test_nu = np.linspace(1e9, 1e12, 1024) + test_nu = np.linspace(1e0, 1e3, 1024) cmb_temperature_samples_K = T_CMB + np.linspace( - self.cmb.T.min(), self.cmb.T.max(), 64 + self.cmb.T.min(), self.cmb.T.max(), 3 + ) # noqa + cmb_brightness = planck_spectrum( + 1e9 * test_nu, cmb_temperature_samples_K[:, None] ) - cmb_brightness = planck_spectrum(test_nu, cmb_temperature_samples_K[:, None]) self.data["cmb"] = np.zeros((self.instrument.dets.n, self.plan.n_time)) - for band in self.instrument.bands: + pbar = tqdm(self.instrument.bands, disable=not self.verbose) + + for band in pbar: + pbar.set_description(f"Sampling CMB ({band.name})") + band_mask = self.instrument.dets.band_name == band.name - band_cmb_power_samples_W = np.trapz( - y=cmb_brightness * band.passband(1e-9 * test_nu), x=test_nu + band_cmb_power_samples_W = ( + 1e12 + * band.efficiency + * np.trapz(y=cmb_brightness * band.passband(test_nu), x=test_nu) ) - self.data["cmb"][band_mask] = np.interp( - T_CMB + cmb_temperatures[band_mask], + # dP_dTCMB = self.instrument.dets.dP_dTCMB[:, None] + # self.data["cmb"][band_mask] = * cmb_temperatures[band_mask] + + self.data["cmb"][band_mask] = sp.interpolate.interp1d( cmb_temperature_samples_K, band_cmb_power_samples_W, - ) + )(T_CMB + cmb_temperatures[band_mask]) diff --git a/maria/coords/__init__.py b/maria/coords/__init__.py deleted file mode 100644 index a33d416..0000000 --- a/maria/coords/__init__.py +++ /dev/null @@ -1,463 +0,0 @@ -import functools -import time as ttime -from typing import Any - -import dask.array as da -import numpy as np -import pandas as pd -import scipy as sp -from astropy import units as u -from astropy.coordinates import EarthLocation, SkyCoord -from astropy.time import Time -from scipy.interpolate import interp1d - - -def now(): - return ttime.time() - - -def dx_dy_to_phi_theta(dx, dy, cphi, ctheta): - """ - A fast and well-conditioned to convert from local dx/dy coordinates to phi/theta coordinates. - """ - - if not dx.shape == dy.shape: - raise ValueError( - f"The shapes of 'dx' and 'dy' must be the same. Got shapes {np.shape(dx)} and {np.shape(dy)}" - ) - - r = np.sqrt(dx**2 + dy**2) # distance from the center - p = np.arctan2(dx, -dy) # 0 at the bottom, increases CCW to pi at the top - - # if we're looking at the north pole, we have (lon, lat) = (p, pi/2 - r) - # a projection looking from the east - proj_from_east = (np.sin(r) * np.cos(p) + 1j * np.cos(r)) * np.exp( - 1j * (ctheta - np.pi / 2) - ) - phi = cphi + np.arctan2(np.sin(r) * np.sin(p), np.real(proj_from_east)) - theta = np.arcsin(np.imag(proj_from_east)) - - return ( - phi, - theta, - ) - - -def phi_theta_to_dx_dy(phi, theta, cphi, ctheta): - """ - A fast and well-conditioned to convert from phi/theta coordinates to local dx/dy coordinates. - """ - - if not phi.shape == theta.shape: - raise ValueError( - f"The shapes of 'phi' and 'theta' must be the same. Got shapes {np.shape(phi)} and {np.shape(theta)}" - ) - - dphi = phi - cphi - proj_from_east = (np.cos(dphi) * np.cos(theta) + 1j * np.sin(theta)) * np.exp( - 1j * (np.pi / 2 - ctheta) - ) - dz = np.sin(dphi) * np.cos(theta) + 1j * np.real(proj_from_east) - r = np.abs(dz) - dz *= np.arcsin(r) / r - - # negative, because we're looking at the observer - return np.real(dz), -np.imag(dz) - - -def phi_theta_to_xyz(phi, theta): - """ - Project a longitude and lattitude onto the unit sphere. - """ - # you can add a newaxis on numpy floats, but not python floats. who knew? - return np.concatenate( - [ - (np.cos(phi) * np.cos(theta))[..., None], - (np.sin(phi) * np.cos(theta))[..., None], - (np.sin(theta))[..., None], - ], - axis=-1, - ) - - -def xyz_to_phi_theta(xyz): - """ - Find the longitude and latitude of a 3-vector. - """ - return np.arctan2(xyz[..., 1], xyz[..., 0]) % (2 * np.pi), np.arcsin( - xyz[..., 2] / np.sqrt(np.sum(xyz**2, axis=-1)) - ) - - -def get_center_phi_theta(phi, theta, keep_last_dim=False): - """ """ - xyz = phi_theta_to_xyz(phi, theta) - - if keep_last_dim: - center_xyz = xyz.mean(axis=tuple(range(xyz.ndim - 2))) - center_xyz /= np.sqrt(np.sum(np.square(center_xyz), axis=-1))[..., None] - else: - center_xyz = xyz.mean(axis=tuple(range(xyz.ndim - 1))) - center_xyz /= np.sqrt(np.sum(np.square(center_xyz))) - - return xyz_to_phi_theta(center_xyz) - - -frames = { - "az_el": { - "astropy_name": "altaz", - "astropy_phi": "az", - "astropy_theta": "alt", - "phi": "az", - "theta": "el", - "phi_name": "RA", - "theta_name": "Dec.", - }, - "ra_dec": { - "astropy_name": "icrs", - "astropy_phi": "ra", - "astropy_theta": "dec", - "phi": "ra", - "theta": "dec", - "phi_name": "RA", - "theta_name": "Dec.", - }, - "galactic": { - "astropy_name": "galactic", - "astropy_phi": "l", - "astropy_theta": "b", - "phi": "l", - "theta": "b", - "phi_name": "RA", - "theta_name": "Dec.", - }, -} - - -class Angle: - def __init__(self, a, units="radians"): - if units == "radians": - self.a = a - elif units == "degrees": - self.a = (np.pi / 180) * a - elif units == "arcmin": - self.a = (np.pi / 180 / 60) * a - elif units == "arcsec": - self.a = (np.pi / 180 / 3600) * a - else: - raise ValueError( - "'units' must be one of ['radians', 'degrees', 'arcmin', 'arcsec']" - ) - - self.is_scalar = len(np.shape(self.a)) == 0 - - if not self.is_scalar: - self.a = np.unwrap(self.a) - - @property - def radians(self): - return self.a - - @property - def rad(self): - return self.radians - - @property - def degrees(self): - return (180 / np.pi) * self.a - - @property - def deg(self): - return self.degrees - - @property - def arcmin(self): - return (60 * 180 / np.pi) * self.a - - @property - def arcsec(self): - return (3600 * 180 / np.pi) * self.a - - def __repr__(self): - return f"Angle(units='radians', value={self.a.__repr__()})" - - @property - def units(self): - # peak-to-peak - max_deg = self.deg if self.is_scalar else self.deg.max() - - if max_deg < 0.5 / 60: - return "arcsec" - if max_deg < 0.5: - return "arcmin" - - return "degrees" - - -class Coordinates: - """ - A class for managing coordinates, allowing us to access different coordinate frames at once. - """ - - def __init__( - self, - time: float, - phi: float, - theta: float, - location: EarthLocation, - frame: str = "ra_dec", - dtype=np.float32, - time_offset: float = 0, - ): - if float(time.min()) != 0: - time_offset = time.min() - time -= time_offset - - self.phi = np.atleast_1d( - (phi if isinstance(phi, da.Array) else da.from_array(phi)) - ).astype(dtype) - self.theta = np.atleast_1d( - (theta if isinstance(theta, da.Array) else da.from_array(theta)) - ).astype(dtype) - self.time = np.atleast_1d( - (time if isinstance(time, da.Array) else da.from_array(time)) - ).astype(dtype) - - self.timestep = np.median(np.gradient(self.time)) - - if self.time.ndim > 1: - raise ValueError("'time' can be at most one-dimensional.") - if ( - len(self.time) != self.phi.shape[-1] - or len(self.time) != self.theta.shape[-1] - ): - raise ValueError( - "The size of the last dimensions of [phi, theta, time] must all match." - ) - - *self.input_shape, self.time_shape = self.phi.shape - - self.location = location - self.frame = frame - self.time_offset = time_offset - self.dtype = dtype - - _center_phi, _center_theta = get_center_phi_theta( - self.phi, self.theta, keep_last_dim=True - ) - - # three fiducial offsets from the boresight to train a transformation matrix - fid_offsets = np.radians([[-1.0, 0.0], [1.0, 0.0], [0.0, 1.0]]) - - # the time-ordered coordinates of those offsets in the frame - time_ordered_fid_phi, time_ordered_fid_theta = dx_dy_to_phi_theta( - *fid_offsets.T, _center_phi[:, None], _center_theta[:, None] - ) - - time_ordered_fid_points = phi_theta_to_xyz( - time_ordered_fid_phi, time_ordered_fid_theta - ) - - # time_offset = _time.min().compute() - # end_time = _time.max().compute() - duration = self.time.max().compute() - - self.sampled_time = np.linspace(0, duration, int(np.maximum(2, duration / 1))) - - self.sampled_fid_points = interp1d(self.time, time_ordered_fid_points, axis=0)( - self.sampled_time - ) - self.sampled_fid_points_inv = np.linalg.inv(self.sampled_fid_points) - sampled_fid_phi, sampled_fid_theta = xyz_to_phi_theta(self.sampled_fid_points) - - # call astropy to compute a literal representation of the sampled fiducial points - self.sampled_fid_skycoords = SkyCoord( - sampled_fid_phi * u.rad, - sampled_fid_theta * u.rad, - obstime=Time( - self.time_offset + np.kron(np.ones((3, 1)), self.sampled_time).T, - format="unix", - ), - frame=frames[frame]["astropy_name"], - location=location, - ) - - self.TRANSFORMS = {} - - def downsample(self, timestep: float = None, factor: int = None): - if timestep is None and factor is None: - raise ValueError("You must supply either 'timestep' or 'factor'.") - - timestep = timestep or factor * self.timestep - - ds_time = np.arange(self.time.min(), self.time.max(), timestep) - ds_phi = sp.interpolate.interp1d(self.time, self.phi, axis=-1)(ds_time) - ds_theta = sp.interpolate.interp1d(self.time, self.theta, axis=-1)(ds_time) - - return Coordinates( - time=ds_time, - phi=ds_phi, - theta=ds_theta, - location=self.location, - frame=self.frame, - dtype=self.dtype, - time_offset=self.time_offset, - ) - - @functools.cached_property - def boresight(self): - cphi, ctheta = get_center_phi_theta(self.phi, self.theta, keep_last_dim=True) - - return Coordinates( - time=self.time, - phi=cphi, - theta=ctheta, - location=self.location, - frame=self.frame, - dtype=self.dtype, - time_offset=self.time_offset, - ) - - @functools.cached_property - def summary(self): - # compute summary for the string repr - summary = pd.DataFrame(columns=["min", "mean", "max"]) - boresight = self.boresight - for attr in ["az", "el", "ra", "dec"]: - for stat in ["min", "mean", "max"]: - summary.loc[ - attr, stat - ] = f"{float(np.degrees(getattr(getattr(boresight, attr), stat)().compute())):.03f}°" - - return summary - - def compute_points(self): - return phi_theta_to_xyz(self.phi, self.theta) - - def to_frame(self, frame): - """ - Convert to a frame. This is expensive to compute, so we cache the output. - """ - - if frame == self.frame: - return self.phi, self.theta - - elif frame == "az_el": - sampled_fid_skycoords_new_frame = self.sampled_fid_skycoords.altaz - - elif frame == "ra_dec": - sampled_fid_skycoords_new_frame = self.sampled_fid_skycoords.icrs - - elif frame == "galactic": - sampled_fid_skycoords_new_frame = self.sampled_fid_skycoords.galactic - - else: - raise ValueError("'frame' must be one of ['az_el', 'ra_dec', 'galactic']") - - # find a rotation matrix between the init frame and the new frame - sampled_fid_phi_new_frame = getattr( - sampled_fid_skycoords_new_frame, frames[frame]["astropy_phi"] - ).rad - sampled_fid_theta_new_frame = getattr( - sampled_fid_skycoords_new_frame, frames[frame]["astropy_theta"] - ).rad - sampled_fid_points_new_frame = phi_theta_to_xyz( - sampled_fid_phi_new_frame, sampled_fid_theta_new_frame - ) - sampled_rotation_matrix = ( - self.sampled_fid_points_inv @ sampled_fid_points_new_frame - ).swapaxes(-2, -1) - - # apply the matrix and convert to coordinates - extra_dims = tuple(range(len(self.input_shape))) - transform = np.expand_dims( - interp1d(self.sampled_time, sampled_rotation_matrix, axis=0)(self.time), - axis=extra_dims, - ) - return xyz_to_phi_theta( - ( - transform.astype(self.dtype) - @ np.expand_dims(self.compute_points(), axis=-1) - ).squeeze() - ) - - @functools.cached_property - def az_el(self): - return self.to_frame(frame="az_el") - - @functools.cached_property - def center_az_el(self): - return get_center_phi_theta(*self.az_el) - - @functools.cached_property - def ra_dec(self): - return self.to_frame(frame="ra_dec") - - @functools.cached_property - def center_ra_dec(self): - return get_center_phi_theta(*self.ra_dec) - - @functools.cached_property - def galactic(self): - return self.to_frame(frame="galactic") - - @functools.cached_property - def center_galactic(self): - return get_center_phi_theta(*self.galactic) - - def __getattr__(self, attr: str) -> Any: - if attr == "az": - return self.az_el[0] - if attr == "el": - return self.az_el[1] - if attr == "center_az": - return self.dtype(self.center_az_el[0]) - if attr == "center_el": - return self.dtype(self.center_az_el[1]) - - if attr == "ra": - return self.ra_dec[0] - if attr == "dec": - return self.ra_dec[1] - if attr == "center_ra": - return self.dtype(self.center_ra_dec[0]) - if attr == "center_dec": - return self.dtype(self.center_ra_dec[1]) - - if attr == "l": - return self.galactic[0] - if attr == "b": - return self.galactic[1] - if attr == "center_l": - return self.dtype(self.center_galactic[0]) - if attr == "center_b": - return self.dtype(self.center_galactic[1]) - - raise AttributeError(f"Coordinates object has no attribute named '{attr}'.") - - def offsets(self, frame, center="auto", units="radians"): - if isinstance(center, str): - if center == "auto": - center = get_center_phi_theta(*self.to_frame(frame)) - - if frame == "az_el": - dx, dy = phi_theta_to_dx_dy(self.az, self.el, *center) - elif frame == "ra_dec": - dx, dy = phi_theta_to_dx_dy(self.ra, self.dec, *center) - elif frame == "galactic": - dx, dy = phi_theta_to_dx_dy(self.l, self.b, *center) - - if units == "radians": - return dx, dy - elif units == "degrees": - return np.degrees(dx), np.degrees(dy) - elif units == "arcmin": - return 60 * np.degrees(dx), 60 * np.degrees(dy) - elif units == "arcsec": - return 3600 * np.degrees(dx), 3600 * np.degrees(dy) - - def __repr__(self): - return self.summary.__repr__() - - def _repr_html_(self): - return self.summary._repr_html_() diff --git a/maria/errors.py b/maria/errors.py index b933bf4..370019a 100644 --- a/maria/errors.py +++ b/maria/errors.py @@ -1,2 +1,6 @@ class PointingError(BaseException): ... + + +class ConfigurationError(Exception): + ... diff --git a/maria/utils/functions.py b/maria/functions.py similarity index 91% rename from maria/utils/functions.py rename to maria/functions.py index 756254b..5dfab2d 100644 --- a/maria/utils/functions.py +++ b/maria/functions.py @@ -4,7 +4,11 @@ from maria.constants import c, h, k_B -def planck_spectrum(nu, T): +def rayleigh_jeans_spectrum(nu: float, T: float): + return 2 * k_B * nu**2 * T / c**2 + + +def planck_spectrum(nu: float, T: float): return 2 * h * nu**3 / (c**2 * np.expm1(h * nu / (k_B * T))) diff --git a/maria/instrument/__init__.py b/maria/instrument/__init__.py index a301c42..8a93355 100644 --- a/maria/instrument/__init__.py +++ b/maria/instrument/__init__.py @@ -10,13 +10,12 @@ import pandas as pd from matplotlib.collections import EllipseCollection from matplotlib.patches import Patch +from todder.coords import Angle -from .. import utils -from ..coords import Angle -from .arrays import generate_array -from .bands import BandList, all_bands, parse_bands_config # noqa F401 -from .beams import compute_angular_fwhm # noqa F401 -from .detectors import Detectors # noqa F401 +from ..io import read_yaml +from .bands import BandList, all_bands, parse_bands # noqa +from .beams import compute_angular_fwhm +from .detectors import Detectors, generate_array HEX_CODE_LIST = [ mpl.colors.to_hex(mpl.colormaps.get_cmap("Paired")(t)) @@ -28,9 +27,7 @@ here, this_filename = os.path.split(__file__) -# all_instrument_params = utils.io.read_yaml(f"{here}/../configs/default_params.yml")["instrument"] - -INSTRUMENT_CONFIGS = utils.io.read_yaml(f"{here}/instruments.yml") +INSTRUMENT_CONFIGS = read_yaml(f"{here}/configs.yml") INSTRUMENT_DISPLAY_COLUMNS = [ "description", @@ -80,12 +77,12 @@ def get_instrument(instrument_name="default", **kwargs): "field_of_view", ] -band_params_to_inherit = { - "time_constant": 0.0, # seconds - "white_noise": 0.0, # Kelvin - "pink_noise": 0.0, # Kelvin / s - "efficiency": 1.0, -} +band_params_to_inherit = [ + "time_constant", # seconds + "white_noise", # Kelvin + "pink_noise", # Kelvin / s + "efficiency", +] passband_params_to_inherit = { "center": 150, # GHz @@ -137,12 +134,16 @@ def get_subarrays(instrument_config): for subarray_name in config["arrays"]: subarray = config["arrays"][subarray_name] - subarray["bands"] = parse_bands_config(subarray.get("bands")) + subarray["bands"] = parse_bands(subarray.get("bands")) if "file" in subarray: # it points to a file: - subarray["file"] = f"{here}/{subarray['file']}" + if not os.path.exists(subarray["file"]): + subarray["file"] = f"{here}/detectors/arrays/{subarray['file']}" df = pd.read_csv(subarray["file"]) + if subarray["bands"] is None: + subarray["bands"] = {} + subarray["n"] = len(df) for band_name in np.unique(df.band_name.values): subarray["bands"][band_name] = all_bands[band_name] @@ -157,13 +158,14 @@ def get_subarrays(instrument_config): else: raise ValueError("You must pass 'bands' for each array.") - for band_name, band_config in subarray["bands"].items(): - for param, default_value in band_params_to_inherit.items(): - band_config[param] = band_config.get(param, default_value) + # for band_name, band_config in subarray["bands"].items(): + # for param in band_params_to_inherit: + # if param in config: + # band_config[param] = band_config.get(param, default_value) - if "passband" not in band_config: - for param, default_value in passband_params_to_inherit.items(): - band_config[param] = band_config.get(param, default_value) + # if "passband" not in band_config: + # for param, default_value in passband_params_to_inherit: + # band_config[param] = band_config.get(param, default_value) subarrays[subarray_name] = subarray diff --git a/maria/instrument/bands.py b/maria/instrument/bands.py deleted file mode 100644 index 5ef7e18..0000000 --- a/maria/instrument/bands.py +++ /dev/null @@ -1,226 +0,0 @@ -import glob -import os -from collections.abc import Mapping -from dataclasses import dataclass -from typing import Sequence - -import numpy as np -import pandas as pd - -from .. import utils -from ..constants import c, k_B -from ..utils.io import flatten_config, read_yaml - -BAND_FIELD_TYPES = { - "center": "float", - "width": "float", - "shape": "str", - "time_constant": "float", - "white_noise": "float", - "pink_noise": "float", - "efficiency": "float", -} - -here, this_filename = os.path.split(__file__) - -all_bands = {} -for path in glob.glob(f"{here}/data/bands/*.yml"): - tag = os.path.split(path)[1].split(".")[0] - all_bands[tag] = read_yaml(path) -all_bands = {} -for path in glob.glob(f"{here}/data/bands/*.yml"): - tag = os.path.split(path)[1].split(".")[0] - all_bands[tag] = read_yaml(path) - -all_bands = flatten_config(all_bands) - - -def validate_band_config(band): - if "passband" not in band: - if any([key not in band for key in ["center", "width"]]): - raise ValueError("The band's center and width must be specified!") - - -def parse_bands_config(bands): - """ - There are many ways to specify bands, and this handles them. - """ - parsed_band_config = {} - - if isinstance(bands, Mapping): - for name, band in bands.items(): - validate_band_config(band) - parsed_band_config[name] = band - return parsed_band_config - - if isinstance(bands, list): - for band in bands: - if isinstance(band, str): - if band not in all_bands: - raise ValueError(f'Band "{band}" is not supported.') - parsed_band_config[band] = all_bands[band] - - if isinstance(band, Mapping): - validate_band_config(band) - name = band.get("name", f'f{int(band["center"]):>03}') - parsed_band_config[name] = band - - return parsed_band_config - - -class BandList(Sequence): - @classmethod - def from_config(cls, config): - bands = [] - - if isinstance(config, str): - config = utils.io.read_yaml(f"{here}/{config}") - - for name in config.keys(): - band_config = config[name] - if "file" in band_config.keys(): - band_config = utils.io.read_yaml(f'{here}/{band_config["file"]}') - - bands.append(Band(name=name, **band_config)) - return cls(bands=bands) - - def add(self, band): - if band.name not in self.names: - self.bands.append(band) - else: - raise ValueError() - - def __init__(self, bands: list = []): - self.bands = bands - - def __getattr__(self, attr): - if attr in self.names: - return self.__getitem__(attr) - if [hasattr(band, attr) for band in self.bands]: - return [getattr(band, attr) for band in self.bands] - raise AttributeError(f"BandList object has no attribute named '{attr}'.") - - def __getitem__(self, index): - if type(index) is int: - return self.bands[index] - elif type(index) is str: - if index not in self.names: - raise ValueError(f"BandList has no band named {index}.") - return self.bands[self.names.index(index)] - else: - raise ValueError( - f"Invalid index {index}. A bandList must be indexed by either an integer or a string." - ) - - def __len__(self): - return len(self.bands) - - def __repr__(self): - return self.summary.T.__repr__() - - def __repr_html__(self): - return self.summary.T.__repr_html__() - - def __short_repr__(self): - return f"BandList([{', '.join(self.names)}])" - - @property - def names(self): - return [band.name for band in self.bands] - - @property - def summary(self) -> pd.DataFrame: - table = pd.DataFrame(columns=list(BAND_FIELD_TYPES.keys()), index=self.names) - - for attr, dtype in BAND_FIELD_TYPES.items(): - for band in self.bands: - table.at[band.name, attr] = getattr(band, attr) - table[attr] = table[attr].astype(dtype) - - return table - - -@dataclass -class Band: - name: str - center: float - width: float - time_constant: float = 0.0 - white_noise: float = 0.0 - pink_noise: float = 0.0 - shape: str = "top_hat" - efficiency: float = 1.0 - - @classmethod - def from_config(cls, name, config): - if "passband" in config: - df = pd.read_csv(f"{here}/{config.pop('passband')}", index_col=0) - - nu, pb = df.nu.values, df.passband.values - - center = np.round(np.sum(pb * nu), 3) - width = np.round(nu[pb > 1e-2 * pb.max()].ptp(), 3) - - band = cls(name=name, center=center, width=width, shape="custom", **config) - band._nu = nu - band._pb = pb - return band - - else: - return cls(name=name, **config) - - @property - def nu_min(self) -> float: - if self.shape == "flat": - return self.center - 0.5 * 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 == "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). - """ - _nu = np.atleast_1d(nu) - - if self.shape == "top_hat": - return np.exp(-np.log(2) * (2 * (_nu - self.center) / self.width) ** 8) - - 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 pW_to_KRJ(self): - """ - Absolute calibration (i.e. temperature per power) in Kelvins per picowatt. - """ - nu = np.linspace(self.nu_min, self.nu_max, 256) - dP_dT = ( - self.efficiency - * np.trapz(2 * k_B * (1e9 * nu) ** 2 * self.passband(nu), 1e9 * nu) - / c**2 - ) - return 1e-12 / dP_dT - - @property - def wavelength(self): - """ - Return the wavelength of the center, in meters. - """ - return c / (1e9 * self.center) diff --git a/maria/instrument/bands/__init__.py b/maria/instrument/bands/__init__.py new file mode 100644 index 0000000..0a6bd16 --- /dev/null +++ b/maria/instrument/bands/__init__.py @@ -0,0 +1,328 @@ +import copy +import glob +import os +from collections.abc import Mapping +from typing import Sequence + +import numpy as np +import pandas as pd + +from ...atmosphere.spectrum import Spectrum +from ...constants import T_CMB, c +from ...functions import planck_spectrum, rayleigh_jeans_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__) + +all_bands = {} +for path in glob.glob(f"{here}/configs/*.yml"): + tag = os.path.split(path)[1].split(".")[0] + all_bands[tag] = read_yaml(path) +all_bands = {} +for path in glob.glob(f"{here}/configs/*.yml"): + tag = os.path.split(path)[1].split(".")[0] + all_bands[tag] = read_yaml(path) + +all_bands = flatten_config(all_bands) + + +def validate_band_config(band): + if "passband" not in band: + if any([key not in band for key in ["center", "width"]]): + raise ValueError("The band's center and width must be specified!") + + +def parse_bands(bands): + """ + There are many ways to specify bands, and this handles all of them. + """ + bands = copy.deepcopy(bands) + + if isinstance(bands, Mapping): + for name, band in bands.items(): + validate_band_config(band) + + # if we get a list of bands, convert them to a labeled dict + if isinstance(bands, list): + bands_mapping = {} + for band in bands: + if isinstance(band, str): + # here 'band' is a name + if band not in all_bands: + raise ValueError(f"Could not find band '{band}'.") + bands_mapping[band] = all_bands[band] + + if isinstance(band, Mapping): + if "center" not in band: + raise RuntimeError("You must specify the band center.") + name = band.get("name", f'f{int(band["center"]):>03}') + bands_mapping[name] = band + + return parse_bands(bands_mapping) + + return bands + + +class Band: + def __init__( + self, + name: str, + center: float, + width: float, + shape: str = "gaussian", + efficiency: float = 1.0, + sensitivity: float = None, + sensitivity_kind: str = "rayleigh-jeans", + NEP: float = None, + NEP_per_loading: float = 0.0, + knee: float = 1.0, + time_constant: float = 0.0, + ): + self.name = name + self.center = center + self.width = width + self.shape = shape + self.efficiency = efficiency + self.NEP_per_loading = NEP_per_loading + + self.knee = knee + self.time_constant = time_constant + + self.spectrum = Spectrum(region="chajnantor") + + if NEP is not None: + if sensitivity is not None: + raise RuntimeError( + "When defining a band, you must specify exactly one of 'NEP' or 'sensitivity'." + ) + self.NEP = NEP + + if sensitivity is not None: + self.set_sensitivity( + sensitivity, kind=sensitivity_kind + ) # this sets the NEP automatically + + def __repr__(self): + 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}" + parts.append(s) + + return f"Band({', '.join(parts)})" + + @property + def sensitivity(self): + if not hasattr(self, "_sensitivity"): + self._sensitivity = self.NEP * self.transmission() / self.dP_dTRJ + return self._sensitivity + + @sensitivity.setter + def sensitivity(self, value): + self.set_sensitivity(value) + + def set_sensitivity( + self, value, kind="rayleigh-jeans", region="chajnantor", pwv=0, elevation=90 + ): + transmission = self.transmission( + region=region, zenith_pwv=pwv, elevation=elevation + ) + + if kind.lower() == "rayleigh-jeans": + self.NEP = self.dP_dTRJ * value / transmission + + elif kind.lower() == "cmb": + self.NEP = self.dP_dTCMB * value / transmission + + self._sensitivity = value + self._sensitivity_kind = kind + + def transmission(self, region="chajnantor", zenith_pwv=1, elevation=90) -> float: + if self.spectrum.region != region: + self.spectrum = Spectrum(region=region) + return self.spectrum.transmission( + nu=self.center, zenith_pwv=zenith_pwv, elevation=elevation + ) + + @classmethod + def from_config(cls, name, config): + if "passband" in config: + df = pd.read_csv(f"{here}/{config.pop('passband')}", index_col=0) + + nu, pb = df.nu.values, df.passband.values + + center = np.round(np.sum(pb * nu), 3) + width = np.round(nu[pb > 1e-2 * pb.max()].ptp(), 3) + + band = cls(name=name, center=center, width=width, shape="custom", **config) + band._nu = nu + band._pb = pb + else: + band = cls(name=name, **config) + + return band + + @property + def nu_min(self) -> float: + if self.shape == "flat": + return self.center - 0.5 * 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 == "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). + """ + _nu = np.atleast_1d(nu) + + if self.shape == "top_hat": + return np.exp(np.log(0.5) * (2 * (_nu - self.center) / self.width) ** 8) + + if self.shape == "gaussian": + return np.exp(np.log(0.5) * (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 dP_dTRJ(self) -> float: + """ + In watts 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) + + return self.efficiency * dP_dTRJ + + @property + def dP_dTCMB(self) -> float: + """ + In watts per kelvin CMB, assuming perfect transmission. + """ + + eps = 1e-3 + + nu = np.linspace(self.nu_min, self.nu_max, 256) + + 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] + / eps + ) + dP_dTCMB = self.efficiency * np.trapz(dI_dTCMB * self.passband(nu), 1e9 * nu) + + return self.efficiency * dP_dTCMB + + @property + def wavelength(self): + """ + Return the wavelength of the center, in meters. + """ + return c / (1e9 * self.center) + + +class BandList(Sequence): + @classmethod + def from_config(cls, config): + bands = [] + + if isinstance(config, str): + config = read_yaml(f"{here}/{config}") + + for name in config.keys(): + band_config = config[name] + if "file" in band_config.keys(): + band_config = read_yaml(f'{here}/{band_config["file"]}') + + bands.append(Band(name=name, **band_config)) + return cls(bands=bands) + + def add(self, band): + if not isinstance(band, Band): + raise ValueError("'band' must be a Band type.") + if band.name in self.names: + raise RuntimeError(f"There is already a band called '{band.name}'.") + self.bands.append(band) + + def __init__(self, bands: list = []): + self.bands = bands + + def __getattr__(self, attr): + if attr in self.names: + return self.__getitem__(attr) + if all([hasattr(band, attr) for band in self.bands]): + return [getattr(band, attr) for band in self.bands] + raise AttributeError(f"BandList object has no attribute named '{attr}'.") + + def __getitem__(self, index): + if type(index) is int: + return self.bands[index] + elif type(index) is str: + if index not in self.names: + raise ValueError(f"BandList has no band named {index}.") + return self.bands[self.names.index(index)] + else: + raise ValueError( + f"Invalid index {index}. A bandList must be indexed by either an integer or a string." + ) + + def __len__(self): + return len(self.bands) + + def __repr__(self): + return self.summary.__repr__() + + def _repr_html_(self): + return self.summary._repr_html_() + + def __short_repr__(self): + return f"BandList([{', '.join(self.names)}])" + + @property + def names(self): + return [band.name for band in self.bands] + + @property + def summary(self) -> pd.DataFrame: + table = pd.DataFrame(columns=list(BAND_FIELDS.keys()), 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) + + return table diff --git a/maria/instrument/bands/configs/abs.yml b/maria/instrument/bands/configs/abs.yml new file mode 100644 index 0000000..f1b7361 --- /dev/null +++ b/maria/instrument/bands/configs/abs.yml @@ -0,0 +1,7 @@ +f150: + center: 150 + width: 30 + time_constant: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz + efficiency: 0.5 diff --git a/maria/instrument/data/bands/act.yml b/maria/instrument/bands/configs/act.yml similarity index 54% rename from maria/instrument/data/bands/act.yml rename to maria/instrument/bands/configs/act.yml index 8bec284..e99a8fc 100644 --- a/maria/instrument/data/bands/act.yml +++ b/maria/instrument/bands/configs/act.yml @@ -2,50 +2,50 @@ pa4: f150: center: 150 width: 50 - shape: top_hat + shape: gaussian time_constant: 0 - white_noise: 266.e-6 - pink_noise: 1.e-1 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 0.5 f220: center: 220 width: 60 - shape: top_hat + shape: gaussian time_constant: 0 - white_noise: 266.e-6 - pink_noise: 1.e-1 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 0.5 pa5: f090: center: 90 width: 30 - shape: top_hat + shape: gaussian time_constant: 0 - white_noise: 266.e-6 - pink_noise: 1.e-1 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 0.5 f150: center: 150 width: 50 - shape: top_hat + shape: gaussian time_constant: 0 - white_noise: 266.e-6 - pink_noise: 1.e-1 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 0.5 pa6: f090: center: 90 width: 30 - shape: top_hat + shape: gaussian time_constant: 0 - white_noise: 266.e-6 - pink_noise: 1.e-1 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 0.5 f150: center: 150 width: 50 - shape: top_hat + shape: gaussian time_constant: 0 - white_noise: 266.e-6 - pink_noise: 1.e-1 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 0.5 diff --git a/maria/instrument/data/bands/alma.yml b/maria/instrument/bands/configs/alma.yml similarity index 53% rename from maria/instrument/data/bands/alma.yml rename to maria/instrument/bands/configs/alma.yml index 3b0a12d..e226057 100644 --- a/maria/instrument/data/bands/alma.yml +++ b/maria/instrument/bands/configs/alma.yml @@ -1,89 +1,89 @@ f043: center: 43 width: 16 - shape: top_hat + shape: gaussian time_constant: 0 - white_noise: 0 - pink_noise: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 1.0 f078: center: 78 width: 22 - shape: top_hat + shape: gaussian time_constant: 0 - white_noise: 0 - pink_noise: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 1.0 f100: center: 100 width: 32 - shape: top_hat + shape: gaussian time_constant: 0 - white_noise: 0 - pink_noise: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 1.0 f144: center: 144 width: 38 - shape: top_hat + shape: gaussian time_constant: 0 - white_noise: 0 - pink_noise: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 1.0 f187: center: 187 width: 48 - shape: top_hat + shape: gaussian time_constant: 0 - white_noise: 0 - pink_noise: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 1.0 f243: center: 243 width: 64 - shape: top_hat + shape: gaussian time_constant: 0 - white_noise: 0 - pink_noise: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 1.0 f324: center: 324 width: 98 - shape: top_hat + shape: gaussian time_constant: 0 - white_noise: 0 - pink_noise: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 1.0 f447: center: 447 width: 114 - shape: top_hat + shape: gaussian time_constant: 0 - white_noise: 0 - pink_noise: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 1.0 f661: center: 661 width: 118 - shape: top_hat + shape: gaussian time_constant: 0 - white_noise: 0 - pink_noise: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 1.0 f869: center: 869 width: 163 - shape: top_hat + shape: gaussian time_constant: 0 - white_noise: 0 - pink_noise: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 1.0 diff --git a/maria/instrument/data/bands/atlast.yml b/maria/instrument/bands/configs/atlast.yml similarity index 53% rename from maria/instrument/data/bands/atlast.yml rename to maria/instrument/bands/configs/atlast.yml index 749c6e6..0a2de4e 100644 --- a/maria/instrument/data/bands/atlast.yml +++ b/maria/instrument/bands/configs/atlast.yml @@ -1,53 +1,53 @@ f027: center: 27 width: 5 - shape: top_hat + shape: gaussian time_constant: 0 efficiency: 1.0 - white_noise: 0 - pink_noise: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz f039: center: 39 width: 5 - shape: top_hat + shape: gaussian time_constant: 0 efficiency: 1.0 - white_noise: 0 - pink_noise: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz f093: center: 93 width: 10 - shape: top_hat + shape: gaussian time_constant: 0 efficiency: 1.0 - white_noise: 0 - pink_noise: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz f150: center: 150 width: 10 - shape: top_hat + shape: gaussian time_constant: 0 efficiency: 1.0 - white_noise: 0 - pink_noise: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz f225: center: 225 width: 30 - shape: top_hat + shape: gaussian time_constant: 0 efficiency: 1.0 - white_noise: 0 - pink_noise: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz f280: center: 280 width: 40 - shape: top_hat + shape: gaussian time_constant: 0 efficiency: 1.0 - white_noise: 0 - pink_noise: 0 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz diff --git a/maria/instrument/data/bands/mustang2.yml b/maria/instrument/bands/configs/mustang2.yml similarity index 57% rename from maria/instrument/data/bands/mustang2.yml rename to maria/instrument/bands/configs/mustang2.yml index e113170..d21c19d 100644 --- a/maria/instrument/data/bands/mustang2.yml +++ b/maria/instrument/bands/configs/mustang2.yml @@ -3,14 +3,14 @@ f093: center: 90 width: 30 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 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 0.5 f093-custom: passband: data/passbands/mustang2/f093.csv 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 + NEP: 1.0 # pW / sqrt(s) + knee: 1.0 # Hz efficiency: 0.2 diff --git a/maria/instrument/data/passbands/mustang2/f093.csv b/maria/instrument/bands/passbands/mustang2/f093.csv similarity index 100% rename from maria/instrument/data/passbands/mustang2/f093.csv rename to maria/instrument/bands/passbands/mustang2/f093.csv diff --git a/maria/instrument/instruments.yml b/maria/instrument/configs.yml similarity index 95% rename from maria/instrument/instruments.yml rename to maria/instrument/configs.yml index 14e7d1d..de0b0a0 100644 --- a/maria/instrument/instruments.yml +++ b/maria/instrument/configs.yml @@ -2,11 +2,8 @@ default: description: A simple test array - documentation: array: - n: 250 field_of_view: 1 - array_shape: circle bands: f150: center: 150 @@ -60,7 +57,7 @@ ALMA: documentation: https://www.eso.org/public/teles-instr/alma/ field_of_view: 0 array: - file: data/arrays/alma/alma.cycle1.total.csv + file: alma/alma.cycle1.total.csv primary_size: 12 diff --git a/maria/instrument/data/bands/abs.yml b/maria/instrument/data/bands/abs.yml deleted file mode 100644 index ea791f4..0000000 --- a/maria/instrument/data/bands/abs.yml +++ /dev/null @@ -1,7 +0,0 @@ -f150: - center: 150 - width: 30 - 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/detectors.py b/maria/instrument/detectors/__init__.py similarity index 97% rename from maria/instrument/detectors.py rename to maria/instrument/detectors/__init__.py index 17b4656..dd99c34 100644 --- a/maria/instrument/detectors.py +++ b/maria/instrument/detectors/__init__.py @@ -4,7 +4,8 @@ import numpy as np import pandas as pd -from .bands import BandList +from ..bands import BandList +from .arrays import generate_array # noqa here, this_filename = os.path.split(__file__) diff --git a/maria/instrument/arrays.py b/maria/instrument/detectors/arrays.py similarity index 87% rename from maria/instrument/arrays.py rename to maria/instrument/detectors/arrays.py index cb60045..ed34393 100644 --- a/maria/instrument/arrays.py +++ b/maria/instrument/detectors/arrays.py @@ -6,8 +6,8 @@ from scipy.spatial import ConvexHull from scipy.spatial.distance import cdist -from .bands import Band -from .beams import compute_angular_fwhm +from ..bands import Band +from ..beams import compute_angular_fwhm here, this_filename = os.path.split(__file__) @@ -21,7 +21,7 @@ def generate_2d_offsets(n, packing="hex", shape="circle", normalize=False): These points are spread such that each is a unit of distance away from its nearest neighbor. """ - n = int(n) + n = int(np.maximum(n, 1)) bigger_n = 2 * n if packing == "square": @@ -59,8 +59,11 @@ def generate_2d_offsets(n, packing="hex", shape="circle", normalize=False): offsets = np.c_[x[subset_index], y[subset_index]] if normalize: - hull_pts = offsets[ConvexHull(offsets).vertices] - offsets /= cdist(hull_pts, hull_pts, metric="euclidean").max() + hull_pts = ( + offsets[ConvexHull(offsets).vertices] if len(offsets) > 16 else offsets + ) + span = cdist(hull_pts, hull_pts, metric="euclidean").max() + offsets /= span if span > 0 else 1.0 return offsets @@ -71,10 +74,11 @@ def generate_2d_offsets_from_diameter( n = np.square(diameter) span = 0 - for i in range(max_iterations): + for _ in range(max_iterations): offsets = generate_2d_offsets(n=n, packing=packing, shape=shape) - ch = ConvexHull(offsets) - hull_pts = ch.points[ch.vertices] + hull_pts = ( + offsets[ConvexHull(offsets).vertices] if len(offsets) > 16 else offsets + ) span = cdist(hull_pts, hull_pts, metric="euclidean").max() n *= np.square(diameter / span) @@ -113,6 +117,15 @@ def generate_array( ] detector_spacing = beam_spacing * np.max(resolutions) + if n is None: + topological_diameter = np.radians(field_of_view) / detector_spacing + if topological_diameter > 1: + offsets = detector_spacing * generate_2d_offsets_from_diameter( + diameter=topological_diameter, packing=array_packing, shape=array_shape + ) + else: + n = 1 # what? + if n is not None: if field_of_view is not None: offsets = np.radians(field_of_view) * generate_2d_offsets( @@ -128,12 +141,6 @@ def generate_array( n=n, packing=array_packing, shape=array_shape ) - else: - topological_diameter = np.radians(field_of_view) / detector_spacing - offsets = detector_spacing * generate_2d_offsets_from_diameter( - diameter=topological_diameter, packing=array_packing, shape=array_shape - ) - baselines = baseline_diameter * generate_2d_offsets( n=len(offsets), packing=baseline_packing, shape=baseline_shape, normalize=True ) diff --git a/maria/instrument/data/arrays/alma/aca.cfg b/maria/instrument/detectors/arrays/alma/aca.cfg similarity index 100% rename from maria/instrument/data/arrays/alma/aca.cfg rename to maria/instrument/detectors/arrays/alma/aca.cfg diff --git a/maria/instrument/data/arrays/alma/alma.cycle1.csv b/maria/instrument/detectors/arrays/alma/alma.cycle1.csv similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle1.csv rename to maria/instrument/detectors/arrays/alma/alma.cycle1.csv diff --git a/maria/instrument/data/arrays/alma/alma.cycle1.f043.csv b/maria/instrument/detectors/arrays/alma/alma.cycle1.f043.csv similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle1.f043.csv rename to maria/instrument/detectors/arrays/alma/alma.cycle1.f043.csv diff --git a/maria/instrument/data/arrays/alma/alma.cycle1.f078.csv b/maria/instrument/detectors/arrays/alma/alma.cycle1.f078.csv similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle1.f078.csv rename to maria/instrument/detectors/arrays/alma/alma.cycle1.f078.csv diff --git a/maria/instrument/data/arrays/alma/alma.cycle1.f100.csv b/maria/instrument/detectors/arrays/alma/alma.cycle1.f100.csv similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle1.f100.csv rename to maria/instrument/detectors/arrays/alma/alma.cycle1.f100.csv diff --git a/maria/instrument/data/arrays/alma/alma.cycle1.f144.csv b/maria/instrument/detectors/arrays/alma/alma.cycle1.f144.csv similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle1.f144.csv rename to maria/instrument/detectors/arrays/alma/alma.cycle1.f144.csv diff --git a/maria/instrument/data/arrays/alma/alma.cycle1.f187.csv b/maria/instrument/detectors/arrays/alma/alma.cycle1.f187.csv similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle1.f187.csv rename to maria/instrument/detectors/arrays/alma/alma.cycle1.f187.csv diff --git a/maria/instrument/data/arrays/alma/alma.cycle1.f243.csv b/maria/instrument/detectors/arrays/alma/alma.cycle1.f243.csv similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle1.f243.csv rename to maria/instrument/detectors/arrays/alma/alma.cycle1.f243.csv diff --git a/maria/instrument/data/arrays/alma/alma.cycle1.f324.csv b/maria/instrument/detectors/arrays/alma/alma.cycle1.f324.csv similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle1.f324.csv rename to maria/instrument/detectors/arrays/alma/alma.cycle1.f324.csv diff --git a/maria/instrument/data/arrays/alma/alma.cycle1.f447.csv b/maria/instrument/detectors/arrays/alma/alma.cycle1.f447.csv similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle1.f447.csv rename to maria/instrument/detectors/arrays/alma/alma.cycle1.f447.csv diff --git a/maria/instrument/data/arrays/alma/alma.cycle1.f661.csv b/maria/instrument/detectors/arrays/alma/alma.cycle1.f661.csv similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle1.f661.csv rename to maria/instrument/detectors/arrays/alma/alma.cycle1.f661.csv diff --git a/maria/instrument/data/arrays/alma/alma.cycle1.f869.csv b/maria/instrument/detectors/arrays/alma/alma.cycle1.f869.csv similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle1.f869.csv rename to maria/instrument/detectors/arrays/alma/alma.cycle1.f869.csv diff --git a/maria/instrument/data/arrays/alma/alma.cycle1.total.csv b/maria/instrument/detectors/arrays/alma/alma.cycle1.total.csv similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle1.total.csv rename to maria/instrument/detectors/arrays/alma/alma.cycle1.total.csv diff --git a/maria/instrument/data/arrays/alma/alma.cycle10.2.cfg b/maria/instrument/detectors/arrays/alma/alma.cycle10.2.cfg similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle10.2.cfg rename to maria/instrument/detectors/arrays/alma/alma.cycle10.2.cfg diff --git a/maria/instrument/data/arrays/alma/alma.cycle10.3.cfg b/maria/instrument/detectors/arrays/alma/alma.cycle10.3.cfg similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle10.3.cfg rename to maria/instrument/detectors/arrays/alma/alma.cycle10.3.cfg diff --git a/maria/instrument/data/arrays/alma/alma.cycle10.4.cfg b/maria/instrument/detectors/arrays/alma/alma.cycle10.4.cfg similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle10.4.cfg rename to maria/instrument/detectors/arrays/alma/alma.cycle10.4.cfg diff --git a/maria/instrument/data/arrays/alma/alma.cycle10.5.cfg b/maria/instrument/detectors/arrays/alma/alma.cycle10.5.cfg similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle10.5.cfg rename to maria/instrument/detectors/arrays/alma/alma.cycle10.5.cfg diff --git a/maria/instrument/data/arrays/alma/alma.cycle10.6.cfg b/maria/instrument/detectors/arrays/alma/alma.cycle10.6.cfg similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle10.6.cfg rename to maria/instrument/detectors/arrays/alma/alma.cycle10.6.cfg diff --git a/maria/instrument/data/arrays/alma/alma.cycle10.7.cfg b/maria/instrument/detectors/arrays/alma/alma.cycle10.7.cfg similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle10.7.cfg rename to maria/instrument/detectors/arrays/alma/alma.cycle10.7.cfg diff --git a/maria/instrument/data/arrays/alma/alma.cycle10.8.cfg b/maria/instrument/detectors/arrays/alma/alma.cycle10.8.cfg similarity index 100% rename from maria/instrument/data/arrays/alma/alma.cycle10.8.cfg rename to maria/instrument/detectors/arrays/alma/alma.cycle10.8.cfg diff --git a/maria/utils/io.py b/maria/io.py similarity index 66% rename from maria/utils/io.py rename to maria/io.py index cc98c51..7b1e100 100644 --- a/maria/utils/io.py +++ b/maria/io.py @@ -4,10 +4,13 @@ from collections.abc import Mapping from datetime import datetime +import astropy as ap import h5py +import pandas as pd import pytz import requests import yaml +from tqdm import tqdm def flatten_config(m: dict, prefix: str = ""): @@ -42,27 +45,43 @@ def read_yaml(path: str): return res if res is not None else {} -def cache_is_ok(path: str, CACHE_MAX_AGE: float = 86400): +def cache_is_ok(path: str, max_age: float = 30 * 86400, verbose: bool = False): """ Check if we need to reload the cache. """ if not os.path.exists(path): + if verbose: + print(f"Cached file at {path} does not exist.") return False cache_age = ttime.time() - os.path.getmtime(path) - if cache_age > CACHE_MAX_AGE: + if cache_age > max_age: + if verbose: + print(f"Cached file at {path} is too old.") return False - extension = path.split(".")[-1] + if not test_file(path): + if verbose: + print(f"Could not open cached file at {path}.") + return False + + return True + +def test_file(path) -> bool: + ext = path.split(".")[-1] try: - if extension == "h5": + if ext in ["h5"]: with h5py.File(path, "r") as f: f.keys() - else: + elif ext in ["csv"]: + pd.read_csv(path) + elif ext in ["txt", "dat"]: with open(path, "r") as f: f.read() + elif ext in ["fits"]: + ap.io.fits.open(path) except Exception: return False @@ -72,9 +91,10 @@ def cache_is_ok(path: str, CACHE_MAX_AGE: float = 86400): def fetch_cache( source_url: str, cache_path: str, - CACHE_MAX_AGE: float = 7 * 86400, + max_age: float = 30 * 86400, refresh: bool = False, chunk_size: int = 8192, + verbose: bool = False, ): """ Download the cache if needed @@ -86,18 +106,24 @@ def fetch_cache( print(f"created cache at {cache_dir}") os.makedirs(cache_dir, exist_ok=True) - if (not cache_is_ok(cache_path, CACHE_MAX_AGE=CACHE_MAX_AGE)) or refresh: - print(f"updating cache from {source_url}") - + if (not cache_is_ok(cache_path, max_age=max_age, verbose=verbose)) or refresh: with requests.get(source_url, stream=True) as r: r.raise_for_status() with open(cache_path, "wb") as f: - for chunk in r.iter_content(chunk_size=chunk_size): + chunks = tqdm( + r.iter_content(chunk_size=chunk_size), + desc=f"Updating cache from {source_url}", + disable=not verbose, + ) + for chunk in chunks: f.write(chunk) cache_size = os.path.getsize(cache_path) print(f"downloaded data ({1e-6 * cache_size:.01f} MB) to {cache_path}") + if not test_file(cache_path): + raise RuntimeError("Could not open cached file.") + def datetime_handler(time): """ diff --git a/maria/map/__init__.py b/maria/map/__init__.py index 872a260..d3d522a 100644 --- a/maria/map/__init__.py +++ b/maria/map/__init__.py @@ -1,182 +1,15 @@ import os -from dataclasses import dataclass, fields -from operator import attrgetter -from typing import List, Tuple import astropy as ap -import matplotlib as mpl -import matplotlib.pyplot as plt import numpy as np import scipy as sp -from astropy.io import fits -from astropy.wcs import WCS -from matplotlib.colors import ListedColormap from tqdm import tqdm -from .. import utils from ..instrument import beams +from .map import Map here, this_filename = os.path.split(__file__) -# from https://gist.github.com/zonca/6515744 -cmb_cmap = ListedColormap( - np.loadtxt(f"{here}/Planck_Parchment_RGB.txt") / 255.0, name="cmb" -) -cmb_cmap.set_bad("white") - -mpl.colormaps.register(cmb_cmap) - - -@dataclass -class Map: - """ - We define height and width, which determines the shape of the - - This means that there might be non-square pixels - """ - - freqs: List[float] - center: Tuple[float, float] - width: float = 1 - height: float = 1 - degrees: bool = True - inbright: float = 1 - header: ap.io.fits.header.Header = None - frame: str = "ra_dec" - units: str = "K_RJ" - data: np.array = None # 3D instrument - weight: np.array = None - - def __post_init__(self): - assert len(self.freqs) == self.n_freqs - - assert self.data is not None - - self.width = self.res * self.n_x - self.height = self.res * self.n_y - - if self.weight is None: - self.weight = np.ones(self.data.shape) - - def __repr__(self): - nodef_f_vals = ( - (f.name, attrgetter(f.name)(self)) - for f in fields(self) - if f.name not in ["data", "header"] - ) - - nodef_f_repr = ", ".join(f"{name}={value}" for name, value in nodef_f_vals) - return f"{self.__class__.__name__}({nodef_f_repr})" - - @property - def res(self): - """ - TODO: don't do this - """ - return self.x_res - - @property - def x_res(self): - return self.width / self.n_x - - @property - def y_res(self): - return self.height / self.n_y - - @property - def n_freqs(self): - return self.data.shape[0] - - @property - def n_x(self): - return self.data.shape[2] - - @property - def n_y(self): - return self.data.shape[1] - - @property - def x_side(self): - x = self.res * np.arange(self.n_x) - return x - x.mean() - - @property - def y_side(self): - y = self.res * np.arange(self.n_y) - return y - y.mean() - - def plot( - self, cmap="cmb", rel_vmin=0.001, rel_vmax=0.999, units="degrees", **kwargs - ): - for i_freq, freq in enumerate(self.freqs): - header = fits.header.Header() - - header["RESTFRQ"] = freq - - res_degrees = self.res if self.degrees else np.degrees(self.res) - center_degrees = self.center if self.degrees else np.degrees(self.center) - - header["CDELT1"] = res_degrees # degree - header["CDELT2"] = res_degrees # degree - - header["CRPIX1"] = self.n_x / 2 - header["CRPIX2"] = self.n_y / 2 - - header["CTYPE1"] = "RA---SIN" - header["CUNIT1"] = "deg " - header["CTYPE2"] = "DEC--SIN" - header["CUNIT2"] = "deg " - header["RADESYS"] = "FK5 " - - header["CRVAL1"] = center_degrees[0] - header["CRVAL2"] = center_degrees[1] - wcs_input = WCS(header, naxis=2) # noqa F401 - - fig = plt.figure() - - ax = fig.add_subplot(1, 1, 1) # , projection=wcs_input) - - ax.set_title(f"{freq} GHz") - - map_extent_radians = [ - -self.width / 2, - self.width / 2, - -self.height / 2, - self.height / 2, - ] - if self.degrees: - map_extent_radians = np.radians(map_extent_radians) - - if units == "degrees": - map_extent = np.degrees(map_extent_radians) - if units == "arcmin": - map_extent = 60 * np.degrees(map_extent_radians) - if units == "arcsec": - map_extent = 3600 * np.degrees(map_extent_radians) - - d = self.data.ravel() - w = self.weight.ravel() - - sort = np.argsort(d) - d, w = d[sort], w[sort] - - vmin, vmax = np.interp([rel_vmin, rel_vmax], np.cumsum(w) / np.sum(w), d) - - map_im = ax.imshow( - self.data.T, - cmap=cmap, - interpolation="none", - extent=map_extent, - vmin=vmin, - vmax=vmax, - ) - - ax.set_xlabel(rf"$\Delta\,\theta_x$ [{units}]") - ax.set_ylabel(rf"$\Delta\,\theta_y$ [{units}]") - - cbar = fig.colorbar(map_im, ax=ax, shrink=1.0) - cbar.set_label("mJy km/s/pixel") - class MapMixin: """ @@ -185,75 +18,17 @@ class MapMixin: TODO: add errors """ - def _initialize_map(self): - if not self.map_file: - return - - self.map_file = self.map_file - hudl = ap.io.fits.open(self.map_file) - - map_data = hudl[0].data - if map_data.ndim < 2: - raise ValueError() - elif map_data.ndim == 2: - map_data = map_data[None] - - map_n_freqs, map_n_y, map_n_x = map_data.shape - - if map_n_freqs != len(self.map_freqs): - raise ValueError() - - map_width = self.map_res * map_n_x - map_height = self.map_res * map_n_y - - self.raw_map_data = map_data.copy() - - res_degrees = self.map_res if self.degrees else np.degrees(self.map_res) - - if self.map_units == "Jy/pixel": - for i, nu in enumerate(self.map_freqs): - map_data[i] = map_data[i] / utils.units.KbrightToJyPix( - 1e9 * nu, res_degrees, res_degrees - ) - - self.map_data = map_data - - self.map = Map( - data=map_data, - header=hudl[0].header, - freqs=np.atleast_1d(self.map_freqs), - width=np.radians(map_width) if self.degrees else map_width, - height=np.radians(map_height) if self.degrees else map_height, - center=np.radians(self.map_center) if self.degrees else map_height, - degrees=False, - frame=self.pointing_frame, - inbright=self.map_inbright, - units=self.map_units, - ) - - self.map.header["HISTORY"] = "History_input_adjustments" - self.map.header["comment"] = "Changed input CDELT1 and CDELT2" - self.map.header["comment"] = ( - "Changed surface brightness units to " + self.map.units - ) - self.map.header["comment"] = "Repositioned the map on the sky" - - if self.map.inbright is not None: - self.map.data *= self.map.inbright / np.nanmax(self.map.data) - self.map.header["comment"] = "Amplitude is rescaled." - def _run(self, **kwargs): self.sample_maps() def _sample_maps(self): - dx, dy = self.coords.offsets(frame=self.map_frame, center=self.map.center) + dx, dy = self.coords.offsets(frame=self.map.frame, center=self.map.center) self.data["map"] = np.zeros((dx.shape)) - freqs = tqdm(self.map.freqs) if self.verbose else self.map.freqs - for i, nu in enumerate(freqs): - if self.verbose: - freqs.set_description("Sampling input map") + 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( @@ -288,4 +63,54 @@ def _sample_maps(self): if hasattr(self, "atmospheric_transmission"): self.data["map"] *= self.atmospheric_transmission - self.data["map"] /= self.instrument.dets.pW_to_KRJ[:, None] + self.data["map"] *= 1e12 * self.instrument.dets.dP_dTRJ[:, None] + + +def from_fits( + filename: str, + resolution: float, + frequency: float, + center: tuple, + index: int = 0, + units: str = "K_RJ", + degrees: bool = True, + frame: str = "az_el", +) -> Map: + if not os.path.exists(filename): + raise FileNotFoundError(filename) + + hudl = ap.io.fits.open(filename) + + frequency = np.atleast_1d(frequency) + + map_data = hudl[index].data + if map_data.ndim < 2: + raise ValueError("Map should have at least 2 dimensions.") + elif map_data.ndim == 2: + map_data = map_data[None] + + *_, map_n_y, map_n_x = map_data.shape + + map_width = resolution * map_n_x + map_height = resolution * map_n_y + + # res_degrees = resolution if degrees else np.degrees(resolution) + + # if self.map_units == "Jy/pixel": + # for i, nu in enumerate(self.map_freqs): + # map_data[i] = map_data[i] / utils.units.KbrightToJyPix( + # 1e9 * nu, res_degrees, res_degrees + # ) + + return Map( + data=map_data, + header=hudl[0].header, + frequency=frequency, + width=np.radians(map_width) if degrees else map_width, + height=np.radians(map_height) if degrees else map_height, + center=np.radians(center) if degrees else center, + degrees=False, + frame=frame, + # inbright=self.map_inbright, + units=units, + ) diff --git a/maria/map/map.py b/maria/map/map.py new file mode 100644 index 0000000..fa058ec --- /dev/null +++ b/maria/map/map.py @@ -0,0 +1,191 @@ +import os +from dataclasses import dataclass, fields +from operator import attrgetter +from typing import List, Tuple + +import astropy as ap +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np +from astropy.io import fits +from astropy.wcs import WCS +from matplotlib.colors import ListedColormap + +here, this_filename = os.path.split(__file__) + +# from https://gist.github.com/zonca/6515744 +cmb_cmap = ListedColormap( + np.loadtxt(f"{here}/Planck_Parchment_RGB.txt") / 255.0, name="cmb" +) +cmb_cmap.set_bad("white") + +mpl.colormaps.register(cmb_cmap) + + +@dataclass +class Map: + """ + We define height and width, which determines the shape of the + + This means that there might be non-square pixels + """ + + data: np.array + frequency: List[float] = 100.0 + center: Tuple[float, float] = () + width: float = 1 + height: float = 1 + degrees: bool = True + inbright: float = 1 + header: ap.io.fits.header.Header = None + frame: str = "ra_dec" + units: str = "K_RJ" + weight: np.array = None + + def __post_init__(self): + self.frequency = np.atleast_1d(self.frequency) + + if len(self.frequency) != self.data.shape[-3]: + raise ValueError( + f"Number of supplied frequencies ({len(self.frequency)}) does not match the " + f"frequency dimension of the supplied map ({self.data.shape[-2]})." + ) + + self.header = ap.io.fits.header.Header() + + self.header["CDELT1"] = np.degrees(self.x_res) # degree + self.header["CDELT2"] = np.degrees(self.y_res) # degree + self.header["CTYPE1"] = "RA---SIN" + self.header["CUNIT1"] = "deg " + self.header["CTYPE2"] = "DEC--SIN" + self.header["CUNIT2"] = "deg " + self.header["CRPIX1"] = self.data.shape[-1] + self.header["CRPIX2"] = self.data.shape[-2] + self.header["CRVAL1"] = self.center[0] + self.header["CRVAL2"] = self.center[1] + self.header["RADESYS"] = "FK5 " + + self.width = self.res * self.n_x + self.height = self.res * self.n_y + + if self.weight is None: + self.weight = np.ones(self.data.shape) + + def __repr__(self): + nodef_f_vals = ( + (f.name, attrgetter(f.name)(self)) + for f in fields(self) + if f.name not in ["data", "weight", "header"] + ) + + nodef_f_repr = ", ".join(f"{name}={value}" for name, value in nodef_f_vals) + return f"{self.__class__.__name__}({nodef_f_repr})" + + @property + def res(self): + """ + TODO: don't do this + """ + return self.x_res + + @property + def x_res(self): + return self.width / self.n_x + + @property + def y_res(self): + return self.height / self.n_y + + @property + def n_freqs(self): + return self.data.shape[0] + + @property + def n_x(self): + return self.data.shape[2] + + @property + def n_y(self): + return self.data.shape[1] + + @property + def x_side(self): + x = self.res * np.arange(self.n_x) + return x - x.mean() + + @property + def y_side(self): + y = self.res * np.arange(self.n_y) + return y - y.mean() + + def plot( + self, cmap="cmb", rel_vmin=0.001, rel_vmax=0.999, units="degrees", **kwargs + ): + for i_freq, freq in enumerate(self.frequency): + header = fits.header.Header() + + header["RESTFRQ"] = freq + + res_degrees = self.res if self.degrees else np.degrees(self.res) + center_degrees = self.center if self.degrees else np.degrees(self.center) + + header["CDELT1"] = res_degrees # degree + header["CDELT2"] = res_degrees # degree + + header["CRPIX1"] = self.n_x / 2 + header["CRPIX2"] = self.n_y / 2 + + header["CTYPE1"] = "RA---SIN" + header["CUNIT1"] = "deg " + header["CTYPE2"] = "DEC--SIN" + header["CUNIT2"] = "deg " + header["RADESYS"] = "FK5 " + + header["CRVAL1"] = center_degrees[0] + header["CRVAL2"] = center_degrees[1] + wcs_input = WCS(header, naxis=2) # noqa F401 + + fig = plt.figure() + + ax = fig.add_subplot(1, 1, 1) # , projection=wcs_input) + + ax.set_title(f"{freq} GHz") + + map_extent_radians = [ + -self.width / 2, + self.width / 2, + -self.height / 2, + self.height / 2, + ] + if self.degrees: + map_extent_radians = np.radians(map_extent_radians) + + if units == "degrees": + map_extent = np.degrees(map_extent_radians) + if units == "arcmin": + map_extent = 60 * np.degrees(map_extent_radians) + if units == "arcsec": + map_extent = 3600 * np.degrees(map_extent_radians) + + d = self.data.ravel() + w = self.weight.ravel() + + sort = np.argsort(d) + d, w = d[sort], w[sort] + + vmin, vmax = np.interp([rel_vmin, rel_vmax], np.cumsum(w) / np.sum(w), d) + + map_im = ax.imshow( + self.data.T, + cmap=cmap, + interpolation="none", + extent=map_extent, + vmin=vmin, + vmax=vmax, + ) + + ax.set_xlabel(rf"$\Delta\,\theta_x$ [{units}]") + ax.set_ylabel(rf"$\Delta\,\theta_y$ [{units}]") + + cbar = fig.colorbar(map_im, ax=ax, shrink=1.0) + cbar.set_label("mJy km/s/pixel") diff --git a/maria/map/mappers.py b/maria/map/mappers.py index 8c8d263..9752673 100644 --- a/maria/map/mappers.py +++ b/maria/map/mappers.py @@ -4,10 +4,10 @@ import numpy as np import scipy as sp from astropy.io import fits +from todder import TOD -from .. import utils -from ..tod import TOD -from . import Map +from .. import units, utils +from .map import Map np.seterr(invalid="ignore") @@ -95,16 +95,16 @@ def save_maps(self, filepath): else: raise ValueError(f"Units {self.map.units} not implemented.") - save_maps = np.zeros((len(self.map.freqs), self.n_x, self.n_y)) + save_maps = np.zeros((len(self.map.frequency), self.n_x, self.n_y)) for i, key in enumerate(self.band_data.keys()): self.header["CRVAL3"] = self.band_data[key]["band_center"] * 1e9 - self.header["CDELT3"] = self.band_data[key]["band_centerwidth"] * 1e9 + self.header["CDELT3"] = self.band_data[key]["band_width"] * 1e9 save_maps[i] = self.map.data[i] if self.map.units == "Jy/pixel": - save_maps[i] *= utils.units.KbrightToJyPix( + save_maps[i] *= units.KbrightToJyPix( self.header["CRVAL3"], self.header["CDELT1"], self.header["CDELT2"] ) @@ -235,7 +235,7 @@ def run(self): self.raw_map_cnts[band] += map_cnt self.band_data[band]["band_center"] = tod.dets.band_center.mean() - self.band_data[band]["band_centerwidth"] = 30 + self.band_data[band]["band_width"] = 30 band_map_numer = self.raw_map_sums[band].copy() band_map_denom = self.raw_map_cnts[band].copy() @@ -264,7 +264,9 @@ def run(self): self.map = Map( data=self.map_data, weight=self.map_weight, - freqs=np.array([self.band_data[band]["band_center"] for band in self.band]), + frequency=np.array( + [self.band_data[band]["band_center"] for band in self.band] + ), width=np.degrees(self.width) if self.degrees else self.width, height=np.degrees(self.height) if self.degrees else self.height, center=np.degrees(self.center) if self.degrees else self.center, diff --git a/maria/noise/__init__.py b/maria/noise/__init__.py index 357c938..e10c03a 100644 --- a/maria/noise/__init__.py +++ b/maria/noise/__init__.py @@ -1,4 +1,5 @@ import numpy as np +from todder.sim.noise import generate_noise_with_knee from ..sim.base import BaseSimulation @@ -12,54 +13,56 @@ 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_name=band.name).index + band_mask = self.instrument.dets.band_name == band.name - if band.white_noise > 0: - self.data["noise"][band_index] += ( - np.sqrt(self.plan.sample_rate) - * band.white_noise - * np.random.standard_normal( - size=(len(band_index), self.plan.n_time) - ) - / band.pW_to_KRJ - ) + self.data["noise"][band_mask] = generate_noise_with_knee( + self.plan.time, n=band_mask.sum(), NEP=band.NEP, knee=band.knee + ) - 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, - ) - / band.pW_to_KRJ - ) + # # if band.white_noise > 0: + # # self.data["noise"][band_index] += ( + # # np.sqrt(self.plan.sample_rate) + # # * np.sqrt(band.white_noise) + # # * np.random.standard_normal( + # # size=(len(band_index), self.plan.n_time) + # # ) + # # ) - def _spectrum_noise(self, spectrum_func, size, dt, amp=2.0): - """ - make noise with a certain spectral density - """ - freqs = np.fft.rfftfreq( - size, dt - ) # real-fft frequencies (not the negative ones) - spectrum = np.zeros_like( - freqs, dtype="complex" - ) # make complex numbers for spectrum - spectrum[1:] = spectrum_func( - freqs[1:] - ) # get s pectrum amplitude for all frequencies except f=0 - phases = np.random.uniform( - 0, 2 * np.pi, len(freqs) - 1 - ) # random phases for all frequencies except f=0 - spectrum[1:] *= np.exp(1j * phases) - noise = np.fft.irfft(spectrum) # return the reverse fourier transform - return noise + # # if band.pink_noise > 0: + # # for i in band_index: + # # self.data["noise"][i] += ( + # # np.sqrt(band.pink_noise) + # # * self._spectrum_noise( + # # spectrum_func=self._pink_spectrum, + # # size=int(self.plan.n_time), + # # dt=self.plan.dt, + # # ) + # # ) - def _pink_spectrum(self, f, f_min=0, f_max=np.inf, amp=1.0): - s = (f / amp) ** -0.5 - s[np.logical_or(f < f_min, f > f_max)] = 0 # apply band pass - return s + # def _spectrum_noise(self, spectrum_func, size, dt, amp=2.0): + # """ + # make noise with a certain spectral density + # """ + # freqs = np.fft.rfftfreq( + # size, dt + # ) # real-fft frequencies (not the negative ones) + # spectrum = np.zeros_like( + # freqs, dtype="complex" + # ) # make complex numbers for spectrum + # spectrum[1:] = spectrum_func( + # freqs[1:] + # ) # get s pectrum amplitude for all frequencies except f=0 + # phases = np.random.uniform( + # 0, 2 * np.pi, len(freqs) - 1 + # ) # random phases for all frequencies except f=0 + # spectrum[1:] *= np.exp(1j * phases) + # noise = np.fft.irfft(spectrum) # return the reverse fourier transform + # return noise + + # def _pink_spectrum(self, f, f_min=0, f_max=np.inf, amp=1.0): + # s = (f / amp) ** -0.5 + # s[np.logical_or(f < f_min, f > f_max)] = 0 # apply band pass + # return s class NoiseSimulation(NoiseMixin, BaseSimulation): diff --git a/maria/plan/__init__.py b/maria/plan/__init__.py index 412db45..eec93b6 100644 --- a/maria/plan/__init__.py +++ b/maria/plan/__init__.py @@ -8,8 +8,9 @@ import numpy as np import pandas as pd import pytz +from todder import coords -from .. import coords, utils +from ..io import datetime_handler, read_yaml from . import patterns MAX_VELOCITY_WARN = 10 # in deg/s @@ -20,7 +21,7 @@ here, this_filename = os.path.split(__file__) -plan_configs = utils.io.read_yaml(f"{here}/plans.yml") +plan_configs = read_yaml(f"{here}/plans.yml") plan_params = set() for key, config in plan_configs.items(): plan_params |= set(config.keys()) @@ -36,7 +37,7 @@ def __init__(self, invalid_plan): ) -def get_plan_config(scan_pattern="stare", **kwargs): +def get_plan_config(scan_pattern="daisy", **kwargs): if scan_pattern not in plan_configs.keys(): raise UnsupportedPlanError(scan_pattern) plan_config = plan_configs[scan_pattern].copy() @@ -45,7 +46,7 @@ def get_plan_config(scan_pattern="stare", **kwargs): return plan_config -def get_plan(scan_pattern="stare", **kwargs): +def get_plan(scan_pattern="daisy", **kwargs): plan_config = get_plan_config(scan_pattern, **kwargs) return Plan(**plan_config) @@ -73,10 +74,10 @@ class Plan: """ description: str = "" - start_time: str = "2022-02-10T06:00:00" + start_time: str or int = "2022-02-10T06:00:00" duration: float = 60.0 sample_rate: float = 20.0 - pointing_frame: str = "ra_dec" + frame: str = "ra_dec" degrees: bool = True scan_center: Tuple[float, float] = (4, 10.5) scan_pattern: str = "daisy_miss_center" @@ -104,7 +105,7 @@ def __post_init__(self): if not hasattr(self, "start_time"): self.start_time = datetime.now().timestamp() - self.start_datetime = utils.io.datetime_handler(self.start_time) + self.start_datetime = datetime_handler(self.start_time) self.end_datetime = self.start_datetime + timedelta(seconds=self.duration) self.time_min = self.start_datetime.timestamp() @@ -168,9 +169,9 @@ def __post_init__(self): self.phi, self.theta = coords.dx_dy_to_phi_theta( *self.scan_offsets_radians, *self.scan_center_radians ) - if self.pointing_frame == "ra_dec": + if self.frame == "ra_dec": self.ra, self.dec = self.phi, self.theta - elif self.pointing_frame == "az_el": + elif self.frame == "az_el": self.az, self.el = self.phi, self.theta else: raise ValueError("Not a valid pointing frame!") @@ -200,8 +201,8 @@ def plot(self): pointing_units = "deg." if self.degrees else "rad." label = ( - f"""{coords.frames[self.pointing_frame]['phi_name']} = {center_phi} {pointing_units}""" - f"""{coords.frames[self.pointing_frame]['theta_name']} = {center_theta} {pointing_units}""" + f"""{coords.frames[self.frame]['phi_name']} = {center_phi} {pointing_units}""" + f"""{coords.frames[self.frame]['theta_name']} = {center_theta} {pointing_units}""" ) ax.plot(dx, dy, lw=5e-1) diff --git a/maria/plan/plans.yml b/maria/plan/plans.yml index a69ed3d..1ed9b7c 100644 --- a/maria/plan/plans.yml +++ b/maria/plan/plans.yml @@ -2,7 +2,7 @@ stare: description: A stare. duration: 60 scan_pattern: stare - pointing_frame: az_el + frame: az_el scan_center: [0, 45] degrees: True sample_rate: 20 @@ -10,7 +10,7 @@ stare: daisy: description: A daisy scan. - pointing_frame: ra_dec + frame: ra_dec degrees: True scan_pattern: daisy scan_center: [150., 10] @@ -28,7 +28,7 @@ back_and_forth: duration: 60 scan_pattern: back_and_forth scan_center: [45, 45] - pointing_frame: az_el + frame: az_el degrees: True sample_rate: 20 scan_options: @@ -38,7 +38,7 @@ back_and_forth: double_circle: description: 'A double circle scan.' - pointing_frame: ra_dec + frame: ra_dec duration: 60 sample_rate: 50 scan_pattern: double_circle diff --git a/maria/sim/__init__.py b/maria/sim/__init__.py index ea471c7..988851c 100644 --- a/maria/sim/__init__.py +++ b/maria/sim/__init__.py @@ -1,20 +1,23 @@ import os +import time as ttime -from maria.instrument import Instrument -from maria.plan import Plan -from maria.site import Site +from tqdm import tqdm -from ..atmosphere import Atmosphere, AtmosphereMixin -from ..cmb import CMBMixin -from ..map import MapMixin +from ..atmosphere import Atmosphere +from ..atmosphere.sim import AtmosphereMixin +from ..cmb import CMB, CMBMixin, generate_cmb, get_cmb +from ..instrument import Instrument +from ..map import Map, MapMixin from ..noise import NoiseMixin -from .base import BaseSimulation, master_params, parse_sim_kwargs +from ..plan import Plan +from ..site import Site +from .base import BaseSimulation here, this_filename = os.path.split(__file__) class Simulation(BaseSimulation, AtmosphereMixin, CMBMixin, MapMixin, NoiseMixin): - """A simulation! This is what users should touch, primarily.""" + """A simulation of a telescope. This is what users should touch, primarily.""" @classmethod def from_config(cls, config: dict = {}, **params): @@ -22,60 +25,87 @@ def from_config(cls, config: dict = {}, **params): def __init__( self, - instrument: str or Instrument = "MUSTANG-2", - plan: str or Plan = "stare", - site: str or Site = "hoagie_haven", + instrument: Instrument or str = "default", + plan: Plan or str = "daisy", + site: Site or str = "hoagie_haven", + map: Map = None, + atmosphere: str = None, + cmb: CMB or str = None, + noise: bool = True, + atmosphere_kwargs: dict = {}, + cmb_kwargs: dict = {}, + noise_kwargs: dict = {}, verbose: bool = True, - **kwargs, ): - self.parsed_sim_kwargs = parse_sim_kwargs(kwargs, master_params, strict=True) + # self.parsed_sim_kwargs = parse_sim_kwargs(kwargs, master_params, strict=True) + + ref_time = ttime.time() super().__init__( instrument, plan, site, verbose=verbose, - **self.parsed_sim_kwargs["instrument"], - **self.parsed_sim_kwargs["plan"], - **self.parsed_sim_kwargs["site"], + # **self.parsed_sim_kwargs["instrument"], + # **self.parsed_sim_kwargs["plan"], + # **self + # .parsed_sim_kwargs["site"], ) - self.noise = True + duration = ttime.time() - ref_time + ref_time = ttime.time() - self.params = {} + if self.verbose: + print(f"Initialized base in {int(1e3 * duration)} ms.") - for sub_type, sub_master_params in master_params.items(): - self.params[sub_type] = {} - if sub_type in ["instrument", "site", "plan"]: - sub_type_dataclass = getattr(self, sub_type) - for k in sub_type_dataclass.__dataclass_fields__.keys(): - v = getattr(sub_type_dataclass, k) - setattr(self, k, v) - self.params[sub_type][k] = v - else: - for k, v in sub_master_params.items(): - setattr(self, k, kwargs.get(k, v)) - self.params[sub_type][k] = v + self.noise = noise + + # self.noise = True + + # self.params = {} - if self.map_file: - if not os.path.exists(self.map_file): - raise FileNotFoundError(self.map_file) - self._initialize_map() + # for sub_type, sub_master_params in master_params.items(): + # self.params[sub_type] = {} + # if sub_type in ["instrument", "site", "plan"]: + # sub_type_dataclass = getattr(self, sub_type) + # for k in sub_type_dataclass.__dataclass_fields__.keys(): + # v = getattr(sub_type_dataclass, k) + # setattr(self, k, v) + # self.params[sub_type][k] = v + # else: + # for k, v in sub_master_params.items(): + # setattr(self, k, kwargs.get(k, v)) + # self.params[sub_type][k] = v - if self.atmosphere_model: - weather_override = {k: v for k, v in {"pwv": self.pwv}.items() if v} + if map: + self.map = map + ... + + if atmosphere: + weather_kwargs = ( + atmosphere_kwargs.pop("weather") + if "weather" in atmosphere_kwargs + else {} + ) self.atmosphere = Atmosphere( t=self.plan.time.mean(), region=self.site.region, altitude=self.site.altitude, - weather_override=weather_override, + weather_kwargs=weather_kwargs, ) - self._initialize_atmosphere() + if atmosphere == "2d": + self._initialize_2d_atmosphere(**atmosphere_kwargs) - if self.cmb_source: - self._initialize_cmb(source=self.cmb_source) + if cmb: + if cmb in ["spectrum", "power_spectrum", "generate", "generated"]: + for _ in tqdm(range(1), desc="Generating CMB"): + self.cmb = generate_cmb(verbose=self.verbose, **cmb_kwargs) + elif cmb in ["real", "planck"]: + self.cmb = get_cmb(**cmb_kwargs) + else: + raise ValueError(f"Invalid value for cmb: '{cmb}'.") def _run(self): # number of bands are lost here @@ -99,7 +129,7 @@ def _run(self): 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] + # self.data[k] *= self.instrument.dets.pW_to_KRJ[:, None] def __repr__(self): object_reprs = [ diff --git a/maria/sim/base.py b/maria/sim/base.py index 9a71475..92d0ed6 100644 --- a/maria/sim/base.py +++ b/maria/sim/base.py @@ -2,13 +2,13 @@ import dask.array as da import numpy as np +from todder import TOD +from todder.coords import Coordinates, dx_dy_to_phi_theta -from .. import utils -from ..coords import Coordinates, dx_dy_to_phi_theta from ..instrument import Instrument, get_instrument +from ..io import read_yaml from ..plan import Plan, get_plan from ..site import Site, get_site -from ..tod import TOD here, this_filename = os.path.split(__file__) @@ -20,7 +20,7 @@ def __init__(self, invalid_keys): ) -master_params = utils.io.read_yaml(f"{here}/params.yml") +master_params = read_yaml(f"{here}/params.yml") def parse_sim_kwargs(kwargs, master_kwargs, strict=False): @@ -90,11 +90,11 @@ def __init__( self.calibration = np.ones((self.instrument.dets.n, self.plan.n_time)) self.boresight = Coordinates( - self.plan.time - self.plan.time.min(), - self.plan.phi, - self.plan.theta, + time=self.plan.time - self.plan.time.min(), + phi=self.plan.phi, + theta=self.plan.theta, location=self.site.earth_location, - frame=self.plan.pointing_frame, + frame=self.plan.frame, time_offset=self.plan.time.min(), ) @@ -119,9 +119,9 @@ def __init__( ) self.coords = Coordinates( - self.boresight.time, - det_az, - det_el, + time=self.boresight.time, + phi=det_az, + theta=det_el, location=self.site.earth_location, frame="az_el", time_offset=self.boresight.time_offset, diff --git a/maria/sim/params.yml b/maria/sim/params.yml index 3407d3d..43463eb 100644 --- a/maria/sim/params.yml +++ b/maria/sim/params.yml @@ -25,7 +25,7 @@ plan: # defaults to a 45 degree stare due north scan_options: mapping site: # default to the ALMA site - site_description: '' + description: '' region: 'chajnantor' latitude: -23.0294 longitude: -67.7548 diff --git a/maria/site/__init__.py b/maria/site/__init__.py index 95059ce..f57997d 100644 --- a/maria/site/__init__.py +++ b/maria/site/__init__.py @@ -4,17 +4,17 @@ import pandas as pd from astropy.coordinates import EarthLocation -from .. import utils +from ..io import read_yaml here, this_filename = os.path.split(__file__) -SITE_CONFIGS = utils.io.read_yaml(f"{here}/sites.yml") +SITE_CONFIGS = read_yaml(f"{here}/sites.yml") SITE_PARAMS = set() for key, config in SITE_CONFIGS.items(): SITE_PARAMS |= set(config.keys()) SITE_DISPLAY_COLUMNS = [ - "site_description", + "description", "region", "latitude", "longitude", @@ -51,7 +51,7 @@ def get_location(site_name): ) -def get_site_config(site_name="default", **kwargs): +def get_site_config(site_name="hoagie_haven", **kwargs): if site_name not in SITE_CONFIGS.keys(): raise InvalidSiteError(site_name) SITE_CONFIG = SITE_CONFIGS[site_name].copy() @@ -60,13 +60,13 @@ def get_site_config(site_name="default", **kwargs): return SITE_CONFIG -def get_site(site_name="default", **kwargs): - return Site(**get_site_config(site_name, **kwargs)) +def get_site(site_name="hoagie_haven", **kwargs): + return Site(**get_site_config(site_name=site_name, **kwargs)) @dataclass class Site: - site_description: str = "" + description: str = "" region: str = "princeton" altitude: float = None # in meters seasonal: bool = True diff --git a/maria/site/sites.yml b/maria/site/sites.yml index 730303e..bb89fb3 100644 --- a/maria/site/sites.yml +++ b/maria/site/sites.yml @@ -1,5 +1,5 @@ hoagie_haven: - site_description: Hoagie Haven + description: Hoagie Haven region: princeton latitude: 40.3522 longitude: -74.6519 @@ -7,7 +7,7 @@ hoagie_haven: site_documentation: http://hoagiehaven.com/ amundsen_scott: - site_description: Amundsen–Scott South Pole Station (BICEP, SPT) + description: Amundsen–Scott South Pole Station (BICEP, SPT) region: south_pole latitude: 0.000 longitude: -90.000 @@ -15,7 +15,7 @@ amundsen_scott: site_documentation: https://www.nsf.gov/geo/opp/support/southp.jsp cerro_chajnantor: - site_description: Cerro Chajnantor (CCAT, TAO) + description: Cerro Chajnantor (CCAT, TAO) region: chajnantor latitude: -22.9862 longitude: -67.7422 @@ -23,7 +23,7 @@ cerro_chajnantor: site_documentation: https://act.princeton.edu/ cerro_toco: - site_description: Cerro Toco (ACT, POLARBEAR, SO) + description: Cerro Toco (ACT, POLARBEAR, SO) region: chajnantor latitude: -22.9586 longitude: -67.7878 @@ -31,7 +31,7 @@ cerro_toco: site_documentation: https://act.princeton.edu/ effelsberg: - site_description: Effelsberg (ERT) + description: Effelsberg (ERT) region: effelsberg latitude: 50.5247 longitude: 6.8828 @@ -39,7 +39,7 @@ effelsberg: site_documentation: https://www.mpifr-bonn.mpg.de/en/effelsberg green_bank: - site_description: Green Bank Observatory (MUSTANG-2) + description: Green Bank Observatory (MUSTANG-2) region: green_bank latitude: 38.4331 longitude: -79.8398 @@ -47,7 +47,7 @@ green_bank: site_documentation: https://greenbankobservatory.org llano_de_chajnantor: - site_description: Llano de Chajnantor (ALMA) + description: Llano de Chajnantor (ALMA) region: chajnantor latitude: -23.0294 longitude: -67.7548 @@ -55,7 +55,7 @@ llano_de_chajnantor: site_documentation: https://www.almaobservatory.org/en/home/ mauna_kea: - site_description: Mauna Kea Observatory (JCMT) + description: Mauna Kea Observatory (JCMT) region: mauna_kea latitude: 19.8228 longitude: -155.4770 @@ -63,7 +63,7 @@ mauna_kea: site_documentation: https://www.maunakeaobservatories.org/ meerkat: - site_description: Meerkat National Park (HERA, MeerKAT) + description: Meerkat National Park (HERA, MeerKAT) region: meerkat latitude: -30.7132 longitude: 21.4431 @@ -71,7 +71,7 @@ meerkat: site_documentation: https://www.protectedplanet.net/555703705?locale=en pic_de_bure: - site_description: Pic de Bure (NOEMA) + description: Pic de Bure (NOEMA) region: pic_de_bure latitude: 44.6339 longitude: 5.9079 @@ -79,7 +79,7 @@ pic_de_bure: site_documentation: https://iram-institute.org/observatories/noema/ pico_veleta: - site_description: IRAM 30m telescope + description: IRAM 30m telescope region: pico_veleta latitude: 37.0662 longitude: -3.39272 @@ -87,7 +87,7 @@ pico_veleta: site_documentation: https://iram-institute.org/observatories/30-meter-telescope/ qitai: - site_description: Qitai Radio Telescope + description: Qitai Radio Telescope region: qitai latitude: 43.6011 longitude: 89.6825 @@ -95,7 +95,7 @@ qitai: site_documentation: https://en.wikipedia.org/wiki/Qitai_Radio_Telescope sierra_negra: - site_description: Sierra Negra (LMT) + description: Sierra Negra (LMT) region: sierra_negra latitude: 18.9858 longitude: -97.3147 @@ -103,7 +103,7 @@ sierra_negra: site_documentation: http://lmtgtm.org/ srt: - site_description: Sardinia Radio Telescope (MISTRAL) + description: Sardinia Radio Telescope (MISTRAL) region: san_basilio latitude: 39.4928 longitude: 9.2450 @@ -111,7 +111,7 @@ srt: site_documentation: https://en.wikipedia.org/wiki/Sardinia_Radio_Telescope summit_station: - site_description: Summit Station (Greenland Telescope) + description: Summit Station (Greenland Telescope) region: summit_camp latitude: 72.5789 longitude: -38.4596 @@ -119,7 +119,7 @@ summit_station: site_documentation: https://geo-summit.org/ teide: - site_description: Teide Observatory (QUIJOTE) + description: Teide Observatory (QUIJOTE) region: teide latitude: 28.3007 longitude: -16.5110 diff --git a/maria/tests/test_atmosphere.py b/maria/tests/test_atmosphere.py index f5725a0..92e18e0 100644 --- a/maria/tests/test_atmosphere.py +++ b/maria/tests/test_atmosphere.py @@ -4,7 +4,7 @@ import maria from maria import Simulation -from maria.atmosphere import Atmosphere +from maria.atmosphere import Atmosphere, Spectrum, Weather @pytest.mark.parametrize("region_name", maria.all_regions) @@ -13,10 +13,13 @@ def test_atmosphere(region_name): @pytest.mark.parametrize("region_name", ["chajnantor", "green_bank", "south_pole"]) -def test_atmosphere_from_cache(region_name): - atmosphere = Atmosphere( - region=region_name, spectrum_from_cache=True, weather_from_cache=True - ) +def test_spectrum_from_cache(region_name): + spectrum = Spectrum(region=region_name, refresh_cache=True) + + +@pytest.mark.parametrize("region_name", ["chajnantor", "green_bank", "south_pole"]) +def test_weather_from_cache(region_name): + weather = Weather(region=region_name, refresh_cache=True) def test_atmosphere_2d(): @@ -24,6 +27,6 @@ def test_atmosphere_2d(): instrument="MUSTANG-2", plan="daisy", site="green_bank", - atmosphere_model="2d", + atmosphere="2d", ) tod = sim.run() diff --git a/maria/tests/test_cmb.py b/maria/tests/test_cmb.py new file mode 100644 index 0000000..7a02a2f --- /dev/null +++ b/maria/tests/test_cmb.py @@ -0,0 +1,10 @@ +import time + +import pytest + +import maria +from maria.cmb import generate_cmb + + +def test_generate_cmb(): + generate_cmb(nside=1024) diff --git a/maria/tests/test_coordinates.py b/maria/tests/test_coordinates.py index 540e0bc..43b663d 100644 --- a/maria/tests/test_coordinates.py +++ b/maria/tests/test_coordinates.py @@ -1,7 +1,6 @@ import numpy as np import pytest - -from ..coords import dx_dy_to_phi_theta, phi_theta_to_dx_dy +from todder.coords import dx_dy_to_phi_theta, phi_theta_to_dx_dy def test_offsets_transform(): diff --git a/maria/tests/test_noise.py b/maria/tests/test_noise.py new file mode 100644 index 0000000..06bc3f6 --- /dev/null +++ b/maria/tests/test_noise.py @@ -0,0 +1,33 @@ +import numpy as np +import pytest + +from maria import Simulation + + +@pytest.mark.noise +def test_linear_angular_model(): + sim = Simulation( + instrument="MUSTANG-2", + plan="daisy", + site="green_bank", + noise=True, + ) + tod = sim.run() + + target_error = sim.instrument.dets.NEP / np.sqrt(sim.plan.duration) + + scaled_residuals = ( + tod.data.compute().mean(axis=1) / target_error + ) # this is should be distributed as a zero-mean unit-variance Gaussian + + min_noise = 0.9 + max_noise = 1.1 + if scaled_residuals.std() < min_noise: + raise RuntimeError( + f"Noise residuals are too low ({scaled_residuals.std():.03f} < {min_noise})" + ) + + if scaled_residuals.std() > max_noise: + raise RuntimeError( + f"Noise residuals are too high ({scaled_residuals.std():.03f} > {max_noise})" + ) diff --git a/maria/tests/test_noise_sims.py b/maria/tests/test_noise_sims.py deleted file mode 100644 index 88f3c34..0000000 --- a/maria/tests/test_noise_sims.py +++ /dev/null @@ -1,14 +0,0 @@ -import pytest - -from maria.noise import NoiseSimulation - - -@pytest.mark.noise -def test_linear_angular_model(): - noise_sim = NoiseSimulation( - instrument="MUSTANG-2", - plan="daisy", - site="green_bank", - white_noise_level=1e0, - ) - tod = noise_sim.run() diff --git a/maria/tests/test_pipeline.py b/maria/tests/test_pipeline.py index c9d3f1c..b7fe078 100644 --- a/maria/tests/test_pipeline.py +++ b/maria/tests/test_pipeline.py @@ -3,21 +3,23 @@ import numpy as np import pytest +import maria from maria import Simulation from maria.map.mappers import BinMapper -from ..utils.io import fetch_cache +from ..io import fetch_cache here, this_filename = os.path.split(__file__) -TEST_MAP_URL = ( - "https://github.com/thomaswmorris/maria-data/raw/master/maps/cluster.fits" +TEST_MAP_SOURCE_URL = ( + "https://github.com/thomaswmorris/maria-data/raw/master/maps/cluster.fits" # noqa ) +TEST_MAP_CACHE_PATH = "/tmp/maria-data/maps/test.fits" @pytest.mark.mock_obs def test_mustang2(): - fetch_cache(TEST_MAP_URL, "/tmp/test_map.fits", refresh=True) + fetch_cache(TEST_MAP_SOURCE_URL, TEST_MAP_CACHE_PATH) pointing_center = (73.5287496858916, 2.961663679507145) pixel_size = 8.71452898559111e-05 @@ -25,27 +27,11 @@ def test_mustang2(): sample_rate = 100 scan_velocity = 38 / 3600 - inputfile = "/tmp/test_map.fits" - outfile_map = "/tmp/test_map_output.fits" + instrument = maria.get_instrument("MUSTANG-2") + site = maria.get_site("green_bank") - atm_model = "2d" - - sim = Simulation( - # Mandatory minimal weither settings - # --------------------- - instrument="MUSTANG-2", # Instrument type - site="green_bank", # Site - plan="daisy", # Scanning strategy - atmosphere_model=atm_model, # atmospheric model - # True sky input - # --------------------- - map_file=inputfile, # Input files must be a fits file. - map_units="Jy/pixel", # Units of the input map in Kelvin Rayleigh Jeans (K, defeault) or Jy/pixel - map_res=pixel_size, # resolution of the map - map_center=pointing_center, # RA & Dec in degree - map_freqs=[93], - # MUSTANG-2 Observational setup - # ----------------------------s + plan = maria.get_plan( + "daisy", scan_options={ "radius": 4.0 / 60.0, # The radius of the Daisy scan in degrees "speed": scan_velocity, # scan velocity in when the scan goes through the center deg/s @@ -53,9 +39,23 @@ def test_mustang2(): duration=duration, # Seconds sample_rate=sample_rate, # Hz scan_center=pointing_center, # Degrees - pointing_frame="ra_dec", # Frame + frame="ra_dec", # Frame start_time="2022-02-11T23:00:00", # observation date - pwv_rms_frac=0.005, # level of atmospheric fluctuations + ) + + map = maria.map.from_fits( + filename=TEST_MAP_CACHE_PATH, + resolution=pixel_size, + frequency=90, + center=pointing_center, + ) + + sim = Simulation( + instrument=instrument, # Instrument type + site=site, # Site + plan=plan, # Scanning strategy + atmosphere="2d", # atmospheric model + map=map, ) tod = sim.run() @@ -81,4 +81,4 @@ def test_mustang2(): mapper.run() - mapper.save_maps(outfile_map) + mapper.save_maps("/tmp/test-output.fits") diff --git a/maria/tests/test_tod.py b/maria/tests/test_tod.py index 868ba1a..8b1fd51 100644 --- a/maria/tests/test_tod.py +++ b/maria/tests/test_tod.py @@ -1,17 +1,17 @@ -import os +# import os -import pytest +# import pytest -from maria.tod.mustang2 import load_mustang2_tod +# from maria.tod.mustang2 import load_mustang2_tod -from ..utils.io import fetch_cache +# from ..io import fetch_cache -here, this_filename = os.path.split(__file__) +# here, this_filename = os.path.split(__file__) -TEST_MUSTANG_TOD_URL = "https://github.com/thomaswmorris/maria-data/raw/master/tods/mustang2/m2_sample_tod_60s.fits" +# TEST_MUSTANG_TOD_URL = "https://github.com/thomaswmorris/maria-data/raw/master/tods/mustang2/m2_sample_tod_60s.fits" -def test_tod_io_fits_mustang2(): - fetch_cache(TEST_MUSTANG_TOD_URL, "/tmp/test_m2_tod.fits", refresh=True) - tod = load_mustang2_tod(fname="/tmp/test_m2_tod.fits") - tod.to_fits(fname="/tmp/test_save_m2_tod.fits") +# def test_tod_io_fits_mustang2(): +# fetch_cache(TEST_MUSTANG_TOD_URL, "/tmp/test_m2_tod.fits", refresh=True) +# tod = load_mustang2_tod(fname="/tmp/test_m2_tod.fits") +# tod.to_fits(fname="/tmp/test_save_m2_tod.fits") diff --git a/maria/tests/test_weather.py b/maria/tests/test_weather.py deleted file mode 100644 index 9bd4ada..0000000 --- a/maria/tests/test_weather.py +++ /dev/null @@ -1,15 +0,0 @@ -import pytest - -import maria -from maria.atmosphere import Weather - - -@pytest.mark.parametrize("region_name", maria.all_regions) -def test_weather(region_name): - weather = Weather(region=region_name) - print(f"pwv={weather.pwv:.03f}mm") - - -@pytest.mark.parametrize("region_name", ["chajnantor", "green_bank", "south_pole"]) -def test_weather_from_cache(region_name): - weather = Weather(region=region_name, from_cache=True) diff --git a/maria/tod/__init__.py b/maria/tod/__init__.py deleted file mode 100644 index 9232c87..0000000 --- a/maria/tod/__init__.py +++ /dev/null @@ -1,268 +0,0 @@ -import functools -import json - -import h5py -import numpy as np -import pandas as pd -import scipy as sp -from astropy.io import fits - -from ..coords import Coordinates - - -class TOD: - """ - Time-ordered data. This has per-detector pointing and data. - """ - - def __init__( - self, - coords: Coordinates, - components: dict = {}, - units: dict = {}, - dets: pd.DataFrame = None, - abscal: float = 1.0, - ): - self.coords = coords - self.dets = dets - self.components = components - self.header = fits.header.Header() - self.abscal = abscal - self.units = units - - @functools.cached_property - def data(self): - return sum(self.components.values()) - - @functools.cached_property - def boresight(self): - return self.coords.boresight - - @property - def dt(self): - return np.diff(self.time).mean() - - @property - def fs(self): - return 1 / self.dt - - @property - def nd(self): - return self.data.shape[0] - - @property - def nt(self): - return self.data.shape[-1] - - def subset(self, det_mask=None, time_mask=None, band: str = None): - if band is not None: - 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) - - if time_mask is not None: - if len(time_mask) != self.nt: - raise ValueError("The detector mask must have shape (n_dets,).") - - subset_coords = Coordinates( - time=self.time[time_mask], - phi=self.coords.az[:, time_mask], - theta=self.coords.el[:, time_mask], - location=self.location, - frame="az_el", - ) - - return TOD( - data=self.data[:, time_mask], - coords=subset_coords, - dets=self.dets, - units=self.units, - ) - - if det_mask is not None: - if not (len(det_mask) == self.nd): - raise ValueError("The detector mask must have shape (n_dets,).") - - subset_dets = self.dets.loc[det_mask] if self.dets is not None else None - - subset_coords = Coordinates( - time=self.time, - phi=self.coords.az[det_mask], - theta=self.coords.el[det_mask], - location=self.location, - frame="az_el", - ) - - return TOD( - data=self.data[det_mask], - coords=subset_coords, - dets=subset_dets, - units=self.units, - ) - - @property - def time(self): - return self.coords.time - - @property - def location(self): - return self.coords.location - - @property - def lat(self): - return np.round(self.location.lat.deg, 6) - - @property - def lon(self): - return np.round(self.location.lon.deg, 6) - - @property - def alt(self): - return np.round(self.location.height.value, 6) - - @property - def az(self): - return self.coords.az - - @property - def el(self): - return self.coords.el - - @property - def ra(self): - return self.coords.ra - - @property - def dec(self): - return self.coords.dec - - @property - def azim_phase(self): - return np.pi * ( - sp.signal.sawtooth(2 * np.pi * self.time / self.azim_scan_period, width=1) - + 1 - ) - - @property - def turnarounds(self): - azim_grad = sp.ndimage.gaussian_filter(np.gradient(self.azim), sigma=16) - return np.where(np.sign(azim_grad[:-1]) != np.sign(azim_grad[1:]))[0] - - def splits(self, target_split_time: float = None): - if target_split_time is None: - return list(zip(self.turnarounds[:-1], self.turnarounds[1:])) - else: - fs = self.fs - splits_list = [] - for s, e in self.splits(target_split_time=None): - split_time = self.time[e] - self.time[s] # total time in the split - n_splits = int( - np.ceil(split_time / target_split_time) - ) # number of new splits - n_split_samples = int( - target_split_time * fs - ) # number of samples per new split - for split_start in np.linspace(s, e - n_split_samples, n_splits).astype( - int - ): - splits_list.append( - (split_start, np.minimum(split_start + n_split_samples, e)) - ) - return splits_list - - def to_fits(self, fname, format="MUSTANG-2"): - """ - Save the TOD to a fits file - """ - - if format.lower() == "mustang-2": - header = fits.header.Header() - - header["AZIM"] = (self.coords.center_az, "radians") - header["ELEV"] = (self.coords.center_el, "radians") - header["BMAJ"] = (8.0, "arcsec") - header["BMIN"] = (8.0, "arcsec") - header["BPA"] = (0.0, "degrees") - - header["SITELAT"] = (self.lat, "Site Latitude") - header["SITELONG"] = (self.lon, "Site Longitude") - header["SITEELEV"] = (self.alt, "Site elevation (meters)") - - col01 = fits.Column( - name="DX ", format="E", array=self.coords.ra.flatten(), unit="radians" - ) - col02 = fits.Column( - name="DY ", - format="E", - array=self.coords.dec.flatten(), - unit="radians", - ) - col03 = fits.Column( - name="FNU ", - format="E", - array=self.data.flatten(), - unit=self.units["data"], - ) - col04 = fits.Column(name="UFNU ", format="E") - col05 = fits.Column( - name="TIME ", - format="E", - array=( - (self.time - self.time[0]) * np.ones_like(self.coords.ra) - ).flatten(), - unit="s", - ) - col06 = fits.Column(name="COL ", format="I") - col07 = fits.Column(name="ROW ", format="I") - col08 = fits.Column( - name="PIXID", - format="I", - array=( - np.arange(len(self.coords.ra), dtype=np.int16).reshape(-1, 1) - * np.ones_like(self.coords.ra) - ).flatten(), - ) - col09 = fits.Column( - name="SCAN ", format="I", array=np.zeros_like(self.coords.ra).flatten() - ) - col10 = fits.Column(name="ELEV ", format="E") - - hdu = fits.BinTableHDU.from_columns( - [col01, col02, col03, col04, col05, col06, col07, col08, col09, col10], - header=header, - ) - - hdu.writeto(fname, overwrite=True) - - def to_hdf(self, fname): - with h5py.File(fname, "w") as f: - f.createcomponentsset(fname) - - -class KeyNotFoundError(Exception): - def __init__(self, invalid_keys): - super().__init__(f"The key '{invalid_keys}' is not in the database.") - - -def check_nested_keys(keys_found, data, keys): - for key in data.keys(): - for i in range(len(keys)): - if keys[i] in data[key].keys(): - keys_found[i] = True - - -def check_json_file_for_key(keys_found, file_path, *keys_to_check): - with open(file_path, "r") as json_file: - data = json.load(json_file) - return check_nested_keys(keys_found, data, keys_to_check) - - -def test_multiple_json_files(files_to_test, *keys_to_find): - keys_found = np.zeros(len(keys_to_find)).astype(bool) - - for file_path in files_to_test: - check_json_file_for_key(keys_found, file_path, *keys_to_find) - - if np.sum(keys_found) != len(keys_found): - raise KeyNotFoundError(np.array(keys_to_find)[~keys_found]) diff --git a/maria/tod/atlast.py b/maria/tod/atlast.py deleted file mode 100644 index 01fd778..0000000 --- a/maria/tod/atlast.py +++ /dev/null @@ -1,48 +0,0 @@ -import numpy as np -from astropy.io import fits - -from .. import site -from ..coords import Coordinates -from ..instrument.array import Detectors -from . import TOD - - -def load_atlast_tod( - fname: str, hdu: int = 1, band_center: int = 93, band_width: int = 52 -): - f = fits.open(fname) - raw = f[hdu].data - - det_uids, det_counts = np.unique(raw["PIXID"], return_counts=True) - - if det_counts.std() > 0: - raise ValueError("Cannot reshape a ragged TOD.") - - n_dets = len(det_uids) - n_samp = det_counts.max() - - data = {"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)) - t = 1.6e9 + raw["time"].astype(float).reshape((n_dets, n_samp)).mean(axis=0) - - coords = Coordinates( - time=t, - phi=ra, - theta=dec, - location=site.get_location("llano_de_chajnantor"), - frame="ra_dec", - ) - - dets = Detectors.generate( - bands_config={ - "f093": { - "n_dets": n_dets, - "band_center": band_center, - "band_width": band_width, - } - } - ) - - return TOD(coords=coords, dets=dets, data=data, units={"data": "K"}) diff --git a/maria/tod/mustang2.py b/maria/tod/mustang2.py deleted file mode 100644 index b3903c0..0000000 --- a/maria/tod/mustang2.py +++ /dev/null @@ -1,40 +0,0 @@ -import numpy as np -from astropy.io import fits - -from .. import instrument, site -from ..coords import Coordinates -from . import TOD - - -def load_mustang2_tod(fname: str, hdu: int = 1): - f = fits.open(fname) - raw = f[hdu].data - - det_uids, det_counts = np.unique(raw["PIXID"], return_counts=True) - - if det_counts.std() > 0: - raise ValueError("Cannot reshape a ragged TOD.") - - n_dets = len(det_uids) - n_samp = det_counts.max() - - components = {"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)) - t = 1.6e9 + raw["time"].astype(float).reshape((n_dets, n_samp)).mean(axis=0) - - coords = Coordinates( - time=t, - phi=ra, - theta=dec, - location=site.get_location("green_bank"), - frame="ra_dec", - ) - - m2_config = instrument.get_instrument_config(instrument_name="MUSTANG-2") - m2_config["array"]["n"] = n_dets - - m2 = instrument.get_instrument(**m2_config) - - return TOD(coords=coords, dets=m2.dets, components=components, units={"data": "K"}) diff --git a/maria/utils/units.py b/maria/units.py similarity index 100% rename from maria/utils/units.py rename to maria/units.py diff --git a/maria/utils/__init__.py b/maria/utils/__init__.py index e0ae9b4..f79677d 100644 --- a/maria/utils/__init__.py +++ b/maria/utils/__init__.py @@ -1,15 +1,8 @@ # this is the junk drawer of functions from datetime import datetime -import numpy as np import pytz -from . import functions # noqa F401 -from . import io # noqa F401 -from . import linalg # noqa F401 -from . import signal # noqa F401 -from . import units # noqa F401 - def get_utc_day_hour(t): dt = datetime.fromtimestamp(t, tz=pytz.utc).replace(tzinfo=pytz.utc) @@ -23,158 +16,3 @@ def get_utc_year_day(t): def get_utc_year(t): return datetime.fromtimestamp(t, tz=pytz.utc).replace(tzinfo=pytz.utc).year - - -def get_daisy_offsets(phase, k=2.11): - r = np.sin(k * phase) - return r * np.cos(phase), r * np.sin(phase) - - -def daisy_pattern_miss_center( - phase, petals=3, radius=1, offset_factor=0.15, k1=1 / 3 * np.pi, k2=2 / 3 * np.pi -): - shifted_phase = phase + np.pi / 2 - z1 = np.sin(phase * k2) * np.exp(1j * phase / petals) - z2 = ( - offset_factor * np.sin(shifted_phase * k1) * np.exp(1j * shifted_phase / petals) - ) - z = z1 - z2 - - max_dist = np.atleast_1d(np.abs(z)).max() - z *= radius / max_dist - - -# daisy_sky_x, daisy_sky_y = get_pointing_offset(time, period, plan_type="daisy") - -# return ( -# centers[0] + throws[0] * r * np.cos(p), -# centers[1] + throws[1] * r * np.sin(p), -# ) - -# if plan_type == "box": -# return ( -# centers[0] -# + throws[0] * np.interp(p % (2 * np.pi), np.linspace(0, 2 * np.pi, 5), [-1, -1, +1, +1, -1]), -# centers[1] + throws[1] * np.interp(p, np.linspace(0, 2 * np.pi, 5), [-1, +1, +1, -1, -1]), -# ) - - -# def get_daisy_offsets(phase, k=2.11): -# r = np.sin(k * phase) -# return r * np.cos(phase), r * np.sin(phase) - - -# def get_pointing(time, **kwargs): -# """ -# Returns azimuth and elevation -# """ -# scan_center = kwargs.get("scan_center", (np.radians(10), np.radians(4.5))) -# scan_pattern = kwargs.get("scan_pattern", "stare") -# scan_speed = kwargs.get("scan_speed", 60) -# scan_radius = kwargs.get("scan_radius", np.radians(2)) - -# phase = 2 * np.pi * time / scan_speed - -# if scan_pattern == "daisy": -# dpox, dpoy = get_daisy_offsets(phase) -# return dx_dy_to_phi_theta( # noqa F401 -# scan_radius * dpox, scan_radius * dpoy, *scan_center -# ) - - -# # COORDINATE TRANSFORM UTILS -# class Planner: -# def __init__(self, Array): -# self.array = Array - -# def make_plans(self, start, end, ra, dec, chunk_time, static_config): -# start_time = datetime.fromisoformat(start).replace(tzinfo=pytz.utc).timestamp() -# end_time = datetime.fromisoformat(end).replace(tzinfo=pytz.utc).timestamp() - -# _unix = np.arange(start_time, end_time, chunk_time) -# _ra = np.radians(np.linspace(ra, ra, len(_unix))) -# _dec = np.radians(np.linspace(dec, dec, len(_unix))) - -# _az, _el = self.array.coordinator.transform( -# _unix, _ra, _dec, in_frame="ra_dec", out_frame="az_el" -# ) - -# min_el = np.degrees(np.minimum(_el[1:], _el[:-1])) - -# ok = (min_el > self.array.el_bounds[0]) & (min_el < self.array.el_bounds[1]) - -# self.time, self.az, self.el = _unix[1:][ok], _az[1:][ok], _el[1:][ok] - -# for start_time in _unix[:-1][ok]: -# yield dict( -# { -# "start_time": start_time, -# "end_time": start_time + chunk_time, -# "coord_center": (ra, dec), -# "coord_throw": (2, 2), -# "coord_frame": "ra_dec", -# }, -# **static_config, -# ) - - -# ================ ARRAY ================ - -# ================ STATS ================ - - -# def get_minimal_bounding_rotation_angle(z): # minimal-area rotation angle - -# H = sp.spatial.ConvexHull(points=np.vstack([np.real(z).ravel(), np.imag(z).ravel()]).T) -# HZ = z.ravel()[H.vertices] - -# HE = np.imag(HZ).max() - np.imag(HZ).min() -# HO = 0 -# # for z1,z2 in zip(HZ,np.roll(HZ,1)): -# for RHO in np.linspace(0, np.pi, 1024 + 1)[1:]: - -# # RHO = np.angle(z2-z1) -# RHZ = HZ * np.exp(1j * RHO) - -# im_width = np.imag(RHZ).max() - np.imag(RHZ).min() -# re_width = np.real(RHZ).max() - np.real(RHZ).min() - -# RHE = im_width # * re_width - -# if RHE < HE and re_width > im_width: -# HE = RHE -# HO = RHO - -# return HO - - -# def smallest_max_error_sample(items, max_error=1e0): - -# k = 1 -# cluster_mids = np.sort(sp.cluster.vq.kmeans(items, k_or_guess=1)[0]) -# while (np.abs(np.subtract.outer(items, cluster_mids)) / cluster_mids[None, :]).min(axis=1).max() > max_error: -# cluster_mids = np.sort(sp.cluster.vq.kmeans(items, k_or_guess=k)[0]) -# k += 1 - -# which_cluster = np.abs(np.subtract.outer(items, cluster_mids)).argmin(axis=1) - -# return cluster_mids, which_cluster - - -# def make_beam_filter(side, window_func, args): - -# beam_X, beam_Y = np.meshgrid(side, side) -# beam_R = np.sqrt(np.square(beam_X) + np.square(beam_Y)) -# beam_W = window_func(beam_R, *args) - -# return beam_W / beam_W.sum() - - -# def get_brightness_temperature(f_pb, pb, f_spec, spec): - -# return sp.integrate.trapz( -# sp.interpolate.interp1d(f_spec, spec, axis=-1)(f_pb) * pb, f_pb, axis=-1 -# ) / sp.integrate.trapz(pb, f_pb) - - -# ================ POINTING ================ diff --git a/maria/utils/weather.py b/maria/utils/weather.py deleted file mode 100644 index 34cf011..0000000 --- a/maria/utils/weather.py +++ /dev/null @@ -1,34 +0,0 @@ -import numpy as np - - -def get_vapor_pressure(air_temp, rel_hum): # units are (°K, %) - T = air_temp - 273.15 # in °C - a, b, c = 611.21, 17.67, 238.88 # units are Pa, ., °C - gamma = np.log(1e-2 * rel_hum) + b * T / (c + T) - return a * np.exp(gamma) - - -def get_saturation_pressure(air_temp): # units are (°K, %) - T = air_temp - 273.15 # in °C - a, b, c = 611.21, 17.67, 238.88 # units are Pa, ., °C - return a * np.exp(b * T / (c + T)) - - -def get_dew_point(air_temp, rel_hum): # units are (°K, %) - a, b, c = 611.21, 17.67, 238.88 # units are Pa, ., °C - p_vap = get_vapor_pressure(air_temp, rel_hum) - return c * np.log(p_vap / a) / (b - np.log(p_vap / a)) + 273.15 - - -def get_relative_humidity(air_temp, dew_point): - T, DP = air_temp - 273.15, dew_point - 273.15 # in °C - b, c = 17.67, 238.88 - return 1e2 * np.exp(b * DP / (c + DP) - b * T / (c + T)) - - -def relative_to_absolute_humidity(air_temp, rel_hum): - return 1e-2 * rel_hum * get_saturation_pressure(air_temp) / (461.5 * air_temp) - - -def absolute_to_relative_humidity(air_temp, abs_hum): - return 1e2 * 461.5 * air_temp * abs_hum / get_saturation_pressure(air_temp) diff --git a/pyproject.toml b/pyproject.toml index 1b3c045..336ef3b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -20,15 +20,16 @@ dependencies = [ "requests", "scipy", "tables", + "todder", "tqdm", ] requires-python = ">=3.9" authors = [ - { name = "Thomas Morris", email = "thomasmorris@princeton.edu" }, + { name = "Thomas Morris", email = "thomas.w.morris@yale.edu" }, { name = "Joshiwa van Marrevijk", email = "joshiwavanmarrewijk@eso.org" } ] maintainers = [ - { name = "Thomas Morris", email = "thomasmorris@princeton.edu" }, + { name = "Thomas Morris", email = "thomas.w.morris@yale.edu" }, { name = "Joshiwa van Marrevijk", email = "joshiwavanmarrewijk@eso.org" } ] license = {file = "LICENSE"} diff --git a/pytest.ini b/pytest.ini index fcc35a9..8e35f5d 100644 --- a/pytest.ini +++ b/pytest.ini @@ -1,6 +1,3 @@ [pytest] markers = - atmosphere: Test only the atmospheric simulation code - noise: Test only the noise simulation code - mock_obs: Test mock observation - sim: Test simulations + internet: These tests need an internet connection. From 02ed77ffd883fd65d3c61d9fbbf51dce9df9a36e Mon Sep 17 00:00:00 2001 From: Thomas Morris Date: Tue, 21 May 2024 00:07:36 +0200 Subject: [PATCH 2/2] all the tests pass --- maria/__init__.py | 14 +++++++------- maria/instrument/bands/__init__.py | 16 ++++++++++------ maria/sim/__init__.py | 23 ++++++----------------- maria/tests/test_atmosphere.py | 4 ++-- maria/utils/__init__.py | 2 ++ 5 files changed, 27 insertions(+), 32 deletions(-) diff --git a/maria/__init__.py b/maria/__init__.py index f0e4d80..0b015a1 100644 --- a/maria/__init__.py +++ b/maria/__init__.py @@ -2,10 +2,10 @@ # Benedicta tu in mulieribus, et benedictus fructus ventris tui, Iesus. # Sancta Maria, Mater Dei, ora pro nobis peccatoribus, nunc, et in hora mortis nostrae. -from . import utils # noqa F401 -from ._version import __version__, __version_tuple__ # noqa: F401 -from .instrument import all_instruments, get_instrument # noqa F401 -from .map import Map, mappers # noqa F401 -from .plan import all_plans, get_plan # noqa F401 -from .sim import Simulation # noqa F401 -from .site import all_regions, all_sites, get_site # noqa F401 +from . import utils # noqa +from ._version import __version__, __version_tuple__ # noqa +from .instrument import all_instruments, get_instrument # noqa +from .map import Map, mappers # noqa +from .plan import all_plans, get_plan # noqa +from .sim import Simulation # noqa +from .site import all_regions, all_sites, get_site # noqa diff --git a/maria/instrument/bands/__init__.py b/maria/instrument/bands/__init__.py index 0a6bd16..606cae6 100644 --- a/maria/instrument/bands/__init__.py +++ b/maria/instrument/bands/__init__.py @@ -99,14 +99,18 @@ def __init__( self.spectrum = Spectrum(region="chajnantor") - if NEP is not None: - if sensitivity is not None: - raise RuntimeError( - "When defining a band, you must specify exactly one of 'NEP' or 'sensitivity'." - ) + if (NEP is None) and (sensitivity is None): + self.NEP = 0 + + elif (NEP is not None) and (sensitivity is not None): + raise RuntimeError( + "When defining a band, you must specify exactly one of 'NEP' or 'sensitivity'." + ) # noqa + + elif NEP is not None: self.NEP = NEP - if sensitivity is not None: + elif sensitivity is not None: self.set_sensitivity( sensitivity, kind=sensitivity_kind ) # this sets the NEP automatically diff --git a/maria/sim/__init__.py b/maria/sim/__init__.py index 988851c..ef82368 100644 --- a/maria/sim/__init__.py +++ b/maria/sim/__init__.py @@ -28,12 +28,13 @@ def __init__( instrument: Instrument or str = "default", plan: Plan or str = "daisy", site: Site or str = "hoagie_haven", - map: Map = None, atmosphere: str = None, cmb: CMB or str = None, + map: Map = None, noise: bool = True, atmosphere_kwargs: dict = {}, cmb_kwargs: dict = {}, + map_kwargs: dict = {}, noise_kwargs: dict = {}, verbose: bool = True, ): @@ -60,22 +61,10 @@ def __init__( self.noise = noise - # self.noise = True - - # self.params = {} - - # for sub_type, sub_master_params in master_params.items(): - # self.params[sub_type] = {} - # if sub_type in ["instrument", "site", "plan"]: - # sub_type_dataclass = getattr(self, sub_type) - # for k in sub_type_dataclass.__dataclass_fields__.keys(): - # v = getattr(sub_type_dataclass, k) - # setattr(self, k, v) - # self.params[sub_type][k] = v - # else: - # for k, v in sub_master_params.items(): - # setattr(self, k, kwargs.get(k, v)) - # self.params[sub_type][k] = v + self.atmosphere_kwargs = atmosphere_kwargs + self.cmb_kwargs = cmb_kwargs + self.map_kwargs = map_kwargs + self.noise_kwargs = noise_kwargs if map: self.map = map diff --git a/maria/tests/test_atmosphere.py b/maria/tests/test_atmosphere.py index 92e18e0..639a2e4 100644 --- a/maria/tests/test_atmosphere.py +++ b/maria/tests/test_atmosphere.py @@ -12,12 +12,12 @@ def test_atmosphere(region_name): atmosphere = Atmosphere(region=region_name) -@pytest.mark.parametrize("region_name", ["chajnantor", "green_bank", "south_pole"]) +@pytest.mark.parametrize("region_name", ["chajnantor"]) def test_spectrum_from_cache(region_name): spectrum = Spectrum(region=region_name, refresh_cache=True) -@pytest.mark.parametrize("region_name", ["chajnantor", "green_bank", "south_pole"]) +@pytest.mark.parametrize("region_name", ["chajnantor"]) def test_weather_from_cache(region_name): weather = Weather(region=region_name, refresh_cache=True) diff --git a/maria/utils/__init__.py b/maria/utils/__init__.py index f79677d..73befb3 100644 --- a/maria/utils/__init__.py +++ b/maria/utils/__init__.py @@ -3,6 +3,8 @@ import pytz +from . import linalg, signal # noqa + def get_utc_day_hour(t): dt = datetime.fromtimestamp(t, tz=pytz.utc).replace(tzinfo=pytz.utc)