Simulating observations with MUSTANG-2

MUSTANG-2 is a bolometric array on the Green Bank Telescope. In this notebook we simulate an observation of a galactic cluster.

import maria
import scipy as sp
import numpy as np
import matplotlib.pyplot as plt

map_filename = maria.io.fetch("maps/cluster.fits")

# load in the map from a fits file
input_map = maria.map.read_fits(
    nu=150,
    filename=map_filename,  # filename
    resolution=8.714e-05,  # pixel size in degrees
    index=0,  # index for fits file
    center=(150, 10),  # position in the sky
    units="Jy/pixel",  # Units of the input map
)

input_map.data *= 1e1
input_map.to(units="K_RJ").plot()

100%|██████████| 12.6M/12.6M [00:00<00:00, 133MB/s]
2024-12-14 03:50:36.642 INFO: Downloading https://github.com/thomaswmorris/maria-data/raw/master/maps/cluster.fits to /tmp/maria-data/maps …
2024-12-14 03:50:37.400 INFO: Caching data from https://github.com/thomaswmorris/maria-data/raw/master/maps/cluster.fits
100%|██████████| 4.00M/4.00M [00:00<00:00, 94.0MB/s]

Observing strategy

Below we define the observing strategy containing the duration, integration time, read out rate, and pointing center

# load the map into maria
plan = maria.get_plan(
    scan_pattern="daisy",  # scanning pattern
    scan_options={"radius": 0.05, "speed": 0.01},  # in degrees
    duration=1200,  # integration time in seconds
    sample_rate=50,  # in Hz
    scan_center=(150, 10),  # position in the sky
    frame="ra_dec",
)

plan.plot()

Define instrument

We have configuration files in place to mimic typical telescopes and instruments. To run a MUSTANG-2 like observations, simply initialize the instrument as follows. Other tutorials will go into more detail on how to adjust it to test out how telescope designs affect the recovered signal.

instrument = maria.get_instrument("MUSTANG-2")

We can plot the instruments which will show the location of the detectors in the FoV with their diffrection limited beam size. The detectors for MUSTANG-2 are spaced at 2 f-lambda.

instrument.plot()

Initialize the simulation

The simulation class combines all the components into one.

sim = maria.Simulation(
    instrument,
    plan=plan,
    site="llano_de_chajnantor",
    map=input_map,
    atmosphere="2d",
    cmb="generate",
)

2024-12-14 03:50:50.017 INFO: Initialized base in 11128 ms.
2024-12-14 03:50:50.037 INFO: Downloading https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/spectra/am/chajnantor.h5 to /tmp/maria-data/atmosphere/spectra/am …
2024-12-14 03:50:50.190 INFO: Caching data from https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/spectra/am/chajnantor.h5
100%|██████████| 11.4M/11.4M [00:00<00:00, 123MB/s]
2024-12-14 03:50:50.648 INFO: Downloading https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/weather/era5/chajnantor.h5 to /tmp/maria-data/atmosphere/weather/era5 …
2024-12-14 03:50:51.201 INFO: Caching data from https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/weather/era5/chajnantor.h5
100%|██████████| 4.05M/4.05M [00:00<00:00, 95.6MB/s]
Building atmosphere: 100%|██████████| 6/6 [00:10<00:00,  1.72s/it]
2024-12-14 03:51:03.680 INFO: Initialized atmosphere in 13642 ms.
Generating CMB:   0%|          | 0/1 [00:00<?, ?it/s]2024-12-14 03:51:03.682 INFO: Downloading https://github.com/thomaswmorris/maria-data/raw/master/cmb/spectra/planck.csv to /tmp/maria-data/cmb/spectra …
2024-12-14 03:51:04.045 INFO: Caching data from https://github.com/thomaswmorris/maria-data/raw/master/cmb/spectra/planck.csv

293kB [00:00, 30.1MB/s]
Generating CMB: 100%|██████████| 1/1 [00:05<00:00,  5.97s/it]

Obtaining time-ordered data (TODs)

To acquire time-ordered data (TODs), you only need to run the script. The TOD object contains time stamps, coordinates (in RA and Dec), and the intensity of each scan for every detector. The intensity units are expressed in power can be converted to surface brightness temperature units: Kelvin Rayleigh-Jeans (Even if the input map is given in Jy/pixel. We convert the units under the hood).

tod = sim.run()

Generating turbulence: 100%|██████████| 6/6 [00:00<00:00, 205.27it/s]
Sampling turbulence: 100%|██████████| 6/6 [00:01<00:00,  3.55it/s]
Computing atmospheric emission: 100%|██████████| 1/1 [00:00<00:00,  2.60it/s, band=m2/f093]
Computing atmospheric opacity: 100%|██████████| 1/1 [00:00<00:00,  2.38it/s, band=m2/f093]
Sampling CMB: 100%|██████████| 1/1 [00:00<00:00,  5.66it/s, band=m2/f093]
Sampling map: 100%|██████████| 1/1 [00:05<00:00,  5.45s/it, band=m2/f093]
Generating noise: 100%|██████████| 1/1 [00:00<00:00,  1.37it/s, band=m2/f093]

tod.fields

['atmosphere', 'cmb', 'map', 'noise']

Useful visualizations:

Now, let’s visualize the TODS. The figure shows the mean powerspectra and a time series of the observations.

tod.plot()

Mapping the TODs

To transform the TODs into images, you’ll require a mapper. While there is an option to use your own mapper, we have also implemented one for your convenience. This built-in mapper follows a straightforward approach: it removes the common mode over all detectors (the atmosphere if you scan fast enough), then detrends and optionally Fourier filters the time-streams to remove any noise smaller than the resolution of the dish. Then, the mapper bins the TOD on the specified grid. Notably, the mapper is designed to effectively eliminate correlated atmospheric noise among different scans, thereby increasing SNR of the output map.

An example of how to run the mapper on the TOD is as follows:

from maria.mappers import BinMapper

mapper = BinMapper(
    center=(150, 10),
    frame="ra_dec",
    width=0.1,
    height=0.1,
    resolution=2e-4,
    tod_preprocessing={
        "window": {"name": "hamming"},
        "remove_modes": {"modes_to_remove": [0]},
        "despline": {"knot_spacing": 10},
    },
    map_postprocessing={
        "gaussian_filter": {"sigma": 1},
        "median_filter": {"size": 1},
    },
)

mapper.add_tods(tod)
output_map = mapper.run()

Running mapper (m2/f093): 100%|██████████| 1/1 [00:14<00:00, 14.66s/it]

Visualizing the maps

As interesting as the detector setup, power spectra, and time series are, intuitively we think in the image plane. so let’s visualize it!

output_map.plot()

You can also save the map to a fits file. Here, we recover the units in which the map was initially specified in.

output_map.to_fits("/tmp/simulated_map.fits")