Custom map simulations¶

In this tutorial we will build a simulation from scratch.

We start by defining a Band that will determine our array’s sensitivity to different spectra. We then generate an array by specifying a field of view, which will be populated by evenly-spaced beams of the given band.

[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)

[2]:

array = {"field_of_view": 0.1,
         "shape": "circle",
         "beam_spacing": 1.5,
         "primary_size": 100,
         "bands": [f090, f150]}

instrument = maria.get_instrument(array=array)

print(instrument)
instrument.plot()

Instrument(1 array)
├ arrays:
│             n           FOV baseline        bands
│  array1  1356  5.953 arcmin      0 m  [f090,f150]
│
└ bands:
          center   width    η         NEP   NET_RJ     NET_CMB          FWHM
   f090   90 GHz  20 GHz  0.5  5.445 aW√s  40 uK√s  49.13 uK√s  8.748 arcsec
   f150  150 GHz  30 GHz  0.5  12.25 aW√s  60 uK√s    104 uK√s  5.249 arcsec

../_images/tutorials_custom-map-simulations_3_1.png

As something to observe, we can download a map and construct a map. We also define a plan to do a daisy scan centered on the center of the map.

[3]:

from maria.io import fetch

map_filename = fetch("maps/tarantula_nebula.h5")

input_map = maria.map.load(
    filename=map_filename,
    nu=150e9,
    width=0.25,
    center=(291.156, -31.23),
    units="uJy/pixel")

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

Downloading https://github.com/thomaswmorris/maria-data/raw/master/maps/tarantula_nebula.h5: 100%|██████████| 514k/514k [00:00<00:00, 110MB/s]

ProjectedMap:
  shape[stokes, nu, t, y, x]: (1, 1, 1, 251, 251)
  stokes: ['I']
  nu: [150.] GHz
  t: [1.74370329e+09] s
  quantity: spectral_flux_density_per_pixel
  units: uJy/pixel
    min: 2.573e+00
    max: 6.945e+02
  center:
    ra: 19ʰ24ᵐ37.44ˢ
    dec: -31°13’48.00”
  size[y, x]: (15 arcmin, 15 arcmin)
  resolution[y, x]: (3.586 arcsec, 3.586 arcsec)
  memory: 0.504 MB

../_images/tutorials_custom-map-simulations_5_2.png

[4]:

site = maria.get_site("llano_de_chajnantor", altitude=5065)

print(site)
site.plot()

Site:
  region: chajnantor
  location: 23°01’45.84”S 67°45’17.28”W
  altitude: 5.065 km
  seasonal: True
  diurnal: True

Downloading https://github.com/thomaswmorris/maria-data/raw/master/world_heightmap.h5: 100%|██████████| 7.34M/7.34M [00:00<00:00, 199MB/s]

../_images/tutorials_custom-map-simulations_6_2.png

[5]:

plan = maria.Plan(
    start_time="2024-08-06T03:00:00",
    scan_pattern="daisy",
    scan_options={"radius": 0.08, "speed": 0.01},  # in degrees
    duration=1200,  # in seconds
    sample_rate=50,  # in Hz
    scan_center=(291.156, -31.23),
    frame="ra_dec")

print(plan)
plan.plot()

Plan:
  start_time: 2024-08-06 03:00:00.000 +00:00
  duration: 1200 s
  sample_rate: 50 Hz
  center:
    ra: 19ʰ24ᵐ37.44ˢ
    dec: -31°13’48.00”
  scan_pattern: daisy
  scan_radius: 9.59 arcmin
  scan_kwargs: {'radius': 0.08, 'speed': 0.01}

../_images/tutorials_custom-map-simulations_7_1.png

[6]:

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

print(sim)

Downloading https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/spectra/am/v2/chajnantor.h5: 100%|██████████| 14.3M/14.3M [00:00<00:00, 85.3MB/s]
Downloading https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/weather/era5/chajnantor.h5: 100%|██████████| 8.00M/8.00M [00:00<00:00, 166MB/s]
Constructing atmosphere: 100%|██████████| 6/6 [00:31<00:00,  5.27s/it]
Generating CMB (nside=1024):   0%|          | 0/1 [00:00<?, ?it/s]
Downloading https://github.com/thomaswmorris/maria-data/raw/master/cmb/spectra/planck.csv: 293kB [00:00, 82.6MB/s]
Generating CMB (nside=1024): 100%|██████████| 1/1 [00:06<00:00,  6.38s/it]

Simulation
├ Instrument(1 array)
│ ├ arrays:
│ │             n           FOV baseline        bands
│ │  array1  1356  5.953 arcmin      0 m  [f090,f150]
│ │
│ └ bands:
│           center   width    η         NEP   NET_RJ     NET_CMB          FWHM
│    f090   90 GHz  20 GHz  0.5  5.445 aW√s  40 uK√s  49.13 uK√s  8.748 arcsec
│    f150  150 GHz  30 GHz  0.5  12.25 aW√s  60 uK√s    104 uK√s  5.249 arcsec
├ 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: 1200 s
│   sample_rate: 50 Hz
│   center:
│     ra: 19ʰ24ᵐ37.44ˢ
│     dec: -31°13’48.00”
│   scan_pattern: daisy
│   scan_radius: 9.59 arcmin
│   scan_kwargs: {'radius': 0.08, 'speed': 0.01}
├ Atmosphere(6 processes with 6 layers):
│ ├ spectrum:
│ │   region: chajnantor
│ └ weather:
│     region: chajnantor
│     altitude: 5.064 km
│     time: Aug 5 23:09:59 -04:00
│     pwv[mean, rms]: (0.5 mm, 15 um)
├ CMB:
│   nside: 1024
│   stokes: ['I', 'Q', 'U']
│   nu: [148.] GHz
│   t: [1.74533873e+09] s
│   quantity: brightness_temperature
│   units: uK_b
│   resolution: 3.435 arcmin
└ ProjectedMap:
    shape[stokes, nu, t, y, x]: (1, 1, 1, 251, 251)
    stokes: ['I']
    nu: [150.] GHz
    t: [1.74370329e+09] s
    quantity: rayleigh_jeans_temperature
    units: K_RJ
      min: 1.232e-05
      max: 3.325e-03
    center:
      ra: 19ʰ24ᵐ37.44ˢ
      dec: -31°13’48.00”
    size[y, x]: (15 arcmin, 15 arcmin)
    resolution[y, x]: (3.586 arcsec, 3.586 arcsec)
    memory: 0.504 MB

[7]:

tod = sim.run()

print(tod)
tod.plot()

Generating turbulence: 100%|██████████| 6/6 [00:00<00:00, 28.72it/s]
Sampling turbulence: 100%|██████████| 6/6 [00:09<00:00,  1.63s/it]
Computing atmospheric emission: 100%|██████████| 2/2 [00:15<00:00,  7.66s/it, band=f150]
Sampling CMB: 100%|██████████| 2/2 [00:59<00:00, 29.50s/it, band=f150]
Sampling map: 100%|██████████| 2/2 [00:23<00:00, 11.73s/it, band=f150]
Generating noise: 100%|██████████| 2/2 [00:05<00:00,  2.99s/it, band=f150]

TOD(shape=(1356, 60000), fields=['atmosphere', 'cmb', 'map', 'noise'], units='pW', start=2024-08-06 03:19:59.978 +00:00, duration=1200.0s, sample_rate=50.0Hz, metadata={'atmosphere': True, 'sim_time': <Arrow [2025-04-22T16:20:53.751089+00:00]>, 'altitude': 5064.0, 'region': 'chajnantor', 'pwv': 0.5, 'base_temperature': 272.519})

../_images/tutorials_custom-map-simulations_9_2.png

[8]:

from maria.mappers import BinMapper

mapper = BinMapper(
    center=(291.156, -31.23),
    frame="ra_dec",
    width=0.25,
    height=0.25,
    resolution=0.25 / 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)
output_map.plot()

2025-04-22 16:21:58.147 INFO: Ran mapper for band f090 in 41.99 s.
2025-04-22 16:22:39.176 INFO: Ran mapper for band f150 in 41.03 s.

ProjectedMap:
  shape[stokes, nu, t, y, x]: (1, 2, 1, 256, 256)
  stokes: ['I']
  nu: [ 90. 150.] GHz
  t: [1.74533896e+09] s
  quantity: rayleigh_jeans_temperature
  units: mK_RJ
    min: -8.865e-01
    max: 1.765e+00
  center:
    ra: 19ʰ24ᵐ37.44ˢ
    dec: -31°13’48.00”
  size[y, x]: (15 arcmin, 15 arcmin)
  resolution[y, x]: (3.516 arcsec, 3.516 arcsec)
  memory: 1.049 MB

../_images/tutorials_custom-map-simulations_10_2.png