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 Crab Nebula (M1).
[1]:
import maria
input_map = maria.map.get("maps/M1.h5").to("K_RJ")
input_map.data[input_map.weight < 0.2 * input_map.weight.max()] = 0
input_map.plot(slices="all")
print(input_map)
2026-07-06 17:29:12.090 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/maps/M1.h5
Downloading: 100%|██████████| 21.8M/21.8M [00:00<00:00, 129MB/s]
ProjectionMap:
data(1, 3, 1205, 1187):
min: -6.150e-02
max: 5.613e-01
units: K_RJ
quantity: rayleigh_jeans_temperature
stokes(1):
components: ['I']
nu(3):
values: [149.8962 214.1375 272.5386] GHz
eta(1205):
height: 20.07’
res: -1”
xi(1187):
width: 19.77’
res: 1”
frame: ra/dec
center:
ra: 05ʰ34ᵐ31.95ˢ
dec: 22°00’52.16”
beam(maj, min, psi): ragged
memory: 34.33 MB
[2]:
from maria import Planner
planner = Planner(target=input_map, site="green_bank", constraints={"el": (60, 90)})
plans = planner.generate_plans(total_duration=600, sample_rate=50)
plans[0].plot()
print(plans)
PlanList(1 plans, 600 s):
start_time duration sample_rate target(ra,dec) center(az,el)
chunk
0 2026-07-07 13:57:25.731 +00:00 600 s 50 Hz (83.63°, 22.02°) (116.6°, 61.17°)
[3]:
instrument = maria.get_instrument("MUSTANG-2")
print(instrument)
instrument.plot()
Instrument(1 array)
├ arrays:
│ n field_of_view max_baseline bands polarized primary_size
│ array1 217 4.2’ 0 m [m2/f093] False 100 m
│
└ bands:
name center width η NEP NET_RJ NET_CMB FWHM
0 m2/f093 86.21 GHz 20.98 GHz 0.1 15 aW√s 571.1 uK_RJ√s 690.5 uK_CMB√s 9.133”
[4]:
sim = maria.Simulation(
instrument,
plans=plans,
site="green_bank",
map=input_map,
atmosphere="2d",
)
print(sim)
Initializing observations: 0%| | 0/1 [00:00<?, ?it/s]2026-07-06 17:29:25.696 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/spectra/am/v3/green_bank.h5
Downloading: 0%| | 0.00/22.0M [00:00<?, ?B/s]
Downloading: 100%|██████████| 22.0M/22.0M [00:00<00:00, 136MB/s]
2026-07-06 17:29:26.451 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/weather/era5/green_bank.h5
Downloading: 100%|██████████| 12.0M/12.0M [00:00<00:00, 151MB/s]
Constructing atmosphere: 0%| | 0/8 [00:00<?, ?it/s]
Constructing atmosphere: 25%|██▌ | 2/8 [00:00<00:00, 13.29it/s]
Constructing atmosphere: 50%|█████ | 4/8 [00:00<00:00, 13.00it/s]
Constructing atmosphere: 75%|███████▌ | 6/8 [00:00<00:00, 13.26it/s]
Constructing atmosphere: 100%|██████████| 8/8 [00:00<00:00, 13.12it/s]
Initializing observations: 100%|██████████| 1/1 [00:03<00:00, 3.97s/it]
Simulation
├ Instrument(1 array)
│ ├ arrays:
│ │ n field_of_view max_baseline bands polarized primary_size
│ │ array1 217 4.2’ 0 m [m2/f093] False 100 m
│ │
│ └ bands:
│ name center width η NEP NET_RJ NET_CMB FWHM
│ 0 m2/f093 86.21 GHz 20.98 GHz 0.1 15 aW√s 571.1 uK_RJ√s 690.5 uK_CMB√s 9.133”
├ Site:
│ region: green_bank
│ timezone: America/New_York
│ location:
│ longitude: 79°50’23.28” W
│ latitude: 38°25’59.16” N
│ altitude: 825 m
│ seasonal: True
│ diurnal: True
├ PlanList(1 plans, 600 s):
│ start_time duration sample_rate target(ra,dec) center(az,el)
│ chunk
│ 0 2026-07-07 13:57:25.731 +00:00 600 s 50 Hz (83.63°, 22.02°) (116.6°, 61.17°)
├ Atmosphere(8 processes with 8 layers):
│ ├ spectrum:
│ │ region: green_bank
│ └ weather:
│ region: green_bank
│ altitude: 825 m
│ time: Jul 7 10:02:25 -04:00
│ pwv[mean, rms]: (27.79 mm, 833.7 um)
[5]:
tods = sim.run()
tods[0].plot()
2026-07-06 17:29:29.025 INFO: Simulating observation 1 of 1
Generating turbulence: 100%|██████████| 8/8 [00:00<00:00, 49.73it/s]
Sampling turbulence: 100%|██████████| 8/8 [00:03<00:00, 2.05it/s]
Computing atmospheric emission: 100%|██████████| 1/1 [00:00<00:00, 1.47it/s, band=m2/f093]
Sampling source 'map': 100%|██████████| 1/1 [00:06<00:00, 6.11s/it, band=m2/f093, message=Sampling channel (74 GHz, 105 GHz)]
Generating noise: 100%|██████████| 1/1 [00:00<00:00, 1.04it/s, band=m2/f093]
2026-07-06 17:29:42.836 INFO: Simulated observation 1 of 1 in 13.8 s
[6]:
from maria.mappers import MaximumLikelihoodMapper
mapper = MaximumLikelihoodMapper(
units="K_RJ",
tods=tods,
resolution=2 * input_map.resolution,
)
mapper.map.plot()
2026-07-06 17:29:47.769 INFO: Inferring center {'ra': '05ʰ34ᵐ32.21ˢ', 'dec': '22°00’56.30”'} for mapper
2026-07-06 17:29:47.781 INFO: Inferring mapper width 23.98’ for mapper from observation patch
2026-07-06 17:29:47.782 INFO: Inferring mapper height 23.98’ to match supplied width
2026-07-06 17:29:50.620 INFO: Inferring stokes parameters 'I' for mapper from detector sensitivities
Preprocessing TODs: 100%|██████████| 1/1 [00:02<00:00, 2.49s/it]
Computing pointing matrices: 100%|██████████| 1/1 [00:01<00:00, 1.55s/it]
[7]:
mapper.fit(epochs=1, max_steps_per_epoch=100, plot=True)
Updating noise model: 100%|██████████| 1/1 [00:01<00:00, 1.57s/it, tod=1/1]
Fitting map (epoch 1/1): 100it [01:42, 1.02s/it, alpha=0.15]
[8]:
from maria.mappers import compute_residual_map
residual_map = compute_residual_map(input_map.to("K_RJ")[:, 0], mapper.map)
residual_map.plot()
