Mock Observations with MUSTANG-2

We use Maria to generate synthetic mock observations for bolometric arrays installed on single-dish telescopes. The tool consists of four parts: An input map, an observing plan, an instrument, and a site. All are then combined in a wrapper named Simulation. As explained in the section “Usage”, we can emulate all kinds of telescopes, like ACT, Toltec, GBT, and AtLAST.

The simulation tool generates time-ordered domain (TOD) data. By simulating distinct atmospheric layers corresponding to the designated site within the simulation, we allow the mock array to sweep across the input celestial sky while scanning through the evolving atmosphere. As a result, the TODs produced combine elements of the real sky and atmospheric and detector noise.

This tutorial will guide you through the process of setting up the simulation tool and creating your own mock observations – In particular for GBT/MUSTANG-2 like observations of a galaxy cluster.

Setting up a MUSTANG-2 Simulation

maria offers flexibility for customizing atmospheric simulations. You can configure it to replicate observatories like ‘AtLAST’ at the APEX site or simulate specific observations, such as MUSTANG-2, under different weather conditions by adjusting the observing time.

Input map

First we specify the input map, our astronomical signal from which we want to make a mock observation. You can give it a fits file and load the map as outlined below. Don’t forgot to put the map at an RA and Dec (in degrees) in the sky. Also, the width or the pixel size (resolution) needs to be specified.

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(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.to(units="K_RJ").plot()

2024-10-17 04:20:26.564 INFO: created cache at /tmp/maria-data/maps
2024-10-17 04:20:28.670 INFO: downloaded data (4.0 MB) to /tmp/maria-data/maps/cluster.fits

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=600, # 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-10-17 04:20:45.085 INFO: Initialized base in 15594 ms.
2024-10-17 04:20:45.096 INFO: created cache at /tmp/maria-data/atmosphere/spectra/am
2024-10-17 04:20:45.820 INFO: downloaded data (11.4 MB) to /tmp/maria-data/atmosphere/spectra/am/chajnantor.h5
2024-10-17 04:20:46.227 INFO: created cache at /tmp/maria-data/atmosphere/weather/era5
2024-10-17 04:20:47.143 INFO: downloaded data (3.8 MB) to /tmp/maria-data/atmosphere/weather/era5/chajnantor.h5
Building atmosphere: 100%|██████████| 6/6 [00:06<00:00,  1.05s/it]
2024-10-17 04:20:54.719 INFO: Initialized atmosphere in 9622 ms.
Generating CMB:   0%|          | 0/1 [00:00<?, ?it/s]2024-10-17 04:20:54.722 INFO: created cache at /tmp/maria-data/cmb/spectra
2024-10-17 04:20:55.039 INFO: downloaded data (0.3 MB) to /tmp/maria-data/cmb/spectra/planck.csv
Generating CMB: 100%|██████████| 1/1 [00:05<00:00,  5.84s/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, 992.70it/s]
Sampling turbulence: 100%|██████████| 6/6 [00:00<00:00,  7.05it/s]
Computing atmospheric emission (m2/f093): 100%|██████████| 1/1 [00:00<00:00,  4.60it/s]
Computing atmospheric transmission (m2/f093): 100%|██████████| 1/1 [00:00<00:00,  4.99it/s]
Sampling CMB (m2/f093): 100%|██████████| 1/1 [00:00<00:00, 12.27it/s]
Sampling map (m2/f093): 100%|██████████| 1/1 [00:03<00:00,  3.22s/it]
Generating noise: 100%|██████████| 1/1 [00:00<00:00,  2.53it/s]

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.map.mappers import BinMapper

mapper = BinMapper(center=(150, 10),
                   frame="ra_dec",
                   width=0.1,
                   height=0.1,
                   resolution=2e-4,
                   tod_preprocessing={
                        "remove_modes": {"n": 1},
                        # "highpass": {"f": 0.01},
                        "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:06<00:00,  6.70s/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')