Observation simulation and thermal noise#

Instrument models#

ps_eor.obssimu provides telescope models for evaluating frequency-dependent SEFD and, where station layouts are available, simulated UV coverage. Sensitivity calculations can use a telescope model directly:

import numpy as np

from ps_eor import obssimu

freqs = np.arange(50, 80, 0.2) * 1e6
telescope = obssimu.Telescope.from_name("nenufar")
stokes_i_sefd = telescope.get_i_sefd(freqs)

Available telescope models#

The following names can be passed to from_name():

Name	Telescope model	Default baseline range	Notes
`"dex"`	`DEx`	0.5–25 \(\lambda\)	Configurable square dipole array; drift scan
`"ska_low"`	`SkaLow` (AA4)	30–250 \(\lambda\)	Full SKA-Low Phase 1 layout
`"ska_low_aastar"`	`SkaLowAAstar`	30–250 \(\lambda\)	SKA-Low AA* layout
`"ska_low_aa2"`	`SkaLowAA2`	30–250 \(\lambda\)	SKA-Low AA2 layout
`"lofar_hba"`	`LofarHBA`	50–250 \(\lambda\)	LOFAR HBA core
`"a12_hba"`	`A12HBA`	10–200 \(\lambda\)	AARTFAAC-12 HBA
`"a12_lba"`	`A12LBA`	20–40 \(\lambda\)	AARTFAAC-12 LBA; drift scan
`"mwa1"`	`MWA1`	18–80 \(\lambda\)	MWA Phase I core
`"hera"`	`HERA`	4–200 \(\lambda\)	Configurable redundant hexagonal array; drift scan
`"hera_56"`	`HERA56`	4–200 \(\lambda\)	56-antenna HERA layout; drift scan
`"hera_120"`	`HERA120`	4–200 \(\lambda\)	120-antenna HERA layout; drift scan
`"hera_208"`	`HERA208`	4–200 \(\lambda\)	208-antenna HERA layout; drift scan
`"hera_320"`	`HERA320`	4–200 \(\lambda\)	320-antenna HERA layout; drift scan
`"nenufar"`	`NenuFAR`	6–60 \(\lambda\)	Full NenuFAR layout
`"nenufar_80"`	`NenuFAR80`	6–60 \(\lambda\)	80-mini-array NenuFAR layout
`"ovro_lwa"`	`OVROLWA`	2–60 \(\lambda\)	OVRO-LWA; drift scan

The default baseline limits are used by TelescopeSimu unless umin or umax is specified explicitly.

UV coverage#

TelescopeSimu projects physical baselines through an observation. Configure the pointing, hour-angle interval, and time resolution, then grid the samples:

telescope = obssimu.DEx(n_antenna_side=16, sep_antenna=6)
observation = obssimu.TelescopeSimu(
    telescope,
    freqs,
    dec_deg=-27,
    hal=-1,
    har=1,
    timeres=60,
)
coverage = observation.image_gridding(
    fov_deg=telescope.fov,
    min_weight=1,
)

Frequencies are in Hz, physical coordinates in metres, UV coordinates in wavelengths, hour angles in hours, and times in seconds.

Noise power and sensitivity#

The gridded simulation produces a NoiseStdCube and a matching power-spectrum estimator:

noise_std = coverage.get_noise_std_cube(
    total_time_sec=100 * 3600,
)
ps_gen = coverage.get_ps_gen(
    filter_kpar_min=0.05,
    filter_wedge_theta=0,
)

kbins = np.logspace(np.log10(ps_gen.kmin), np.log10(0.5), 8)
noise_spectra = ps_gen.get_all(kbins, noise_std)

The .data arrays in these products are the expected thermal-noise power. They are not the sensitivity. The corresponding one-sigma sensitivity is in .err; for example, the spherical sensitivity is noise_spectra.ps3d.err.

A random realization is needed only for testing realization-dependent processing:

noise_cube = noise_std.generate_noise_cube()

Foreground-mode filtering#

filter_kpar_min and filter_wedge_theta configure modes excluded from the spherical average:

ps_gen = coverage.get_ps_gen(
    filter_kpar_min=0.05,    # h cMpc^-1
    filter_wedge_theta=15,   # degrees
)

The first removes low-|k_parallel| modes; the second removes modes below the wedge line for the requested angle. They change the sensitivity and effective mode count but do not alter the cylindrical power array or simulate foreground subtraction.

Coherent and incoherent accumulation#

Increase total_time_sec for observations that coherently measure the same sky modes, such as repeated nights of one LST block. The longer coherent integration reduces the visibility-noise amplitude and hence its expected noise power:

noise_std = coverage.get_noise_std_cube(total_time_sec=100 * 3600)

Independent sky fields or LST blocks are combined incoherently at the power-spectrum level. Their average has the same expected thermal-noise power, but its statistical uncertainty falls as 1 / sqrt(n_incoherent). For a returned spherical product this can be represented explicitly as:

n_incoherent = 20
ps3d = ps_gen.get_ps3d(kbins, noise_std)
ps3d.err /= np.sqrt(n_incoherent)
ps3d.w *= n_incoherent
ps3d.n_eff *= n_incoherent

Command-line simulation#

The same workflow is available through pstool. First simulate one UV track, then reuse it for a chosen accumulated time and filtering selection:

$ pstool simu_uv nenufar 8.5 --total_time 10 --int_time 60 \
    --bandwidth 10 --df 195.3 --out_filename coverage
$ pstool simu_noise_ps coverage_nenufar_z8.5_10h.h5 100 \
    --filter_kpar_min 0.05 --filter_wedge_theta 15 \
    --n_incoherent_avg 20 --out_filename sensitivity

simu_noise_img draws an image-domain realization, while simu_noise_ps_zrange evaluates a selected spherical scale across redshift. See Simulation and sensitivity commands for all arguments and output products.