Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor to fix bugs and other problems #43

Merged
merged 5 commits into from
Dec 20, 2023
+1,902 −953
Merged

Refactor to fix bugs and other problems #43

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
06aa859
finished refactor
Dec 18, 2023
4fafa7c
refactored maps
Dec 18, 2023
11b5ce7
pointing plots
Dec 19, 2023
36ed52d
fix passband glitch
Dec 19, 2023
381d641
maps work yay
Dec 20, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
107 changes: 107 additions & 0 deletions Untitled.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

50 changes: 43 additions & 7 deletions docs/source/tutorials/Synthetic_AtLAST_observations.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -21,9 +21,21 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 2,
"metadata": {},
"outputs": [],
"outputs": [
{
"ename": "ImportError",
"evalue": "cannot import name 'mappers' from 'maria' (/Users/tom/repos/maria/maria/__init__.py)",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mImportError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[2], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mmaria\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m Simulation\n\u001b[0;32m----> 2\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mmaria\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m mappers\n",
"\u001b[0;31mImportError\u001b[0m: cannot import name 'mappers' from 'maria' (/Users/tom/repos/maria/maria/__init__.py)"
]
}
],
"source": [
"from maria import Simulation\n",
"from maria import mappers"
Expand All @@ -46,9 +58,21 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 3,
"metadata": {},
"outputs": [],
"outputs": [
{
"ename": "NameError",
"evalue": "name 'map_file' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[3], line 8\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[38;5;66;03m# - Input figure\u001b[39;00m\n\u001b[1;32m 7\u001b[0m cmap \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mRdBu_r\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m----> 8\u001b[0m inputfile \u001b[38;5;241m=\u001b[39m \u001b[43mmap_file\u001b[49m\n\u001b[1;32m 10\u001b[0m hdu \u001b[38;5;241m=\u001b[39m fits\u001b[38;5;241m.\u001b[39mopen(inputfile)\n\u001b[1;32m 11\u001b[0m hdu[\u001b[38;5;241m0\u001b[39m]\u001b[38;5;241m.\u001b[39mdata \u001b[38;5;241m=\u001b[39m hdu[\u001b[38;5;241m0\u001b[39m]\u001b[38;5;241m.\u001b[39mdata\n",
"\u001b[0;31mNameError\u001b[0m: name 'map_file' is not defined"
]
}
],
"source": [
"from astropy.io import fits\n",
"from astropy.wcs import WCS\n",
Expand Down Expand Up @@ -89,9 +113,21 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 4,
"metadata": {},
"outputs": [],
"outputs": [
{
"ename": "NameError",
"evalue": "name 'map_file' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"Cell \u001b[0;32mIn[4], line 11\u001b[0m\n\u001b[1;32m 1\u001b[0m sim \u001b[38;5;241m=\u001b[39m Simulation(\n\u001b[1;32m 2\u001b[0m \u001b[38;5;66;03m# Mandatory minimal weither settings\u001b[39;00m\n\u001b[1;32m 3\u001b[0m \u001b[38;5;66;03m# ---------------------\u001b[39;00m\n\u001b[1;32m 4\u001b[0m array\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mAtLAST\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;66;03m# Array type\u001b[39;00m\n\u001b[1;32m 5\u001b[0m pointing\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdaisy\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;66;03m# Scanning strategy\u001b[39;00m\n\u001b[1;32m 6\u001b[0m site\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mAPEX\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;66;03m# Site\u001b[39;00m\n\u001b[1;32m 7\u001b[0m \u001b[38;5;66;03m# atm_model = 'single_layer', # The atmospheric model, set to None if you want a noiseless observation.\u001b[39;00m\n\u001b[1;32m 8\u001b[0m atm_model\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m, \u001b[38;5;66;03m# The atmospheric model, set to None if you want a noiseless observation.\u001b[39;00m\n\u001b[1;32m 9\u001b[0m \u001b[38;5;66;03m# True sky input\u001b[39;00m\n\u001b[1;32m 10\u001b[0m \u001b[38;5;66;03m# ---------------------\u001b[39;00m\n\u001b[0;32m---> 11\u001b[0m map_file\u001b[38;5;241m=\u001b[39m\u001b[43mmap_file\u001b[49m, \u001b[38;5;66;03m# Input files must be a fits file.\u001b[39;00m\n\u001b[1;32m 12\u001b[0m \u001b[38;5;66;03m# map_file can also be set to None if are only interested in the noise\u001b[39;00m\n\u001b[1;32m 13\u001b[0m scan_center\u001b[38;5;241m=\u001b[39mmap_center,\n\u001b[1;32m 14\u001b[0m map_center\u001b[38;5;241m=\u001b[39mmap_center, \u001b[38;5;66;03m# RA & Dec in degree \u001b[39;00m\n\u001b[1;32m 15\u001b[0m integration_time \u001b[38;5;241m=\u001b[39m integration_time,\n\u001b[1;32m 16\u001b[0m \u001b[38;5;66;03m# Defeault Observational setup\u001b[39;00m\n\u001b[1;32m 17\u001b[0m \u001b[38;5;66;03m# ----------------------------\u001b[39;00m\n\u001b[1;32m 18\u001b[0m pointing_frame\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mra_dec\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;66;03m# frame\u001b[39;00m\n\u001b[1;32m 19\u001b[0m field_of_view\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m2.0\u001b[39m,\n\u001b[1;32m 20\u001b[0m dets\u001b[38;5;241m=\u001b[39m{\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mf090\u001b[39m\u001b[38;5;124m\"\u001b[39m: {\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mn\u001b[39m\u001b[38;5;124m\"\u001b[39m: \u001b[38;5;241m2000\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mband_center\u001b[39m\u001b[38;5;124m\"\u001b[39m: \u001b[38;5;241m92.0\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mband_width\u001b[39m\u001b[38;5;124m\"\u001b[39m: \u001b[38;5;241m52.0\u001b[39m}},\n\u001b[1;32m 21\u001b[0m start_time\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m2022-08-10T06:00:00\u001b[39m\u001b[38;5;124m\"\u001b[39m,\n\u001b[1;32m 22\u001b[0m scan_options\u001b[38;5;241m=\u001b[39m{\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mspeed\u001b[39m\u001b[38;5;124m\"\u001b[39m: \u001b[38;5;241m1\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mradius\u001b[39m\u001b[38;5;124m\"\u001b[39m: \u001b[38;5;241m1.3\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mpetals\u001b[39m\u001b[38;5;124m\"\u001b[39m: \u001b[38;5;241m2.11\u001b[39m},\n\u001b[1;32m 23\u001b[0m sample_rate\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m100.0\u001b[39m,\n\u001b[1;32m 24\u001b[0m \u001b[38;5;66;03m# Additional inputs:\u001b[39;00m\n\u001b[1;32m 25\u001b[0m \u001b[38;5;66;03m# ----------------------\u001b[39;00m\n\u001b[1;32m 26\u001b[0m weather_quantiles\u001b[38;5;241m=\u001b[39m{\n\u001b[1;32m 27\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mcolumn_water_vapor\u001b[39m\u001b[38;5;124m\"\u001b[39m: \u001b[38;5;241m0.1\u001b[39m\n\u001b[1;32m 28\u001b[0m }, \u001b[38;5;66;03m# Weather conditions specific for that site\u001b[39;00m\n\u001b[1;32m 29\u001b[0m map_units\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mJy/pixel\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;66;03m# Kelvin Rayleigh Jeans (K, defeault) or Jy/pixel\u001b[39;00m\n\u001b[1;32m 30\u001b[0m map_res\u001b[38;5;241m=\u001b[39mpixel_size, \u001b[38;5;66;03m# degree, overwrites header information\u001b[39;00m\n\u001b[1;32m 31\u001b[0m )\n",
"\u001b[0;31mNameError\u001b[0m: name 'map_file' is not defined"
]
}
],
"source": [
"sim = Simulation(\n",
" # Mandatory minimal weither settings\n",
Expand Down
352 changes: 290 additions & 62 deletions docs/source/tutorials/Synthetic_MUSTANG2_observations.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

30 changes: 17 additions & 13 deletions maria/array.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,14 +27,17 @@

ARRAY_CONFIGS = utils.io.read_yaml(f"{here}/configs/arrays.yml")

DISPLAY_COLUMNS = ["description", "field_of_view", "primary_size", "bands"]
DISPLAY_COLUMNS = ["array_description", "field_of_view", "primary_size", "bands"]
supported_arrays_table = pd.DataFrame(ARRAY_CONFIGS).T

band_lists = []
for array_name, config in ARRAY_CONFIGS.items():
band_lists.append(
"/".join(
[str(band_config["band_center"]) for band_config in config["dets"].values()]
[
str(band_config["band_center"])
for band_config in config["detector_config"].values()
]
)
)
supported_arrays_table.loc[:, "bands"] = band_lists
Expand Down Expand Up @@ -97,19 +100,20 @@ def get_array_config(array_name=None, **kwargs):
if array_name not in ARRAY_CONFIGS.keys():
raise InvalidArrayError(array_name)
array_config = ARRAY_CONFIGS[array_name].copy()
for k, v in kwargs.items():
if k in all_array_params.keys():
array_config[k] = v
for key, value in kwargs.items():
if key in all_array_params.keys():
array_config[key] = value
else:
raise ValueError(f"'{k}' is not a valid argument for an array!")
raise ValueError(f"'{key}' is not a valid argument for an array!")
return array_config


def get_array(array_name="default", **kwargs):
"""
Get an array from a pre-defined config.
"""
return Array.from_config(get_array_config(array_name, **kwargs))
array_config = get_array_config(array_name=array_name, **kwargs)
return Array.from_config(array_config)


REQUIRED_DET_CONFIG_KEYS = ["n", "band_center", "band_width"]
Expand Down Expand Up @@ -177,7 +181,7 @@ class Array:
An array.
"""

description: str = ""
array_description: str = ""
primary_size: float = 5 # in meters
field_of_view: float = 1 # in deg
geometry: str = "hex"
Expand All @@ -188,7 +192,7 @@ class Array:
max_el_acc: float = np.inf # in deg/s^2
az_bounds: Tuple[float, float] = (0, 360) # in degrees
el_bounds: Tuple[float, float] = (0, 90) # in degrees
documentation: str = ""
array_documentation: str = ""
dets: pd.DataFrame = None # dets, it's complicated

def __repr__(self):
Expand Down Expand Up @@ -229,19 +233,19 @@ def baselines(self):

@classmethod
def from_config(cls, config):
if isinstance(config["dets"], Mapping):
if isinstance(config["detector_config"], Mapping):
field_of_view = config.get("field_of_view", 1)
geometry = config.get("geometry", "hex")
baseline = config.get("baseline", 0) # default to zero baseline
dets = generate_dets_from_config(
config["dets"],
bands=config["detector_config"],
field_of_view=field_of_view,
geometry=geometry,
baseline=baseline,
)

return cls(
description=config["description"],
array_description=config["array_description"],
primary_size=config["primary_size"],
field_of_view=field_of_view,
baseline=baseline,
Expand All @@ -252,7 +256,7 @@ def from_config(cls, config):
max_el_acc=config["max_el_acc"],
az_bounds=tuple(config["az_bounds"]),
el_bounds=tuple(config["el_bounds"]),
documentation=config["documentation"],
array_documentation=config["array_documentation"],
dets=dets,
)

Expand Down
200 changes: 196 additions & 4 deletions maria/atmosphere/__init__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,15 +1,207 @@
from .linear_angular import LinearAngularSimulation # noqa F401
from .single_layer import SingleLayerSimulation # noqa F401
# from .linear_angular import LinearAngularSimulation # noqa F401
import os
import time as ttime

# how do we do the bands? this is a great question.
# because all practical telescope instrumentation assume a constant band
import h5py
import numpy as np
import scipy as sp

# from .kolmogorov_taylor import KolmogorovTaylorSimulation
from tqdm import tqdm

from .. import utils
from ..base import BaseSimulation
from .turbulent_layer import TurbulentLayer # noqa F401

# how do we do the bands? this is a great question.
# because all practical telescope instrumentation assume a constant band


ATMOSPHERE_PARAMS = {
"min_depth": 500,
"max_depth": 3000,
"n_layers": 4,
"min_beam_res": 4,
"turbulent_outer_scale": 500,
}


here, this_filename = os.path.split(__file__)


class AtmosphericSpectrum:
def __init__(self, filepath):
"""
A dataclass to hold spectra as attributes
"""
with h5py.File(filepath, "r") as f:
self.side_nu_GHz = f["side_nu_GHz"][:].astype(float)
self.side_elevation_deg = f["side_elevation_deg"][:].astype(float)
self.side_line_of_sight_pwv_mm = f["side_line_of_sight_pwv_mm"][:].astype(
float
)
self.temperature_rayleigh_jeans_K = f["temperature_rayleigh_jeans_K"][
:
].astype(float)
self.phase_delay_um = f["phase_delay_um"][:].astype(float)


class AtmosphereMixin:
def _initialize_atmosphere(self):
"""
This assume that BaseSimulation.__init__() has been called.
"""

print("Initializing atmosphere")

utils.validate_pointing(self.det_coords.az, self.det_coords.el)

# self.min_atmosphere_height = min_atmosphere_height
# self.max_atmosphere_height = max_atmosphere_height

self._spectrum_filepath = f"{here}/spectra/{self.site.region}.h5"
self.spectrum = (
AtmosphericSpectrum(filepath=self._spectrum_filepath)
if os.path.exists(self._spectrum_filepath)
else None
)

if self.atmosphere_model == "2d":
self.n_layers = 6

self.layer_depths = np.linspace(500, 3000, self.n_layers)
self.layers = []

for layer_index, layer_depth in enumerate(self.layer_depths):
layer = TurbulentLayer(
array=self.array,
boresight=self.boresight,
depth=layer_depth,
weather=self.weather,
)

self.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_layers, self.array.n_dets, self.pointing.n_time))

for layer_index, layer in enumerate(self.layers):
layer.generate()
layer_data[layer_index] = layer.sample()

return layer_data

def _simulate_2d_atmospheric_fluctuations(self):
"""
Simulate layers of two-dimensional turbulence.
"""

turbulence = self._simulate_2d_turbulence()

rel_layer_scaling = np.interp(
self.site.altitude + self.layer_depths[:, None, None] * np.sin(self.EL),
self.weather.altitude_levels,
self.weather.absolute_humidity,
)
rel_layer_scaling /= np.sqrt(np.square(rel_layer_scaling).sum(axis=0)[None])

self.layer_scaling = self.pwv_rms_frac * self.weather.pwv * rel_layer_scaling

self.line_of_sight_pwv = (
self.weather.pwv + (self.layer_scaling * turbulence).sum(axis=0)
) / np.sin(self.EL)

# layer_boundaries = np.linspace(self.min_layer_height, self.max_layer_height, n_layers + 1)
# self.layer_heights = 0.5 * (layer_boundaries[1:] + layer_boundaries[:-1])
# self.n_layers = n_layers

# self.compute_autoregression_boundaries_2d(layer_depth)
# self.min_beam_res = min_atmosphere_beam_res

return

# def _simulate_integrated_water_vapor(self):
# detector_values = self.simulate_normalized_effective_water_vapor()

# # this is "zenith-scaled"
# self.line_of_sight_pwv = (
# self.weather.pwv
# * (1.0 + self.pwv_rms_frac * detector_values)
# / np.sin(self.EL)
# )

@property
def EL(self):
return utils.coords.dx_dy_to_phi_theta(
self.array.offset_x[:, None],
self.array.offset_y[:, None],
self.boresight.az,
self.boresight.el,
)[1]

def _simulate_atmospheric_emission(self, units="K_RJ"):
start_time = ttime.monotonic()

if units == "K_RJ": # Kelvin Rayleigh-Jeans
self._simulate_atmospheric_fluctuations()
self.data["atmosphere"] = np.empty(
(self.array.n_dets, self.pointing.n_time), dtype=np.float32
)

for uband in tqdm(self.array.ubands, desc="Sampling atmosphere"):
# for uband in self.array.ubands:
band_mask = self.array.dets.band.values == uband

det_nu_samples = np.linspace(
self.array.band_min[band_mask], self.array.band_max[band_mask], 64
).mean(axis=-1)

det_temperature_grid = sp.interpolate.interp1d(
self.spectrum.side_nu_GHz,
self.spectrum.temperature_rayleigh_jeans_K,
axis=-1,
)(det_nu_samples).mean(axis=-1)

band_T_RJ_interpolator = sp.interpolate.RegularGridInterpolator(
(
self.spectrum.side_line_of_sight_pwv_mm,
self.spectrum.side_elevation_deg,
),
det_temperature_grid,
)

self.data["atmosphere"][band_mask] = band_T_RJ_interpolator(
(
self.line_of_sight_pwv[band_mask],
np.degrees(self.det_coords.el[band_mask]),
)
)

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

if self.verbose:
print(
f"ran atmospheric simulation in {ttime.monotonic() - start_time:.01f} seconds"
)


class AtmosphereSimulation(BaseSimulation, AtmosphereMixin):
def __init__(self, *args, **kwargs):
super(AtmosphereSimulation, self).__init__(*args, **kwargs)
self._initialize_atmosphere()
Loading
Loading