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

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)
instrument.arrays
[1]:
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:

[2]:
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√s  49.13 uK√s   1.458’
   1  f150  150 GHz  30 GHz  0.5  12.25 aW√s  60 uK√s    104 uK√s  0.8748’
../_images/tutorials_polarized-observations_4_1.png

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]:
from maria.io import fetch

input_map = maria.map.load(fetch("maps/einstein.h5"))
input_map.plot(stokes=[["I", "Q"],
                       ["U", "V"]])
Downloading https://github.com/thomaswmorris/maria-data/raw/master/maps/einstein.h5: 100%|██████████| 931k/931k [00:00<00:00, 95.7MB/s]
../_images/tutorials_polarized-observations_6_1.png
[4]:
plan = maria.Plan(
    start_time="2024-08-06T03:00:00",
    scan_pattern="daisy",
    scan_options={"radius": 0.5, "speed": 0.1},  # in degrees
    duration=900,  # in seconds
    sample_rate=50,  # in Hz
    scan_center=(0, -23),
    frame="ra_dec")

print(plan)
plan.plot()
Plan:
  start_time: 2024-08-06 03:00:00.000 +00:00
  duration: 900 s
  sample_rate: 50 Hz
  center:
    ra: 00ʰ00ᵐ0.00ˢ
    dec: -23°00’0.00”
  scan_pattern: daisy
  scan_radius: 0.9989°
  scan_kwargs: {'radius': 0.5, 'speed': 0.1}
../_images/tutorials_polarized-observations_7_1.png
[5]:
sim = maria.Simulation(
    instrument,
    plan=plan,
    site="llano_de_chajnantor",
    atmosphere="2d",
    atmosphere_kwargs={"weather": {"pwv": 0.5}},
    map=input_map
)

print(sim)
Constructing atmosphere: 100%|██████████| 10/10 [00:33<00:00,  3.37s/it]
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√s  49.13 uK√s   1.458’
│    1  f150  150 GHz  30 GHz  0.5  12.25 aW√s  60 uK√s    104 uK√s  0.8748’
├ Site:
│   region: chajnantor
│   location: 23°01’45.84”S 67°45’17.28”W
│   altitude: 5.064 km
│   seasonal: True
│   diurnal: True
├ Plan:
│   start_time: 2024-08-06 03:00:00.000 +00:00
│   duration: 900 s
│   sample_rate: 50 Hz
│   center:
│     ra: 00ʰ00ᵐ0.00ˢ
│     dec: -23°00’0.00”
│   scan_pattern: daisy
│   scan_radius: 0.9989°
│   scan_kwargs: {'radius': 0.5, 'speed': 0.1}
├ Atmosphere(10 processes with 10 layers):
│ ├ spectrum:
│ │   region: chajnantor
│ └ weather:
│     region: chajnantor
│     altitude: 5.064 km
│     time: Aug 5 23:07:29 -04:00
│     pwv[mean, rms]: (0.5 mm, 15 um)
└ ProjectedMap:
    shape(stokes, nu, y, x): (4, 1, 685, 685)
    stokes: IQUV
    nu: [90.] GHz
    t: naive
    z: naive
    quantity: rayleigh_jeans_temperature
    units: K_RJ
      min: -1.000e-02
      max: 2.540e-01
    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, psi): (91.73 marcsec, 91.73 marcsec, 0°)
    memory: 30.03 MB
[6]:
tod = sim.run()

print(tod)
tod.plot()
Generating turbulence: 100%|██████████| 10/10 [00:01<00:00,  7.81it/s]
Sampling turbulence: 100%|██████████| 10/10 [00:08<00:00,  1.23it/s]
Computing atmospheric emission: 100%|██████████| 2/2 [00:05<00:00,  2.78s/it, band=f150]
Sampling map: 100%|██████████| 2/2 [00:10<00:00,  5.49s/it, channel=[ 0. inf] Hz]
Generating noise: 100%|██████████| 2/2 [00:01<00:00,  1.78it/s, band=f150]
TOD(shape=(652, 45000), fields=['atmosphere', 'map', 'noise'], units='pW', start=2024-08-06 03:14:59.979 +00:00, duration=900.0s, sample_rate=50.0Hz, metadata={'atmosphere': True, 'sim_time': <Arrow [2025-05-30T18:53:24.796257+00:00]>, 'altitude': 5064.0, 'region': 'chajnantor', 'pwv': 0.5, 'base_temperature': 272.523})
../_images/tutorials_polarized-observations_9_2.png
[7]:
from maria.mappers import BinMapper

mapper = BinMapper(
    center=(0, -23),
    stokes="IQU",
    frame="ra_dec",
    width=1.0,
    height=1.0,
    resolution=1.0 / 256,
    tod_preprocessing={
        "window": {"name": "tukey", "kwargs": {"alpha": 0.1}},
        "remove_spline": {"knot_spacing": 30, "remove_el_gradient": True},
        "remove_modes": {"modes_to_remove": [0]},
    },
    map_postprocessing={
        "gaussian_filter": {"sigma": 1},
        # "median_filter": {"size": 1},
    },
    units="mK_RJ",
)

mapper.add_tods(tod)

output_map = mapper.run()

print(output_map)
Mapping band f090: 100%|██████████| 3/3 [00:00<00:00,  7.48it/s, band=f090, stokes=U]
2025-05-30 18:53:36.035 INFO: Ran mapper for band f090 in 6.912 s.
Mapping band f150: 100%|██████████| 3/3 [00:00<00:00,  7.37it/s, band=f150, stokes=U]
2025-05-30 18:53:42.937 INFO: Ran mapper for band f150 in 6.898 s.
ProjectedMap:
  shape(stokes, nu, y, x): (3, 2, 256, 256)
  stokes: IQU
  nu: [ 90. 150.] GHz
  t: naive
  z: naive
  quantity: rayleigh_jeans_temperature
  units: mK_RJ
    min: -1.184e+02
    max: 1.202e+02
  center:
    ra: 00ʰ00ᵐ0.00ˢ
    dec: -23°00’0.00”
  size(y, x): (1°, 1°)
  resolution(y, x): (14.06”, 14.06”)
  beam(maj, min, psi): (14.06”, 14.06”, 0°)
  memory: 6.291 MB

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

[8]:
output_map.plot(stokes=["I", "Q", "U"], nu_index=[[0], [1]])
../_images/tutorials_polarized-observations_12_0.png
[ ]: