diff --git a/maria/atmosphere/__init__.py b/maria/atmosphere/__init__.py index bfdd1ab..081c67f 100644 --- a/maria/atmosphere/__init__.py +++ b/maria/atmosphere/__init__.py @@ -8,6 +8,7 @@ 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 @@ -92,7 +93,7 @@ def _initialize_atmosphere(self): This assume that BaseSimulation.__init__() has been called. """ - utils.validate_pointing(self.det_coords.az, self.det_coords.el) + validate_pointing(self.det_coords.az, self.det_coords.el) if self.atmosphere_model == "2d": self.turbulent_layer_depths = np.linspace( diff --git a/maria/atmosphere/turbulent_layer.py b/maria/atmosphere/turbulent_layer.py index f195800..c278d99 100644 --- a/maria/atmosphere/turbulent_layer.py +++ b/maria/atmosphere/turbulent_layer.py @@ -94,9 +94,9 @@ def __init__( ] # find the detector offsets which form a convex hull - self.detector_offsets = np.c_[ - self.instrument.offset_x, self.instrument.offset_y - ] + self.detector_offsets = np.radians( + np.c_[self.instrument.sky_x, self.instrument.sky_y] + ) # add a small circle of offsets to account for pesky zeros unit_circle_complex = np.exp(1j * np.linspace(0, 2 * np.pi, 64 + 1)[:-1]) diff --git a/maria/coords/coords.c b/maria/coords/coords.c new file mode 100644 index 0000000..fd54433 --- /dev/null +++ b/maria/coords/coords.c @@ -0,0 +1,15 @@ +// Simple C program to display "Hello World" + +// Header file for input output functions +#include + +// main function - +// where the execution of program begins +int main() +{ + + // prints hello world + printf("Hello World"); + + return 0; +} diff --git a/maria/instrument/__init__.py b/maria/instrument/__init__.py index 236d53c..83e8d79 100644 --- a/maria/instrument/__init__.py +++ b/maria/instrument/__init__.py @@ -12,7 +12,7 @@ from .. import utils from .bands import BandList # noqa F401 from .beams import compute_angular_fwhm, construct_beam_filter # noqa F401 -from .dets import Detectors # noqa F401 +from .detectors import Detectors # noqa F401 HEX_CODE_LIST = [ mpl.colors.to_hex(mpl.colormaps.get_cmap("Paired")(t)) @@ -63,11 +63,10 @@ class Instrument: An instrument. """ - description: str = "" + description: str = "An instrument." primary_size: float = None # in meters field_of_view: float = None # in deg baseline: float = None - bath_temp: float = None bands: BandList = None dets: pd.DataFrame = None # dets, it's complicated documentation: str = "" @@ -84,12 +83,8 @@ def from_config(cls, config): return cls(bands=dets.bands, dets=dets, **config) - # def __post_init__(self): - - # if self.baseline is None: - # unique_baselines = sorted(np.unique(self.dets.baseline)) - - # ... + def __post_init__(self): + self.field_of_view = np.round(self.dets.sky_x.ptp(), 3) def __repr__(self): nodef_f_vals = ( @@ -110,16 +105,16 @@ def ubands(self): return self.dets.bands.names @property - def offset_x(self): - return self.dets.offset_x + def sky_x(self): + return self.dets.sky_x @property - def offset_y(self): - return self.dets.offset_y + def sky_y(self): + return self.dets.sky_y @property def offsets(self): - return np.c_[self.offset_x, self.offset_y] + return np.c_[self.sky_x, self.sky_y] @property def baseline_x(self): diff --git a/maria/instrument/detectors.py b/maria/instrument/detectors.py new file mode 100644 index 0000000..11a1509 --- /dev/null +++ b/maria/instrument/detectors.py @@ -0,0 +1,388 @@ +import os +import warnings +from collections.abc import Mapping + +import matplotlib as mpl +import numpy as np +import pandas as pd + +from .. import utils +from ..constants import c +from .bands import BandList, all_bands + +here, this_filename = os.path.split(__file__) + +HEX_CODE_LIST = [ + mpl.colors.to_hex(mpl.colormaps.get_cmap("Paired")(t)) + for t in [*np.linspace(0.05, 0.95, 12)] +] + +DET_COLUMN_TYPES = { + "array": "str", + "uid": "str", + "band_name": "str", + "band_center": "float", + "sky_x": "float", + "sky_y": "float", + "baseline_x": "float", + "baseline_y": "float", + "baseline_z": "float", + "pol_angle": "float", + "pol_label": "str", + "primary_size": "float", + "bath_temp": "float", + "abs_cal_rj": "float", + "time_constant": "float", + "white_noise": "float", + "pink_noise": "float", + "efficiency": "float", + "abs_cal_rj": "float", +} + +SUPPORTED_ARRAY_PACKINGS = ["hex", "square", "sunflower"] +SUPPORTED_ARRAY_SHAPES = ["hex", "square", "circle"] + + +def generate_array_offsets( + n: int = None, + diameter: float = None, + spacing: float = None, + packing: str = "hex", + shape: str = "hex", +) -> np.array: + n_params = (n is not None) + (diameter is not None) + (spacing is not None) + if n_params < 2: + raise ValueError("You must specify at least two of (n, diameter, spacing)") + + if n_params == 3: + warnings.warn( + "'n', 'diameter', and 'spacing' were all supplied for array generation. Ignoring spacing parameter." + ) + + if packing not in SUPPORTED_ARRAY_PACKINGS: + raise ValueError(f"'packing' must be one of {SUPPORTED_ARRAY_PACKINGS}") + + if shape not in SUPPORTED_ARRAY_SHAPES: + raise ValueError(f"'shape' must be one of {SUPPORTED_ARRAY_SHAPES}") + + if n is None: + r = int(np.ceil(diameter / spacing / 2)) + + if shape == "hex": + # a packed hexagon of radius $r$ has $n = 3r^2 - 3r + 1$ points in it + n = int(np.ceil(3 * r**2 - 3 * r + 1)) + + if shape == "circle": + n = int(np.ceil(np.pi * r**2)) + + if shape == "square": + n = 4 * r**2 + + if spacing is None: + if shape == "hex": + # a packed hexagon of radius $r$ has $n = 3r^2 - 3r + 1$ points in it + r = int((np.sqrt(12 * n - 3) + 3) / 6) + spacing = diameter / (2 * r) + + if shape == "circle": + spacing = np.pi * diameter / (4 * np.sqrt(n)) + + if shape == "square": + spacing = diameter / np.sqrt(n) + + super_n = n * 2 + + if packing == "square": + s = int(np.ceil(np.sqrt(super_n))) + side = np.arange(-s, s + 1, dtype=float) + x, y = [foo.ravel() for foo in np.meshgrid(side, side)] + + if packing == "hex": + s = int(np.ceil((np.sqrt(12 * super_n - 3) + 3) / 6)) + side = np.arange(-s, s + 1, dtype=float) + x, y = np.meshgrid(side, side) + y[:, ::2] -= 0.5 + x *= np.sqrt(3) / 2 + x, y = x.ravel(), y.ravel() + + if packing == "sunflower": + i = np.arange(super_n) + golden_angle = np.pi * (3.0 - np.sqrt(5.0)) + x = 0.586 * np.sqrt(i) * np.cos(golden_angle * i) + y = 0.586 * np.sqrt(i) * np.sin(golden_angle * i) + + n_sides = {"square": 4, "hex": 6, "circle": 256}[shape] + + r = np.sqrt(x**2 + y**2) + p = np.arctan2(y, x) + ngon_distance = r * np.cos(np.arcsin(np.sin(n_sides / 2 * p)) * 2 / n_sides) + + subset_index = np.argsort(ngon_distance)[:n] + + return spacing * np.c_[x[subset_index], y[subset_index]] + + +def generate_dets( + bands: list, + n: int = None, + field_of_view: float = None, + array_spacing: float = None, + array_packing: tuple = "hex", + array_shape: tuple = "circle", + array_offset: tuple = (0.0, 0.0), + baseline_diameter: float = 0, + baseline_packing: str = "sunflower", + baseline_shape: str = "circle", + baseline_offset: tuple = (0.0, 0.0, 0.0), + polarized: bool = False, + bath_temp: float = 0, +): + dets = pd.DataFrame() + + detector_offsets = generate_array_offsets( + n=n, + diameter=field_of_view, + spacing=array_spacing, + packing=array_packing, + shape=array_shape, + ) + + baselines = generate_array_offsets( + n=len(detector_offsets), + diameter=baseline_diameter, + packing=baseline_packing, + shape=baseline_shape, + ) + + if polarized: + pol_angles = np.random.uniform(low=0, high=360, size=len(detector_offsets)) + pol_labels = np.r_[["A" for _ in pol_angles], ["B" for _ in pol_angles]] + pol_angles = np.r_[pol_angles, (pol_angles + 90) % 360] + detector_offsets = np.r_[detector_offsets, detector_offsets] + baselines = np.r_[baselines, baselines] + + else: + pol_angles = np.zeros(len(detector_offsets)) + pol_labels = ["A" for i in pol_angles] + + for band in bands: + band_dets = pd.DataFrame( + index=np.arange(len(detector_offsets)), + columns=["band_name", "sky_x", "sky_y"], + dtype=float, + ) + + band_dets.loc[:, "band_name"] = band + band_dets.loc[:, "sky_x"] = array_offset[0] + detector_offsets[:, 0] + band_dets.loc[:, "sky_y"] = array_offset[1] + detector_offsets[:, 1] + + band_dets.loc[:, "baseline_x"] = baseline_offset[0] + baselines[:, 0] + band_dets.loc[:, "baseline_y"] = baseline_offset[1] + baselines[:, 1] + band_dets.loc[:, "baseline_z"] = baseline_offset[2] + + band_dets.loc[:, "bath_temp"] = bath_temp + band_dets.loc[:, "pol_angle"] = pol_angles + band_dets.loc[:, "pol_label"] = pol_labels + + dets = pd.concat([dets, band_dets]) + + return dets + + +def validate_band_config(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 Detectors: + @classmethod + def from_config(cls, config): + """ + Instantiate detectors from a config. We pass the whole config and not just config["dets"] so + that the detectors can inherit instrument parameters if need be (e.g. the size of the primary aperture). + """ + + config["dets"] = utils.io.flatten_config(config["dets"]) + + bands_config = {} + + df = pd.DataFrame( + {col: pd.Series(dtype=dtype) for col, dtype in DET_COLUMN_TYPES.items()} + ) + + for array in config["dets"]: + array_dets_config = config["dets"][array] + + if "band" in array_dets_config: + array_dets_config["bands"] = [array_dets_config.pop("band")] + + if "file" in array_dets_config: + # if a file is supplied, assume it's a csv and read it in + array_df = pd.read_csv( + f'{here}/{array_dets_config["file"]}', index_col=0 + ) + + # if no bands were supplied, get them from the table and hope they're registered + if "bands" not in array_dets_config: + array_dets_config["bands"] = sorted(np.unique(array_df.band.values)) + + array_dets_config["bands"] = parse_bands_config( + array_dets_config["bands"] + ) + + else: + array_dets_config["bands"] = parse_bands_config( + array_dets_config["bands"] + ) + + if "beam_spacing" in array_dets_config: + if len(array_dets_config["bands"]) > 1: + raise ValueError( + "'beam_spacing' parameter is unhandled for an detector array with multiple bands." + ) + + nu = ( + 1e9 + * array_dets_config["bands"][ + list(array_dets_config["bands"].keys())[0] + ]["center"] + ) + beam_resolution_degrees = np.degrees( + 1.22 * (c / nu) / config["primary_size"] + ) + array_dets_config[ + "array_spacing" + ] = beam_resolution_degrees * array_dets_config.pop("beam_spacing") + + array_df = generate_dets(**array_dets_config) + + # add leading zeros to detector uids + fill_level = int(np.log(np.maximum(len(array_df) - 1, 1)) / np.log(10) + 1) + + uid_predix = f"{array}-" if array else "" + uids = [ + f"{uid_predix}{str(i).zfill(fill_level)}" for i in range(len(array_df)) + ] + + array_df.insert(0, "uid", uids) + array_df.insert(1, "array", array) + + df = pd.concat([df, array_df]) + + for band, band_config in array_dets_config["bands"].items(): + bands_config[band] = band_config + + df.index = np.arange(len(df)) + + for col in ["primary_size"]: + if df.loc[:, col].isna().any(): + df.loc[:, col] = config[col] + + for col in [ + "sky_x", + "sky_y", + "baseline_x", + "baseline_y", + "baseline_z", + "pol_angle", + "bath_temp", + ]: + if df.loc[:, col].isna().any(): + df.loc[:, col] = 0 + + bands = BandList.from_config(bands_config) + + for band in bands: + df.loc[df.band_name == band.name, "band_center"] = band.center + df.loc[df.band_name == band.name, "efficiency"] = band.efficiency + + return cls(df=df, bands=bands) + + def __repr__(self): + return self.df.__repr__() + + def _repr_html_(self): + return self.df._repr_html_() + + def __init__(self, df: pd.DataFrame, bands: dict = {}): + self.df = df + self.bands = bands + + for attr in [ + "time_constant", + "white_noise", + "pink_noise", + "efficiency", + "abs_cal_rj", + ]: + values = np.zeros(shape=self.n) + for band in self.bands: + values[self.band_name == band.name] = getattr(band, attr) + self.df.loc[:, attr] = values + + def __getattr__(self, attr): + if attr in self.df.columns: + return self.df.loc[:, attr].values.astype(DET_COLUMN_TYPES[attr]) + + raise AttributeError(f"'Detectors' object has no attribute '{attr}'") + + def subset(self, band_name=None): + bands = BandList([self.bands[band_name]]) + return Detectors(bands=bands, df=self.df.loc[self.band_name == band_name]) + + @property + def n(self): + return len(self.df) + + @property + def offset(self): + return np.c_[self.sky_x, self.sky_y] + + @property + def __len__(self): + return len(self.df) + + @property + def index(self): + return self.df.index.values + + @property + def ubands(self): + return list(self.bands.keys()) + + def passband(self, nu): + _nu = np.atleast_1d(nu) + + PB = np.zeros((len(self.df), len(_nu))) + + for band in self.bands: + PB[self.band_name == band.name] = band.passband(_nu) + + return PB diff --git a/maria/instrument/dets.py b/maria/instrument/dets.py deleted file mode 100644 index 066797e..0000000 --- a/maria/instrument/dets.py +++ /dev/null @@ -1,303 +0,0 @@ -import os -from collections.abc import Mapping - -import matplotlib as mpl -import numpy as np -import pandas as pd - -from .. import utils -from .bands import BandList, all_bands - -here, this_filename = os.path.split(__file__) - -HEX_CODE_LIST = [ - mpl.colors.to_hex(mpl.colormaps.get_cmap("Paired")(t)) - for t in [*np.linspace(0.05, 0.95, 12)] -] - -DET_COLUMN_TYPES = { - "tag": "str", - "uid": "str", - "band_name": "str", - "band_center": "float", - "offset_x": "float", - "offset_y": "float", - "baseline_x": "float", - "baseline_y": "float", - "baseline_z": "float", - "pol_angle": "float", - "primary_size": "float", - "bath_temp": "float", - "abs_cal_rj": "float", - "time_constant": "float", - "white_noise": "float", - "pink_noise": "float", - "efficiency": "float", - "abs_cal_rj": "float", -} - -SUPPORTED_GEOMETRIES = ["flower", "hex", "square", "circle"] - - -def hex_packed_circle_offsets(n): - """ - Returns an array of $n$ hexagonal offsets from the origin with diameter 1. - """ - - h = int(np.ceil((np.sqrt(12 * n - 3) + 3) / 6)) - - side = np.arange(-h, h + 1, dtype=float) - x, y = np.meshgrid(side, side) - - x[1::2] -= 0.5 - y *= np.sqrt(3) / 2 - - offsets = np.c_[x.ravel(), y.ravel()] - - distance_squared = np.sum(offsets**2, axis=1) - - subset_index = np.argsort(distance_squared)[:n] - - return 0.5 * offsets[subset_index] / np.sqrt(distance_squared[subset_index[-1]]) - - -def generate_offsets(n, shape="hex"): - if shape not in SUPPORTED_GEOMETRIES: - raise ValueError(f"'shape' must be one of {SUPPORTED_GEOMETRIES}.") - - if shape == "circle": - return hex_packed_circle_offsets(n).T - - if shape == "hex": - angles = np.linspace(0, 2 * np.pi, 6 + 1)[1:] + np.pi / 2 - z = np.array([0]) - layer = 0 - while len(z) < n: - for angle in angles: - for _z in layer * np.exp(1j * angle) + np.arange(layer) * np.exp( - 1j * (angle + 2 * np.pi / 3) - ): - z = np.r_[z, _z] - layer += 1 - z -= z.mean() - z *= 0.5 / np.abs(z).max() - - return np.c_[np.real(np.array(z[:n])), np.imag(np.array(z[:n]))].T - - if shape == "flower": - golden_ratio = np.pi * (3.0 - np.sqrt(5.0)) # golden angle in radians - z = np.zeros(n).astype(complex) - for i in range(n): - z[i] = np.sqrt((i / (n - 1)) * 2) * np.exp(1j * golden_ratio * i) - z *= 0.5 / np.abs(z.max()) - return np.c_[np.real(z), np.imag(z)].T - - if shape == "square": - side = np.linspace(-0.5, 0.5, int(np.ceil(np.sqrt(n)))) - DX, DY = np.meshgrid(side, side) - return np.c_[DX.ravel()[:n], DY.ravel()[:n]].T - - -def generate_dets( - n: int, - bands: list, - field_of_view: float = 0.0, - array_offset: tuple = (0.0, 0.0), - array_shape: tuple = "hex", - baseline_offset: tuple = (0.0, 0.0, 0.0), - baseline_diameter: float = 0.0, - baseline_shape: str = "flower", - **kwargs, -): - dets = pd.DataFrame() - - detector_offsets = field_of_view * generate_offsets(n=n, shape=array_shape) - - baselines = baseline_diameter * generate_offsets(n=n, shape=baseline_shape) - - for band in bands: - band_dets = pd.DataFrame( - index=np.arange(n), - columns=["band_name", "offset_x", "offset_y"], - dtype=float, - ) - - band_dets.loc[:, "band_name"] = band - band_dets.loc[:, "offset_x"] = np.radians(array_offset[0] + detector_offsets[0]) - band_dets.loc[:, "offset_y"] = np.radians(array_offset[1] + detector_offsets[1]) - - band_dets.loc[:, "baseline_x"] = baseline_offset[0] + baselines[0] - band_dets.loc[:, "baseline_y"] = baseline_offset[1] + baselines[1] - band_dets.loc[:, "baseline_z"] = baseline_offset[2] - - dets = pd.concat([dets, band_dets]) - - return dets - - -def validate_band_config(band): - if any([key not in band for key in ["center"]]): - raise ValueError("The band center 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, 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 - - if isinstance(bands, Mapping): - for name, band in bands.items(): - validate_band_config(band) - parsed_band_config[name] = band - - return parsed_band_config - - -class Detectors: - @classmethod - def from_config(cls, config): - """ - Instantiate detectors from a config. We pass the whole config and not just config["dets"] so - that the detectors can inherit instrument parameters if need be. - """ - - config["dets"] = utils.io.flatten_config(config["dets"]) - - bands_config = {} - - df = pd.DataFrame( - {col: pd.Series(dtype=dtype) for col, dtype in DET_COLUMN_TYPES.items()} - ) - - for tag in config["dets"]: - tag_dets_config = config["dets"][tag] - - if "file" in tag_dets_config: - # if a file is supplied, assume it's a csv and read it in - tag_df = pd.read_csv(f'{here}/{tag_dets_config["file"]}', index_col=0) - - # if no bands were supplied, get them from the table and hope they're registered - if "bands" not in tag_dets_config: - tag_dets_config["bands"] = sorted(np.unique(tag_df.band.values)) - - tag_dets_config["bands"] = parse_bands_config(tag_dets_config["bands"]) - - else: - tag_dets_config["bands"] = parse_bands_config(tag_dets_config["bands"]) - tag_df = generate_dets(**tag_dets_config) - - fill_level = int(np.log(len(tag_df) - 1) / np.log(10) + 1) - - uid_predix = f"{tag}-" if tag else "" - uids = [ - f"{uid_predix}{str(i).zfill(fill_level)}" for i in range(len(tag_df)) - ] - - tag_df.insert(0, "uid", uids) - tag_df.insert(1, "tag", tag) - - df = pd.concat([df, tag_df]) - - for band, band_config in tag_dets_config["bands"].items(): - bands_config[band] = band_config - - df.index = np.arange(len(df)) - - for col in ["primary_size"]: - if df.loc[:, col].isna().any(): - df.loc[:, col] = config[col] - - for col in [ - "offset_x", - "offset_y", - "baseline_x", - "baseline_y", - "baseline_z", - "pol_angle", - "bath_temp", - ]: - if df.loc[:, col].isna().any(): - df.loc[:, col] = 0 - - bands = BandList.from_config(bands_config) - - for band in bands: - df.loc[df.band_name == band.name, "band_center"] = band.center - df.loc[df.band_name == band.name, "efficiency"] = band.efficiency - - return cls(df=df, bands=bands) - - def __repr__(self): - return self.df.__repr__() - - def __repr_html__(self): - return self.df.__repr_html__() - - def __init__(self, df: pd.DataFrame, bands: dict = {}): - self.df = df - self.bands = bands - - for attr in [ - "time_constant", - "white_noise", - "pink_noise", - "efficiency", - "abs_cal_rj", - ]: - values = np.zeros(shape=self.n) - for band in self.bands: - values[self.band_name == band.name] = getattr(band, attr) - self.df.loc[:, attr] = values - - def __getattr__(self, attr): - if attr in self.df.columns: - return self.df.loc[:, attr].values.astype(DET_COLUMN_TYPES[attr]) - - raise AttributeError(f"'Detectors' object has no attribute '{attr}'") - - def subset(self, band_name=None): - bands = BandList([self.bands[band_name]]) - return Detectors(bands=bands, df=self.df.loc[self.band_name == band_name]) - - @property - def n(self): - return len(self.df) - - @property - def offset(self): - return np.c_[self.offset_x, self.offset_y] - - @property - def __len__(self): - return len(self.df) - - @property - def index(self): - return self.df.index.values - - @property - def ubands(self): - return list(self.bands.keys()) - - def passband(self, nu): - _nu = np.atleast_1d(nu) - - PB = np.zeros((len(self.df), len(_nu))) - - for band in self.bands: - PB[self.band_name == band.name] = band.passband(_nu) - - return PB diff --git a/maria/instrument/instruments.yml b/maria/instrument/instruments.yml index 0e54172..3b00f37 100644 --- a/maria/instrument/instruments.yml +++ b/maria/instrument/instruments.yml @@ -30,23 +30,26 @@ AdvACT: documentation: https://act.princeton.edu/overview/camera-specifications/advact dets: pa4: - n: 750 + n: 397 array_offset: [-0.8, -0.5] field_of_view: 1.0 array_shape: hex bands: [act/pa4/f150, act/pa4/f220] + polarized: True pa5: - n: 750 + n: 397 array_offset: [0.0, 1.0] field_of_view: 1.0 array_shape: hex bands: [act/pa5/f090, act/pa5/f150] + polarized: True pa6: - n: 750 + n: 397 array_offset: [0.8, -0.5] field_of_view: 1.0 array_shape: hex bands: [act/pa6/f090, act/pa6/f150] + polarized: True primary_size: 6 ALMA: @@ -67,7 +70,7 @@ AtLAST: bath_temp: 100.e-3 # in K bands: [atlast/f027, atlast/f039, atlast/f093, atlast/f150, atlast/f225, atlast/f280] field_of_view: 0.07 # in degrees - primary_size: 100 # in meters + primary_size: 50 # in meters documentation: https://greenbankobservatory.org/science/gbt-observers/mustang-2/ MUSTANG-2: diff --git a/maria/map/mappers.py b/maria/map/mappers.py index 103399e..494187c 100644 --- a/maria/map/mappers.py +++ b/maria/map/mappers.py @@ -148,7 +148,9 @@ def __init__( def run(self): self.band = sorted( - [band for tod in self.tods for band in np.unique(tod.dets.band_name)] + np.unique( + [band for tod in self.tods for band in np.unique(tod.dets.band_name)] + ) ) self.band_data = {} diff --git a/maria/plan/__init__.py b/maria/plan/__init__.py index c52b30f..8a75b30 100644 --- a/maria/plan/__init__.py +++ b/maria/plan/__init__.py @@ -12,8 +12,11 @@ from .. import coords, utils from . import patterns -MAX_VELOCITY_WARN = np.radians(10) # in deg/s -MAX_ACCELERATION_WARN = np.radians(10) # in deg/s +MAX_VELOCITY_WARN = 10 # in deg/s +MAX_ACCELERATION_WARN = 10 # in deg/s + +MIN_ELEVATION_WARN = 20 # in deg +MIN_ELEVATION_ERROR = 10 # in deg here, this_filename = os.path.split(__file__) @@ -47,6 +50,22 @@ def get_plan(scan_pattern="stare", **kwargs): return Plan(**plan_config) +class PointingError(Exception): + pass + + +def validate_pointing(azim, elev): + el_min = np.atleast_1d(elev).min() + if el_min < np.radians(MIN_ELEVATION_WARN): + warnings.warn( + f"Some detectors come within {MIN_ELEVATION_WARN} degrees of the horizon (el_min = {np.degrees(el_min):.01f}°)" + ) + if el_min <= np.radians(MIN_ELEVATION_ERROR): + raise PointingError( + f"Some detectors come within {MIN_ELEVATION_ERROR} degrees of the horizon (el_min = {np.degrees(el_min):.01f}°)" + ) + + @dataclass class Plan: """ diff --git a/maria/tests/test_instruments.py b/maria/tests/test_instruments.py index 984a85e..127cc3a 100644 --- a/maria/tests/test_instruments.py +++ b/maria/tests/test_instruments.py @@ -2,7 +2,7 @@ import maria -dets1 = { +subarray_dets = { "subarray-1": { "n": 500, "field_of_view": 2, @@ -17,14 +17,32 @@ }, } -dets2 = { +predefined_bands = { "n": 500, "field_of_view": 2, + "array_packing": "square", "array_shape": "hex", "bands": ["alma/f043", "alma/f078"], } -dets3 = {"file": "data/alma/alma.cycle1.total.csv"} +predefined_dets = {"file": "data/alma/alma.cycle1.total.csv"} + +beam_packing_dets = { + "subarray-1": { + "field_of_view": 0.2, + "beam_spacing": 1, + "array_shape": "hex", + "array_packing": "sunflower", + "band": {"center": 30, "width": 5}, + }, + "subarray-2": { + "n": 500, + "field_of_view": 0.2, + "array_shape": "circle", + "array_packing": "hex", + "bands": [{"center": 90, "width": 5}, {"center": 150, "width": 5}], + }, +} @pytest.mark.parametrize("instrument_name", maria.all_instruments) @@ -32,6 +50,8 @@ def test_get_instrument(instrument_name): instrument = maria.get_instrument(instrument_name) -@pytest.mark.parametrize("dets", [dets1, dets2, dets3]) +@pytest.mark.parametrize( + "dets", [subarray_dets, predefined_bands, predefined_dets, beam_packing_dets] +) def test_get_custom_array(dets): instrument = maria.get_instrument(dets=dets) diff --git a/maria/tests/test_pipeline.py b/maria/tests/test_pipeline.py index 8e4c9e1..c9d3f1c 100644 --- a/maria/tests/test_pipeline.py +++ b/maria/tests/test_pipeline.py @@ -18,7 +18,6 @@ @pytest.mark.mock_obs def test_mustang2(): fetch_cache(TEST_MAP_URL, "/tmp/test_map.fits", refresh=True) - map_size = 0.1 pointing_center = (73.5287496858916, 2.961663679507145) pixel_size = 8.71452898559111e-05 @@ -30,8 +29,6 @@ def test_mustang2(): outfile_map = "/tmp/test_map_output.fits" atm_model = "2d" - white_noise_level = 1.3e-2 - pink_noise_level = 2.4 sim = Simulation( # Mandatory minimal weither settings @@ -47,7 +44,6 @@ def test_mustang2(): map_res=pixel_size, # resolution of the map map_center=pointing_center, # RA & Dec in degree map_freqs=[93], - dets={"f093": {"n": 217, "bands": ["mustang2/f093"]}}, # MUSTANG-2 Observational setup # ----------------------------s scan_options={ diff --git a/maria/tod/__init__.py b/maria/tod/__init__.py index f4ee3b9..9260066 100644 --- a/maria/tod/__init__.py +++ b/maria/tod/__init__.py @@ -8,7 +8,7 @@ from .. import instrument, site from ..coords import Coordinates, get_center_phi_theta -from ..instrument.dets import Detectors +from ..instrument.detectors import Detectors class TOD: diff --git a/maria/utils/__init__.py b/maria/utils/__init__.py index d3d1ce9..e0ae9b4 100644 --- a/maria/utils/__init__.py +++ b/maria/utils/__init__.py @@ -1,5 +1,4 @@ # this is the junk drawer of functions -import warnings from datetime import datetime import numpy as np @@ -26,22 +25,6 @@ def get_utc_year(t): return datetime.fromtimestamp(t, tz=pytz.utc).replace(tzinfo=pytz.utc).year -class PointingError(Exception): - pass - - -def validate_pointing(azim, elev): - el_min = np.atleast_1d(elev).min() - if el_min < np.radians(10): - warnings.warn( - f"Some detectors come within 10 degrees of the horizon (el_min = {np.degrees(el_min):.01f}°)" - ) - if el_min <= 0: - raise PointingError( - f"Some detectors are pointing below the horizon (el_min = {np.degrees(el_min):.01f}°)" - ) - - def get_daisy_offsets(phase, k=2.11): r = np.sin(k * phase) return r * np.cos(phase), r * np.sin(phase) @@ -61,7 +44,7 @@ def daisy_pattern_miss_center( z *= radius / max_dist -# daisy_offset_x, daisy_offset_y = get_pointing_offset(time, period, plan_type="daisy") +# daisy_sky_x, daisy_sky_y = get_pointing_offset(time, period, plan_type="daisy") # return ( # centers[0] + throws[0] * r * np.cos(p), diff --git a/pyproject.toml b/pyproject.toml index 0d73809..1b3c045 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,6 +10,7 @@ readme = { file = "README.rst", content-type = "text/x-rst" } dependencies = [ "astropy", "cmocean", + "dask", "h5py", "healpy", "matplotlib",