Polarized observations¶
This tutorial covers working with polarized instrument and maps, and recovering polarized maps from observations.
We start with a normal instrument, and create two orthogonally polarized copies of each detector by setting polarized: True in the Array config:
[1]:
import maria
from maria.instrument import Band
f150 = Band(
center=150e9,
width=30e9,
NET_RJ=60e-6,
knee=1e0,
gain_error=5e-2)
array = {"field_of_view": 0.2,
"shape": "circle",
"beam_spacing": 1.8,
"primary_size": 10,
"polarized": True,
"packing": "sunflower",
"bands": [f150]}
instrument = maria.get_instrument(array=array)
print(instrument.arrays)
n field_of_view max_baseline bands polarized primary_size
array1 88 11.98’ 0 m [f150] True 10 m
We can see the resulting polarization footprint in the instrument plot:
[2]:
print(instrument)
instrument.plot()
Instrument(1 array)
├ arrays:
│ n field_of_view max_baseline bands polarized primary_size
│ array1 88 11.98’ 0 m [f150] True 10 m
│
└ bands:
name center width η NEP NET_RJ NET_CMB FWHM
0 f150 150 GHz 30 GHz 0.5 13.22 aW√s 60 uK_RJ√s 104 uK_CMB√s 52.49”
Let’s observe the use the Einstein map, which has a faint polarization signature underneath the unpolarized signal of Einstein’s face. Remember that all maps are five dimensional (stokes, frequency, time, y, x); this map has four channels in the stokes dimensions (the I, Q, U, and V Stokes parameters). We can plot all the channels by giving plot a shaped set of stokes parameters.
[3]:
input_map = maria.map.get("maps/einstein.h5")
print(input_map)
input_map.plot(slices={"stokes": [["I", "Q"],
["U", "V"]]})
2026-07-06 17:31:53.010 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/maps/einstein.h5
Downloading: 100%|██████████| 932k/932k [00:00<00:00, 46.8MB/s]
ProjectionMap:
data(4, 1, 685, 685):
min: -1.000e-02
max: 2.540e-01
units: K_RJ
quantity: rayleigh_jeans_temperature
stokes(4):
components: ['I' 'Q' 'U' 'V']
nu(1):
values: [90.] GHz
eta(685):
height: 60’
res: -5.263”
xi(685):
width: 60’
res: 5.263”
frame: ra/dec
center:
ra: 00ʰ00ᵐ0.00ˢ
dec: -23°00’0.00”
beam(maj, min, psi): (0 rad, 0 rad, 0 rad)
memory: 7.508 MB
[4]:
from maria import Planner
planner = Planner(target=input_map, site="mauna_kea", constraints={"el": (45, 90)})
plans = planner.generate_plans(
total_duration=1800,
max_chunk_duration=1800,
sample_rate=50,
scan_type="raster",
scan_parameters={
"n": [(15, 1), (1, 16)],
"radius": 0.5 * input_map.width.deg,
"speed": 0.5,
})
plans[0].plot()
print(plans)
PlanList(1 plans, 1800 s):
start_time duration sample_rate target(ra,dec) center(az,el)
chunk
0 2026-07-07 14:22:34.762 +00:00 1800 s 50 Hz (360°, -23°) (165.5°, 46.02°)
[5]:
sim = maria.Simulation(
instrument,
plans=plans,
site="mauna_kea",
atmosphere="2d",
map=input_map,
)
print(sim)
Initializing observations: 0%| | 0/1 [00:00<?, ?it/s]2026-07-06 17:32:04.379 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/spectra/am/v3/mauna_kea.h5
Downloading: 0%| | 0.00/22.0M [00:00<?, ?B/s]
Downloading: 100%|██████████| 22.0M/22.0M [00:00<00:00, 167MB/s]
2026-07-06 17:32:05.773 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/weather/era5/mauna_kea.h5
Downloading: 100%|██████████| 8.62M/8.62M [00:00<00:00, 122MB/s]
Constructing atmosphere: 0%| | 0/8 [00:00<?, ?it/s]
Constructing atmosphere: 12%|█▎ | 1/8 [00:06<00:42, 6.14s/it]
Constructing atmosphere: 25%|██▌ | 2/8 [00:12<00:38, 6.43s/it]
Constructing atmosphere: 38%|███▊ | 3/8 [00:18<00:31, 6.29s/it]
Constructing atmosphere: 50%|█████ | 4/8 [00:23<00:22, 5.63s/it]
Constructing atmosphere: 62%|██████▎ | 5/8 [00:29<00:17, 5.93s/it]
Constructing atmosphere: 75%|███████▌ | 6/8 [00:33<00:10, 5.21s/it]
Constructing atmosphere: 88%|████████▊ | 7/8 [00:38<00:04, 4.96s/it]
Constructing atmosphere: 100%|██████████| 8/8 [00:41<00:00, 5.24s/it]
Initializing observations: 100%|██████████| 1/1 [00:45<00:00, 45.96s/it]
Simulation
├ Instrument(1 array)
│ ├ arrays:
│ │ n field_of_view max_baseline bands polarized primary_size
│ │ array1 88 11.98’ 0 m [f150] True 10 m
│ │
│ └ bands:
│ name center width η NEP NET_RJ NET_CMB FWHM
│ 0 f150 150 GHz 30 GHz 0.5 13.22 aW√s 60 uK_RJ√s 104 uK_CMB√s 52.49”
├ Site:
│ region: mauna_kea
│ timezone: Pacific/Honolulu
│ location:
│ longitude: 155°28’37.20” W
│ latitude: 19°49’22.08” N
│ altitude: 4.092 km
│ seasonal: True
│ diurnal: True
├ PlanList(1 plans, 1800 s):
│ start_time duration sample_rate target(ra,dec) center(az,el)
│ chunk
│ 0 2026-07-07 14:22:34.762 +00:00 1800 s 50 Hz (360°, -23°) (165.5°, 46.02°)
├ Atmosphere(8 processes with 8 layers):
│ ├ spectrum:
│ │ region: mauna_kea
│ └ weather:
│ region: mauna_kea
│ altitude: 4.092 km
│ time: Jul 7 04:37:34 -10:00
│ pwv[mean, rms]: (2.204 mm, 66.12 um)
[6]:
tods = sim.run()
print(tods)
tods[0].plot()
2026-07-06 17:32:49.616 INFO: Simulating observation 1 of 1
Generating turbulence: 100%|██████████| 8/8 [00:02<00:00, 3.21it/s]
Sampling turbulence: 100%|██████████| 8/8 [00:04<00:00, 1.66it/s]
Computing atmospheric emission: 100%|██████████| 1/1 [00:01<00:00, 1.09s/it, band=f150]
Sampling source 'map': 100%|██████████| 1/1 [00:06<00:00, 6.65s/it, band=f150, message=Sampling channel (105 GHz, 195 GHz)]
Generating noise: 100%|██████████| 1/1 [00:00<00:00, 1.03it/s, band=f150]
2026-07-06 17:33:07.978 INFO: Simulated observation 1 of 1 in 18.35 s
[TOD(shape=(88, 90000), fields=['atmosphere', 'map', 'noise'], units='uK_RJ', start=2026-07-07 14:52:34.740 +00:00, duration=1800.0s, sample_rate=50.0Hz, metadata={'atmosphere': True, 'sim_time': <Arrow [2026-07-06T17:33:06.040030+00:00]>, 'altitude': 4092.0, 'region': 'mauna_kea', 'pwv': 2.204, 'base_temperature': 279.407})]
[7]:
from maria.mapping import MaximumLikelihoodMapper
mapper = MaximumLikelihoodMapper(
frame="ra/dec",
resolution=2 * input_map.resolution.deg,
tod_preprocessing={
"remove_spline": {"knot_spacing": 300, "remove_el_gradient_order": 1},
},
tods=tods,
bilinear=False,
)
mapper.map.plot(slices="all")
2026-07-06 17:33:13.510 INFO: Inferring center {'ra': '23ʰ59ᵐ59.39ˢ', 'dec': '-22°59’58.15”'} for mapper
2026-07-06 17:33:13.522 INFO: Inferring mapper width 1.198° for mapper from observation patch
2026-07-06 17:33:13.523 INFO: Inferring mapper height 1.198° to match supplied width
2026-07-06 17:33:16.917 INFO: Inferring stokes parameters 'IQU' for mapper from detector sensitivities
Preprocessing TODs: 100%|██████████| 1/1 [00:03<00:00, 3.25s/it]
Computing pointing matrices: 100%|██████████| 1/1 [00:03<00:00, 3.67s/it]
[8]:
mapper.fit(epochs=2,
steps_per_epoch=25,
plot=True,
plot_kwargs={
"slices": "all",
"center_zero": True,
"contrast": 1e-3})
Updating noise model: 100%|██████████| 1/1 [00:05<00:00, 5.15s/it, tod=1/1]
Fitting map (epoch 1/2): 231it [13:53, 3.80s/it, alpha=0.0348]2026-07-06 17:47:34.144 INFO: Finished conjugate gradient descent, terminating
Fitting map (epoch 1/2): 231it [13:57, 3.62s/it, alpha=0.0348]
Updating noise model: 100%|██████████| 1/1 [00:05<00:00, 5.57s/it, tod=1/1]
Fitting map (epoch 2/2): 254it [14:58, 3.68s/it, alpha=0.00806]2026-07-06 18:02:43.712 INFO: Finished conjugate gradient descent, terminating
Fitting map (epoch 2/2): 254it [15:02, 3.55s/it, alpha=0.00806]
Note that we can’t see any of the circular polarization, since our instrument isn’t sensitive to it.