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:

import maria
from maria.instrument import Band

f090 = Band(
    center=90e9,  # in Hz
    width=20e9,  # in Hz
    NET_RJ=40e-6,  # in K sqrt(s)
    knee=1e0,    # in Hz
    gain_error=5e-2)

f150 = Band(
    center=150e9,
    width=30e9,
    NET_RJ=60e-6,
    knee=1e0,
    gain_error=5e-2)

array = {"field_of_view": 0.5,
         "shape": "circle",
         "beam_spacing": 1.5,
         "primary_size": 10,
         "polarized": True,
         "bands": [f090, f150]}

instrument = maria.get_instrument(array=array)

print(instrument.arrays)

          n     FOV baseline        bands polarized
array1  652  28.68’      0 m  [f090,f150]      True

We can see the resulting polarization footprint in the instrument plot:

print(instrument)
instrument.plot()

Instrument(1 array)
├ arrays:
│            n     FOV baseline        bands polarized
│  array1  652  28.68’      0 m  [f090,f150]      True
│
└ bands:
      name   center   width    η         NEP      NET_RJ         NET_CMB    FWHM
   0  f090   90 GHz  20 GHz  0.5  5.445 aW√s  40 uK_RJ√s  49.13 uK_CMB√s  1.458’
   1  f150  150 GHz  30 GHz  0.5  12.25 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.

input_map = maria.map.get("maps/einstein.h5")
input_map.plot(stokes=[["I", "Q"],
                       ["U", "V"]])

2026-03-20 19:42:59.948 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/maps/einstein.h5
Downloading: 100%|████████████████| 931k/931k [00:00<00:00, 38.3MB/s]

from maria import Planner

planner = Planner(target=input_map, site="llano_de_chajnantor", constraints={"el": (60, 90)})
plans = planner.generate_plans(total_duration=900,  # in seconds
                               sample_rate=50)  # in Hz

plans[0].plot()
print(plans)

PlanList(1 plans, 900 s):
                           start_time duration target(ra,dec)     center(az,el)
chunk
0      2026-03-21 14:26:11.406 +00:00    900 s   (360°, -23°)  (95.74°, 61.92°)

sim = maria.Simulation(
    instrument,
    plans=plans,
    site="llano_de_chajnantor",
    atmosphere="2d",
    atmosphere_kwargs={"weather": {"pwv": 0.5}},
    map=input_map
)

print(sim)

Simulation
├ Instrument(1 array)
│ ├ arrays:
│ │            n     FOV baseline        bands polarized
│ │  array1  652  28.68’      0 m  [f090,f150]      True
│ │
│ └ bands:
│       name   center   width    η         NEP      NET_RJ         NET_CMB    FWHM
│    0  f090   90 GHz  20 GHz  0.5  5.445 aW√s  40 uK_RJ√s  49.13 uK_CMB√s  1.458’
│    1  f150  150 GHz  30 GHz  0.5  12.25 aW√s  60 uK_RJ√s    104 uK_CMB√s  52.49”
├ Site:
│   region: chajnantor
│   timezone: America/Santiago
│   location:
│     longitude: 67°45’17.28” W
│     latitude:  23°01’45.84” S
│     altitude: 5.064 km
│   seasonal: True
│   diurnal: True
├ PlanList(1 plans, 900 s):
│                            start_time duration target(ra,dec)     center(az,el)
│ chunk
│ 0      2026-03-21 14:26:11.406 +00:00    900 s   (360°, -23°)  (95.74°, 61.92°)
├ '2d'
└ ProjectionMap:
    shape(stokes, nu, t, y, x): (4, 1, 1, 685, 685)
    stokes: IQUV
    nu: [90.] GHz
    t: [1.77403578e+09]
    z: naive
    quantity: rayleigh_jeans_temperature
    units: K_RJ
      min: -1.000e-02
      max: 2.540e-01
      rms: 5.754e-02
    center:
      ra: 00ʰ00ᵐ0.00ˢ
      dec: -23°00’0.00”
    size(y, x): (1°, 1°)
    resolution(y, x): (5.255”, 5.255”)
    beam(maj, min, rot): (0 rad, 0 rad, 0 rad)
    memory: 30.03 MB

tods = sim.run()

print(tods)
tods[0].plot()

2026-03-20 19:43:12.271 INFO: Simulating observation 1 of 1
Constructing atmosphere: 100%|████████████████| 8/8 [00:11<00:00,  1.40s/it]
Generating turbulence: 100%|████████████████| 8/8 [00:01<00:00,  5.51it/s]
Sampling turbulence: 100%|████████████████| 8/8 [00:06<00:00,  1.17it/s]
Computing atmospheric emission: 100%|████████████████| 2/2 [00:02<00:00,  1.07s/it, band=f150]
Sampling map:   0%|                | 0/2 [00:00<?, ?it/s, band=f090]               /opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:4596: RuntimeWarning: invalid value encountered in scalar subtract
  diff_b_a = b - a
Sampling map:  50%|████████        | 1/2 [00:07<00:07,  7.76s/it, band=f150]                        /opt/hostedtoolcache/Python/3.12.13/x64/lib/python3.12/site-packages/numpy/lib/_function_base_impl.py:4596: RuntimeWarning: invalid value encountered in scalar subtract
  diff_b_a = b - a
Sampling map: 100%|████████████████| 2/2 [00:14<00:00,  7.14s/it, band=f150, channel=(0 Hz, inf Hz)]
Generating noise: 100%|████████████████| 2/2 [00:01<00:00,  1.61it/s, band=f150]
2026-03-20 19:44:08.405 INFO: Simulated observation 1 of 1 in 56.12 s

[TOD(shape=(652, 45000), fields=['atmosphere', 'map', 'noise'], units='K_RJ', start=2026-03-21 14:41:11.385 +00:00, duration=900.0s, sample_rate=50.0Hz, metadata={'atmosphere': True, 'sim_time': <Arrow [2026-03-20T19:43:52.165263+00:00]>, 'altitude': 5064.0, 'region': 'chajnantor', 'pwv': 0.5, 'base_temperature': 279.05})]

from maria.mappers import BinMapper

mapper = BinMapper(
    stokes="IQUV",
    frame="ra/dec",
    resolution=0.25 / 60,
    tod_preprocessing={
        "remove_spline": {"knot_spacing": 60, "remove_el_gradient": True},
        "remove_modes": {"modes_to_remove": 1},
    },
    map_postprocessing={
        "gaussian_filter": {"sigma": 1},
    },
    units="mK_RJ",
    tods=tods,
)

output_map = mapper.run()

print(output_map)

2026-03-20 19:44:18.646 INFO: Inferring center {'ra': '23ʰ59ᵐ59.71ˢ', 'dec': '-23°00’8.15”'} for mapper.
2026-03-20 19:44:18.665 INFO: Inferring mapper width 1.459° for mapper from observation patch.
2026-03-20 19:44:18.666 INFO: Inferring mapper height 1.459° to match supplied width.
Preprocessing TODs: 100%|████████████████| 1/1 [00:02<00:00,  2.14s/it]
Mapping: 100%|██████████| 1/1 [00:04<00:00,  4.69s/it, tod=1/1]
2026-03-20 19:44:26.471 WARNING: No counts for map (stokes=V, nu=90 GHz)
2026-03-20 19:44:26.486 WARNING: No counts for map (stokes=V, nu=150 GHz)

ProjectionMap:
  shape(stokes, nu, t, y, x): (4, 2, 1, 350, 350)
  stokes: IQUV
  nu: [ 90. 150.] GHz
  t: [1.77410362e+09]
  z: naive
  quantity: rayleigh_jeans_temperature
  units: mK_RJ
    min: -1.135e+02
    max: 1.404e+02
    rms: 2.670e+01
  center:
    ra: 23ʰ59ᵐ59.71ˢ
    dec: -23°00’8.15”
  size(y, x): (1.458°, 1.458°)
  resolution(y, x): (15”, 15”)
  beam(maj, min, rot): ragged
  memory: 15.68 MB

Note that we can’t see any of the circular polarization, since our instrument isn’t sensitive to it.

output_map.plot(stokes=["I", "Q", "U"], nu_index=[[0], [1]])