Simulating observations with MUSTANG-2

MUSTANG-2 is a bolometric array on the Green Bank Telescope. In this notebook we simulate an observation of the Whirlpool Galaxy (M51).

import maria
from maria.io import fetch

input_map = maria.map.load(fetch("maps/M51HA.fits"), nu=93e9)

print(input_map)
input_map.to("K_RJ").plot()

Downloading https://github.com/thomaswmorris/maria-data/raw/master/maps/M51HA.fits: 100%|██████████| 1.06M/1.06M [00:00<00:00, 66.3MB/s]

ProjectedMap:
  shape(nu, y, x): (1, 512, 512)
  stokes: naive
  nu: [93.] GHz
  t: naive
  z: naive
  quantity: spectral_flux_density_per_beam
  units: Jy/beam
    min: 0.000e+00
    max: 1.361e+01
  center:
    ra: 13ʰ29ᵐ5.31ˢ
    dec: 47°11’43.00”
  size(y, x): (4.267’, 4.267’)
  resolution(y, x): (0.5”, 0.5”)
  beam(maj, min, psi): (3.36”, 3.181”, -89.2°)
  memory: 4.194 MB

plan = maria.get_plan(
    scan_pattern="daisy",  # scanning pattern
    scan_options={"radius": 2 / 60, "speed": 0.5 / 60},  # in degrees
    duration=600,  # integration time in seconds
    sample_rate=50,  # in Hz
    scan_center=(202.27211, 47.195277),  # position in the sky
    frame="ra_dec",
)

plan.plot()

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

Instrument(1 array)
├ arrays:
│            n   FOV baseline      bands polarized
│  array1  217  4.2’      0 m  [m2/f093]     False
│
└ bands:
         name     center      width    η      NEP       NET_RJ      NET_CMB  \
   0  m2/f093  86.21 GHz  20.98 GHz  0.1  15 aW√s  0.5711 mK√s  0.6905 mK√s

        FWHM
   0  9.133”

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

print(sim)

Downloading https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/spectra/am/v2/green_bank.h5: 100%|██████████| 11.1M/11.1M [00:00<00:00, 195MB/s]
Downloading https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/weather/era5/green_bank.h5: 100%|██████████| 12.0M/12.0M [00:00<00:00, 169MB/s]
Constructing atmosphere: 100%|██████████| 10/10 [00:04<00:00,  2.38it/s]

Simulation
├ Instrument(1 array)
│ ├ arrays:
│ │            n   FOV baseline      bands polarized
│ │  array1  217  4.2’      0 m  [m2/f093]     False
│ │
│ └ bands:
│          name     center      width    η      NEP       NET_RJ      NET_CMB  \
│    0  m2/f093  86.21 GHz  20.98 GHz  0.1  15 aW√s  0.5711 mK√s  0.6905 mK√s
│
│         FWHM
│    0  9.133”
├ Site:
│   region: green_bank
│   location: 38°25’59.16”N 79°50’23.28”W
│   altitude: 0.825 km
│   seasonal: True
│   diurnal: True
├ Plan:
│   start_time: 2024-02-10 06:00:00.000 +00:00
│   duration: 600 s
│   sample_rate: 50 Hz
│   center:
│     ra: 13ʰ29ᵐ5.31ˢ
│     dec: 47°11’43.00”
│   scan_pattern: daisy
│   scan_radius: 3.996’
│   scan_kwargs: {'radius': 0.03333333333333333, 'speed': 0.008333333333333333}
├ Atmosphere(10 processes with 10 layers):
│ ├ spectrum:
│ │   region: green_bank
│ └ weather:
│     region: green_bank
│     altitude: 0.825 km
│     time: Feb 10 01:04:59 -05:00
│     pwv[mean, rms]: (8.091 mm, 242.7 um)
└ ProjectedMap:
    shape(stokes, nu, y, x): (1, 1, 512, 512)
    stokes: I
    nu: [93.] GHz
    t: naive
    z: naive
    quantity: rayleigh_jeans_temperature
    units: K_RJ
      min: 0.000e+00
      max: 6.487e+01
    center:
      ra: 13ʰ29ᵐ5.31ˢ
      dec: 47°11’43.00”
    size(y, x): (4.267’, 4.267’)
    resolution(y, x): (0.5”, 0.5”)
    beam(maj, min, psi): (3.36”, 3.181”, -89.2°)
    memory: 4.194 MB

tod = sim.run()
tod.plot()

Generating turbulence: 100%|██████████| 10/10 [00:00<00:00, 41.12it/s]
Sampling turbulence: 100%|██████████| 10/10 [00:03<00:00,  3.28it/s]
Computing atmospheric emission: 100%|██████████| 1/1 [00:01<00:00,  1.67s/it, band=m2/f093]
Sampling map: 100%|██████████| 1/1 [00:02<00:00,  2.41s/it, channel=[ 0. inf] Hz]
Generating noise: 100%|██████████| 1/1 [00:00<00:00,  2.79it/s, band=m2/f093]

from maria.mappers import BinMapper

mapper = BinMapper(
    center=(202.27211, 47.195277),
    frame="ra_dec",
    width=6 / 60,
    height=6 / 60,
    resolution=0.03 / 60,
    tod_preprocessing={
        "window": {"name": "hamming"},
        "remove_modes": {"modes_to_remove": [0]},
        "remove_spline": {"knot_spacing": 30, "remove_el_gradient": True},
    },
    map_postprocessing={
        "gaussian_filter": {"sigma": 1},
        "median_filter": {"size": 1},
    },
    units="uK_RJ",
)

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

Mapping band m2/f093: 100%|██████████| 1/1 [00:00<00:00, 24.56it/s, band=m2/f093, stokes=I]
2025-05-30 18:52:07.444 INFO: Ran mapper for band m2/f093 in 2.847 s.

output_map.plot()
output_map.to_fits("/tmp/simulated_mustang_map.fits")