Skip to content

Commit

Permalink
rename pointing to plan
Browse files Browse the repository at this point in the history
  • Loading branch information
Thomas Morris committed Feb 21, 2024
1 parent b333cba commit d9eb82c
Show file tree
Hide file tree
Showing 23 changed files with 203 additions and 164 deletions.
29 changes: 22 additions & 7 deletions docs/source/instruments.rst
Original file line number Diff line number Diff line change
Expand Up @@ -46,12 +46,27 @@ These are all valid:

.. code-block:: python
f090 = {"center": 90, "width": 30}
f150 = {"center": 150, "width": 30}
dets = {"n": 500,
"field_of_view": 2
dets = {
"subarray-1": {
"n": 500,
"field_of_view": 2,
"array_shape": "hex",
"bands": [f090, f150]}
my_custom_array = maria.get_instrument(dets=dets)
"bands": [{"center": 30, "width": 5}, {"center": 40, "width": 5}],
},
"subarray-2": {
"n": 500,
"field_of_view": 2,
"array_shape": "hex",
"bands": [{"center": 90, "width": 5}, {"center": 150, "width": 5}],
},
}
dets = {
"n": 500,
"field_of_view": 2,
"array_shape": "hex",
"bands": ["alma/f043", "alma/f078"],
}
dets = {"file": "path_to_some_dets_file.csv"}
4 changes: 3 additions & 1 deletion docs/source/papers.rst
Original file line number Diff line number Diff line change
@@ -1,4 +1,6 @@
Papers
======

`Morris et al. (2022) <https://arxiv.org/abs/2111.01319>`_
`T. W. Morris et al. (2022) <https://arxiv.org/abs/2111.01319>`_

`J. van Marrewijk, T. W. Morris et al. (2024) <https://arxiv.org/abs/2402.10731>`_
81 changes: 67 additions & 14 deletions docs/source/pointings.rst
Original file line number Diff line number Diff line change
@@ -1,19 +1,72 @@
Customizing pointings
+++++++++++++++++++++
===========
Pointings
===========

We can simulate observing (and later mapping) some celestial signal by supplying the simulation with a map to use as the ground truth.
The observing pointing is instantiated as an ``Pointing``.
The simplest way to get an pointing is to grab a pre-defined one, e.g.::

# The Atacama Cosmology Telescope
act = maria.get_pointing('ACT')


Simulations are made up of an array (defining the instrument), a site (where on earth it is), and a pointing (defining how it observes). Simulations for a few different telescopes might be instantiated as::
# The Atacama Large Millimeter Array
alma = maria.get_pointing('ALMA')

# MUSTANG-2
mustang2_sim = Simulation(array='MUSTANG-2',
site='green_bank',
pointing='daisy',
map_file=path_to_some_fits_file, # 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,
map_center=(10, 4.5), # RA & Dec. in degrees
map_freqs=[90],
degrees=True)
m2 = maria.get_pointing('MUSTANG-2')


+++++++++++++++++++++++
Customizing Pointings
+++++++++++++++++++++++

One way to customize pointings is to load a pre-defined pointing and pass extra parameters.
For example, we can give ACT twice the resolution with

.. code-block:: python
act = maria.get_pointing('ACT', primary_size=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,
"field_of_view": 2
"array_shape": "hex",
"bands": [f090, f150]}
my_custom_array = maria.get_pointing(dets=dets)
Actually, there are several valid ways to define an array of detectors.
These are all valid:

.. code-block:: python
dets = {
"subarray-1": {
"n": 500,
"field_of_view": 2,
"array_shape": "hex",
"bands": [{"center": 30, "width": 5}, {"center": 40, "width": 5}],
},
"subarray-2": {
"n": 500,
"field_of_view": 2,
"array_shape": "hex",
"bands": [{"center": 90, "width": 5}, {"center": 150, "width": 5}],
},
}
dets = {
"n": 500,
"field_of_view": 2,
"array_shape": "hex",
"bands": ["alma/f043", "alma/f078"],
}
dets = {"file": "path_to_some_dets_file.csv"}
2 changes: 1 addition & 1 deletion maria/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@
from . import utils # noqa F401
from .instrument import all_instruments, get_instrument # noqa F401
from .map import Map, mappers # noqa F401
from .pointing import all_pointings, get_pointing # 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 .tod import TOD # noqa F401
4 changes: 2 additions & 2 deletions maria/atmosphere/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -68,7 +68,7 @@ def _simulate_2d_turbulence(self):
"""

layer_data = np.zeros(
(self.n_atmosphere_layers, self.instrument.n_dets, self.pointing.n_time)
(self.n_atmosphere_layers, self.instrument.n_dets, self.plan.n_time)
)

layers = tqdm(self.turbulent_layers) if self.verbose else self.turbulent_layers
Expand Down Expand Up @@ -106,7 +106,7 @@ 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.pointing.n_time), dtype=np.float32
(self.instrument.n_dets, self.plan.n_time), dtype=np.float32
)

bands = (
Expand Down
2 changes: 1 addition & 1 deletion maria/atmosphere/turbulent_layer.py
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ def __init__(
self.res = res

# this is approximately correct
# self.depth = self.depth / np.sin(np.mean(self.pointing.el))
# self.depth = self.depth / np.sin(np.mean(self.plan.el))

# this might change
self.angular_outer_scale = turbulent_outer_scale / self.depth
Expand Down
10 changes: 7 additions & 3 deletions maria/cmb/__init__.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,7 @@
class CMB:
...


class CMBMixin:
...

Expand All @@ -6,8 +10,8 @@ class CMBMixin:
# """
# This simulates scanning over celestial sources.
# """
# def __init__(self, instrument, pointing, site, map_file, add_cmb=False, **kwargs):
# super().__init__(instrument, pointing, site)
# def __init__(self, instrument, plan, site, map_file, add_cmb=False, **kwargs):
# super().__init__(instrument, plan, site)

# pass

Expand Down Expand Up @@ -46,7 +50,7 @@ class CMBMixin:
# np.array([self.CMB_PS[bandnumber]]),
# [0],
# beam=None,
# seed=self.pointing.seed,
# seed=self.plan.seed,
# )[0]

# self.CMB_map += utils.Tcmb
Expand Down
38 changes: 16 additions & 22 deletions maria/map/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -164,7 +164,7 @@ def _initialize_map(self):
if not self.map_file:
return

self.input_map_file = self.map_file
self.map_file = self.map_file
hudl = ap.io.fits.open(self.map_file)

map_data = hudl[0].data
Expand Down Expand Up @@ -193,7 +193,7 @@ def _initialize_map(self):

self.map_data = map_data

self.input_map = Map(
self.map = Map(
data=map_data,
header=hudl[0].header,
freqs=np.atleast_1d(self.map_freqs),
Expand All @@ -206,30 +206,26 @@ def _initialize_map(self):
units=self.map_units,
)

self.input_map.header["HISTORY"] = "History_input_adjustments"
self.input_map.header["comment"] = "Changed input CDELT1 and CDELT2"
self.input_map.header["comment"] = (
"Changed surface brightness units to " + self.input_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.input_map.header["comment"] = "Repositioned the map on the sky"
self.map.header["comment"] = "Repositioned the map on the sky"

if self.input_map.inbright is not None:
self.input_map.data *= self.input_map.inbright / np.nanmax(
self.input_map.data
)
self.input_map.header["comment"] = "Amplitude is rescaled."
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.det_coords.offsets(
frame=self.map_frame, center=self.input_map.center
)
dx, dy = self.det_coords.offsets(frame=self.map_frame, center=self.map.center)

self.data["map"] = np.zeros((dx.shape))

freqs = tqdm(self.input_map.freqs) if self.verbose else self.input_map.freqs
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")
Expand All @@ -238,15 +234,13 @@ def _sample_maps(self):
nu_fwhm = beams.compute_angular_fwhm(
fwhm_0=self.instrument.primary_size, z=np.inf, f=1e9 * nu
)
nu_map_filter = beams.construct_beam_filter(
fwhm=nu_fwhm, res=self.input_map.res
)
nu_map_filter = beams.construct_beam_filter(fwhm=nu_fwhm, res=self.map.res)
filtered_nu_map_data = beams.separably_filter(
self.input_map.data[i], nu_map_filter
self.map.data[i], nu_map_filter
)

# band_res_radians = 1.22 * (299792458 / (1e9 * nu)) / self.instrument.primary_size
# band_res_pixels = band_res_radians / self.input_map.res
# band_res_pixels = band_res_radians / self.map.res
# FWHM_TO_SIGMA = 2.355
# band_beam_sigma_pixels = band_res_pixels / FWHM_TO_SIGMA

Expand All @@ -257,7 +251,7 @@ def _sample_maps(self):
det_mask = det_freq_response > -np.inf # -1e-3

samples = sp.interpolate.RegularGridInterpolator(
(self.input_map.x_side, self.input_map.y_side),
(self.map.x_side, self.map.y_side),
filtered_nu_map_data,
bounds_error=False,
fill_value=0,
Expand Down
10 changes: 5 additions & 5 deletions maria/noise/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,26 +8,26 @@ def _run(self):
self._simulate_noise()

def _simulate_noise(self):
self.data["noise"] = np.zeros((self.instrument.n_dets, self.pointing.n_time))
self.data["noise"] = np.zeros((self.instrument.n_dets, self.plan.n_time))

for band in self.instrument.dets.bands:
band_index = self.instrument.dets.subset(band=band.name).index

if band.white_noise > 0:
self.data["noise"][band_index] += (
np.sqrt(self.pointing.sample_rate)
np.sqrt(self.plan.sample_rate)
* band.white_noise
* np.random.standard_normal(
size=(len(band_index), self.pointing.n_time)
size=(len(band_index), self.plan.n_time)
)
)

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.pointing.n_time),
dt=self.pointing.dt,
size=int(self.plan.n_time),
dt=self.plan.dt,
)

def _spectrum_noise(self, spectrum_func, size, dt, amp=2.0):
Expand Down
Loading

0 comments on commit d9eb82c

Please sign in to comment.