Convolve Alms with a Beam to fill a TOD#

The framework provides the function add_convolved_sky(), which performs harmonic-space convolution of sky and beam alms, filling the detector timestreams. It supports both cases, with and without the HWP. The relevant mathematical details are found in [8] and [4].

To populate an existing TOD with a signal, use the function add_convolved_sky_to_observations() as demonstrated in the following example:

import litebird_sim as lbs
import numpy as np
import healpy as hp

start_time_s = 0
time_span_s = 1

lmax = 64
mmax = lmax - 4

# Create a simulation
sim = lbs.Simulation(
   base_path="./output",
   start_time=start_time_s,
   duration_s=time_span_s,
   random_seed=12345,
)

# Define the scanning strategy
sim.set_scanning_strategy(
   lbs.SpinningScanningStrategy(
       spin_sun_angle_rad=0.785_398_163_397_448_3,
       precession_rate_hz=8.664_850_513_998_931e-05,
       spin_rate_hz=0.000_833_333_333_333_333_4,
       start_time=start_time_s,
   ),
   delta_time_s=1,
)

sim.set_instrument(
   lbs.InstrumentInfo(
       boresight_rotangle_rad=0.0,
       spin_boresight_angle_rad=0.872_664_625_997_164_8,
       spin_rotangle_rad=3.141_592_653_589_793,
   ),
)

# Create a detector object
det = lbs.DetectorInfo(
    name="Detector",
    sampling_rate_hz=10.0,
    quat=[0.0, 0.0, 0.0, 1.0],
)

# Initialize the observation
(obs,) = sim.create_observations(detectors=[det])

# Prepare the quaternions used to compute the pointings
sim.prepare_pointings()

ncoeff = lbs.SphericalHarmonics.num_of_alm_from_lmax(lmax,mmax)

np.random.seed(1234)
blms = lbs.SphericalHarmonics(
    values=(np.random.normal(0,1,3*ncoeff)+1j*np.random.normal(0,1,3*ncoeff)).reshape(3,ncoeff),
    lmax=lmax,
    mmax=mmax,
)

alms = lbs.SphericalHarmonics(
    values=hp.synalm(np.ones((4,lmax+1)),lmax=lmax, mmax=lmax,),
    lmax=lmax,
    mmax=lmax,
    coordinates=lbs.CoordinateSystem.Galactic,
)

Convparams = lbs.BeamConvolutionParameters(
    lmax = lmax,
    mmax = mmax,
    single_precision = True,
    epsilon = 1e-5,
)

# Here scan the map and fill tod
lbs.add_convolved_sky_to_observations(
    obs,
    alms,
    blms,
    convolution_params=Convparams,
    pointings_dtype=np.float32,
)

for i in range(obs.n_samples):
   value = np.round(obs.tod[0][i], 3)
   print(f"{value:.3f}")

Input Data Format #

The input sky alms and beam alms must be encapsulated in instances of the SphericalHarmonics class. These can be stored in a dictionary, using either the channel or detector name as a keyword, or passed directly to add_convolved_sky_to_observations(). The routines in Synthetic Sky Maps already provide correctly formatted inputs. See below for a desctription of the dataclass SphericalHarmonics.

Pointing Information #

Pointing data can either be embedded in the observation or passed explicitly via the pointings parameter. If both observations and pointings are provided, they must be consistent—either as a single observation and a single numpy array or as lists of equal length.

The HWP effect can be incorporated via pointing data (see Scanning strategy) or by using the hwp argument. The polarization angles of detectors are derived from the observation attributes.

The option nside_centering=NSIDE shifts the detector pointings to the centers of the corresponding HEALPix pixels at the given NSIDE resolution. This option is useful for debugging and for reducing sub-pixel effects.

Convolution Parameters #

The convolution parameters must be specified via the BeamConvolutionParameters dataclass. Allowed parameters:

lmax (int): Maximum ℓ value for sky and beam coefficients.
mmax (int): Maximum m (azimuthal moment) for beam coefficients, constrained to mmax ≤ lmax - 4.
single_precision (bool): Set to False for 64-bit floating-point calculations. Default: True.
epsilon (float): Desired relative accuracy of interpolation. Default: 1e-5.
strict_typing (bool): If True (default), a TypeError is raised if pointing types do not
match single_precision. If False, the code silently converts types at the expense of memory.

If convolution parameters are omitted, defaults are inferred from sky and beam alms, triggering a warning.

Note

Type Coherence and Memory Efficiency

To optimize performance and minimize memory footprint, the convolution engine (based on ducc0) requires specific floating-point precision. If the data types of the input sky_alms and beam_alms (instances of SphericalHarmonics) do not match the precision specified in convolution_params (e.g., complex64 for single_precision=True), the code will perform a type conversion using astype(..., copy=False).

This means that if the types are already compatible, no copy is made, and the data is processed in-place. However, if a conversion is necessary, a new array might be created or the existing one casted, but the framework is designed to avoid unnecessary deep copies whenever possible.

Container for Spherical Harmonics #

The SphericalHarmonics class stores spherical harmonic coefficients (a_{ℓm}) together with their associated metadata, such as lmax, mmax, and the number of Stokes parameters (nstokes).

In libraries like HealPy, the a_{ℓm} coefficients are stored in plain NumPy arrays. However, this representation is ambiguous: the values of lmax and mmax cannot be uniquely inferred from the array size unless one assumes lmax == mmax.

This class provides an explicit and robust interface to spherical harmonic data, allowing mmax ≤ lmax and supporting polarized (3-Stokes) and intensity-only (1-Stokes) cases. Internally, the coefficients are stored in a NumPy array of shape (nstokes, ncoeff) — this shape is always enforced, even for scalar (intensity-only) data.

It supports: - Shape validation against lmax and mmax - Safe algebraic operations (addition, scaling, convolution) - Resizing via zero-padding or truncation - I/O through Healpy-compatible FITS files

Full documentation can be found in Maps and Spherical Harmonics.

Example usage:

import litebird_sim as lbs
import numpy as np

lmax = 10
mmax = 3
ncoeff = lbs.SphericalHarmonics.num_of_alm_from_lmax(lmax,mmax)

np.random.seed(12345)
coeff = np.random.normal(0,1,ncoeff)+1j*np.random.normal(0,1,ncoeff)

alms = lbs.SphericalHarmonics(values=coeff,lmax=lmax, mmax=mmax)

alms_resized = alms.resize_alm(lmax_out=3,mmax_out=2)

alms_resized_real_part = np.real(alms_resized.values[0])

for r in alms_resized_real_part:
    value = np.round(r,5)
    print(f"{value:.5f}")

Elliptical Gaussian Beam Spherical Harmonics #

The framework provides a function gauss_beam_to_alm(), which analytically computes the spherical harmonic coefficients (a_ℓm) representing a 2D elliptical Gaussian beam, with optional polarization and cross-polar leakage. The parameters are:

lmax: Maximum spherical harmonic degree.
mmax: Maximum harmonic order.
fwhm_rad: Full width at half maximum of the beam, defined as fwhm = sqrt(fwhm_max*fwhm_min) (in radians).
ellipticity: Ellipticity of the beam defined as fwhm_max/fwhm_min (1.0 for circular).
psi_ell_rad: Orientation of the beam’s major axis with respect to the x-axis (radians).
psi_pol_rad: Polarization reference angle with respect to the x-axis (radians). If None, only intensity is computed.
cross_polar_leakage: Cross-polarization leakage factor.

The function returns a SphericalHarmonics object with the intensity and (if requested) polarized components of the beam in harmonic space.

The function generate_gauss_beam_alms() provides a convenient way to compute Gaussian beam harmonics for all detectors in an Observation. It wraps around gauss_beam_to_alm() and automatically pulls relevant beam parameters from the Observation object. The parameters are:

observation: An Observation object containing per-detector syntetic beam properties.
lmax: Maximum spherical harmonic degree.
mmax: Maximum harmonic order. Defaults to lmax.
store_in_observation: If True, the result is stored in the observation.blms attribute. Default False

The function generate_gauss_beam_alms() returns a dictionary mapping each detector name to its corresponding SphericalHarmonics object. This is simple example of usage:

blms = lbs.generate_gauss_beam_alms(
    observation=my_obs,
    lmax=512,
    store_in_observation=True
)

Loading TICRA GRASP files #

LBS incorporates a patched version of the grasp2alm library by Yusuke Takase, which is a Python port of the Fortran code found in Planck Level-S. This library lets to load .grd and .cut files created using TICRA GRASP, a simulation tool for reflector antennas, and convert them into spherical harmonics coefficients that can be passed to the 4π beam convolution module. The two core functionalities are provided by ticra_grid_to_alm() and ticra_cut_to_alm(), which convert a grid/cut file into a SphericalHarmonics object. Both functions require a file-like object as input:

with open("my_grasp_file.grd", "rt") as file_obj:
    harmonics = lbs.ticra_grid_to_alm(file_obj)

These functions rely on the following procedure:

The GRASP file is loaded in memory, using either a BeamCut or a BeamGrid instance.
The electric field E is converted into Stokes I/Q/U parameters and saved in a BeamStokesPolar object, which uses (θ, φ) coordinates to store the directions on the 4π sphere where the Stokes parameters have been calculated.
The BeamStokesPolar object is converted into a Healpix map via the class BeamHealpixMap, which employs a simple binning procedure.
The spherical harmonic coefficients are computed from the Healpix map using Ducc.

Methods of the Simulation class #

The Simulation class offers various functions to streamline convolution:

Simulation.get_gauss_beam_alms(): Generates Gaussian beam alms for all detectors using DetectorInfo.
Simulation.get_sky(): Produces sky alms based on an instance of SkyGenerationParams.
Simulation.convolve_sky(): Convolves sky and beam alms for all observations in the simulation.

These methods are MPI-compatible, distributing inputs based on the job’s detector configuration without requiring broadcast operations.

For a single-task execution, refer to the following example:

import litebird_sim as lbs
import numpy as np

start_time_s = 0
time_span_s = 10

nside = 256

lmax = 2 * nside
mmax = lmax - 4

# Create a simulation
sim = lbs.Simulation(
   base_path="./output",
   start_time=start_time_s,
   duration_s=time_span_s,
   random_seed=12345,
)

# Define the scanning strategy
sim.set_scanning_strategy(
   lbs.SpinningScanningStrategy(
       spin_sun_angle_rad=0.785_398_163_397_448_3,
       precession_rate_hz=8.664_850_513_998_931e-05,
       spin_rate_hz=0.000_833_333_333_333_333_4,
       start_time=start_time_s,
   ),
   delta_time_s=1,
)

sim.set_instrument(
   lbs.InstrumentInfo(
       boresight_rotangle_rad=0.0,
       spin_boresight_angle_rad=0.872_664_625_997_164_8,
       spin_rotangle_rad=3.141_592_653_589_793,
   ),
)

# Create a detector object
det = lbs.DetectorInfo(
    name="Detector",
    sampling_rate_hz=10.0,
    bandcenter_ghz=100.0,
    quat=[0.0, 0.0, 0.0, 1.0],
    fwhm_arcmin=30.0,
    ellipticity=1.0,
    psi_rad=0.0,
    pol_angle_rad=0.0,
)

# Initialize the observation
sim.create_observations(detectors=[det])

# Prepare the quaternions used to compute the pointings
sim.prepare_pointings()

# Gaussian beam alms
blms = sim.get_gauss_beam_alms(lmax=lmax)

# Create the alms to convolve
sky_params = lbs.SkyGenerationParams(
    make_cmb=True,
    make_fg=False,
    nside=nside,
    units="K_CMB",
    apply_beam=False,
    bandpass_integration=False,
    output_type="alm",
    lmax=lmax,
    seed_cmb=12345,
)

alms = sim.get_sky(parameters = sky_params)

Convparams = lbs.BeamConvolutionParameters(
    lmax = lmax,
    mmax = mmax,
    single_precision = True,
    epsilon = 1e-5,
)

sim.convolve_sky(sky_alms=alms,
                 beam_alms=blms,
                 convolution_params=Convparams,
                 pointings_dtype=np.float32,
                 nthreads = 0)

API reference #

class litebird_sim.beam_convolution.BeamConvolutionParameters(lmax: int, mmax: int, single_precision: bool = True, epsilon: float = 1e-05, strict_typing: bool = True)#

Bases: object

Parameters used by the 4π beam convolution code

Fields:

lmax (int): Maximum value for ℓ for the sky and beam coefficients
mmax (int): Maximum value for m (azimuthal moment) for beam coefficients, must be ≦ lmax - 4
single_precision (bool): Set it to False to use 64-bit floating points in the calculation
epsilon (float): The desired relative accuracy of the interpolation
strict_typing (bool): If True (the default), a TypeError` exception will be raised if the type of the pointings does not match ``single_precision. If False, the type of the pointings will be converted silently to make the code works at the expense of additional memory consumption.

epsilon: float = 1e-05#

lmax: int#

mmax: int#

single_precision: bool = True#

strict_typing: bool = True#

litebird_sim.beam_convolution.add_convolved_sky(tod, pointings, sky_alms: ~litebird_sim.maps_and_harmonics.SphericalHarmonics | dict[str, ~litebird_sim.maps_and_harmonics.SphericalHarmonics], beam_alms: ~litebird_sim.maps_and_harmonics.SphericalHarmonics | dict[str, ~litebird_sim.maps_and_harmonics.SphericalHarmonics], hwp_angle: ~numpy.ndarray | None = None, mueller_hwp: ~numpy.ndarray | None = None, input_sky_names: str | None = None, input_beam_names: str | None = None, convolution_params: ~litebird_sim.beam_convolution.BeamConvolutionParameters | None = None, pointings_dtype=<class 'numpy.float64'>, nside_centering: int | None = None, nthreads: int = 0)#

Convolve a set of sky maps with detector beams and add the resulting signals to the time-ordered data (TOD) for multiple detectors.

Parameters:

tod (np.ndarray) – Time-ordered data (TOD) for multiple detectors, to which the convolved signal is added.
pointings (np.ndarray or callable) – Pointing information for each detector. If an array, it should have shape (n_detectors, n_samples, 3). If a callable, it should return pointing data when passed a detector index.
sky_alms (SphericalHarmonics | Dict[str, SphericalHarmonics]) – Spherical harmonic coefficients representing the sky maps. If a dictionary, keys should correspond to detector or channel names.
beam_alms (SphericalHarmonics | Dict[str, SphericalHarmonics]) – Spherical harmonic coefficients representing the beam functions. If a dictionary, keys should correspond to detector or channel names.
hwp_angle (np.ndarray or None, default=None) – Half-wave plate (HWP) angle values for each detector. If None, the classic 4π convolution is used.
mueller_hwp (np.ndarray or None, default=None) – Mueller matrices of the HWP. If None, the classic 4π convolution is used.
input_sky_names (str or None, default=None) – Names of the sky maps to use for each detector. If None, all detectors use the same sky.
input_beam_names (str or None, default=None) – Names of the beam maps to use for each detector. If None, all detectors use the same beam.
convolution_params (BeamConvolutionParameters, optional) – Parameters controlling the convolution, such as resolution and precision. If None, reasonable defaults are chosen based on the sky and beam properties, but a warning will be raised.
pointings_dtype (dtype, optional) – Data type for pointings generated on the fly. If the pointing is passed or already precomputed this parameter is ineffective. Default is np.float64.
nside_centering (int, default=None) – If set, shifts the detector pointings to the centers of the corresponding HEALPix pixels at the given NSIDE resolution. If None, no centering is applied.
nthreads (int, default=0) – Number of threads to use for convolution and in case for HEALPix operations. If set to 0, all available CPU cores will be used.

Raises:

ValueError – If dictionary keys in sky_alms or beam_alms do not match detector names.
AssertionError – If tod and pointings shapes are inconsistent.

Notes

This function loops over all detectors and applies add_convolved_sky_to_one_detector to each one.
If input_sky_names and input_beam_names are provided, they must match the structure of the tod array.
The function modifies tod in place by adding the convolved signals for all detectors.

litebird_sim.beam_convolution.add_convolved_sky_to_observations(observations: ~litebird_sim.observations.Observation | list[~litebird_sim.observations.Observation], sky_alms: ~litebird_sim.maps_and_harmonics.SphericalHarmonics | dict[str, ~litebird_sim.maps_and_harmonics.SphericalHarmonics] | None, beam_alms: ~litebird_sim.maps_and_harmonics.SphericalHarmonics | dict[str, ~litebird_sim.maps_and_harmonics.SphericalHarmonics] | None, pointings: ~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy._typing._array_like._ScalarT]] | list[~numpy.ndarray[tuple[~typing.Any, ...], ~numpy.dtype[~numpy._typing._array_like._ScalarT]]] | None = None, hwp: ~litebird_sim.hwp.HWP | None = None, convolution_params: ~litebird_sim.beam_convolution.BeamConvolutionParameters | None = None, component: str = 'tod', pointings_dtype=<class 'numpy.float64'>, nside_centering: int | None = None, nthreads: int | None = None)#

Applies beam convolution to sky maps and adds the resulting signal to the TOD of one or more observations.

This function processes one or multiple Observation objects, handling sky and beam convolution for all detectors within each observation. It supports Galactic coordinate transformations and HWP effects.

Parameters:

observations (Observation or list of Observation) – A single Observation object or a list of them, containing detector names, pointings, and TOD data.
sky_alms (SphericalHarmonics or dict) – The spherical harmonics representation of the sky signal, either as a single object or a dictionary keyed by detector/channel names.
beam_alms (SphericalHarmonics or dict) – The spherical harmonics representation of the detector beams, either as a single object or a dictionary keyed by detector/channel names.
pointings (np.ndarray, list of np.ndarray, or None, default=None) – Detector pointing matrices. If None, the function extracts pointings from the Observation objects. Pointing information is always assumed to be in Ecliptic coordinates
hwp (HWP | None, default=None) – Half-Wave Plate (HWP) parameters. If None, the function either assumes the information stored in the Observation objects, or, if they are absent, assumes no HWP.
convolution_params (BeamConvolutionParameters | None, default=None) – Parameters controlling the beam convolution, including resolution limits and numerical precision.
component (str, default="tod") – The name of the TOD component to which the computed data is added.
pointings_dtype (dtype, optional) – Data type for pointings generated on the fly. If the pointing is passed or already precomputed this parameter is ineffective. Default is np.float64.
nside_centering (int, default=None) – If set, shifts the detector pointings to the centers of the corresponding HEALPix pixels at the given NSIDE resolution. If None, no centering is applied.
nthreads (int, default=None) – Number of threads to use in the convolution. If None, the function reads from the OMP_NUM_THREADS environment variable.

Notes

If pointings is not provided, it is inferred from the Observation objects. If the pointing is passed
or already precomputed this parameter is ineffective. Default is np.float32.
The function determines the correct sky and beam harmonics from the detector names or channels.
Calls add_convolved_sky to process the TOD for all detectors.

litebird_sim.beam_convolution.add_convolved_sky_to_one_detector(tod_det, sky_alms_det: SphericalHarmonics, beam_alms_det: SphericalHarmonics, pointings_det, mueller_matrix, hwp_angle, convolution_params: BeamConvolutionParameters | None = None, nthreads: int = 0)#

Convolve given sky alms with a detector beam alms and add the result to the TOD of a single detector.

Parameters:

tod_det (np.ndarray) – Time-ordered data (TOD) for the given detector, to which the convolved signal is added.
sky_alms_det (SphericalHarmonics) – Spherical harmonic coefficients representing the sky map for this detector.
beam_alms_det (SphericalHarmonics) – Spherical harmonic coefficients representing the beam function for this detector.
pointings_det (np.ndarray) – Pointing information for the given detector.
mueller_matrix (np.ndarray or None) – Mueller matrix of the HWP for the given detector. If None, the classic 4π convolution is used.
hwp_angle (np.ndarray or None) – Half-wave plate (HWP) angle values for the given detector. If None, the classic 4π convolution is used.
convolution_params (BeamConvolutionParameters, optional) – Parameters controlling the convolution, such as resolution and precision. If None, reasonable defaults are chosen based on the sky and beam properties.
nthreads (int, default=0) – Number of threads to use for convolution. If set to 0, all available CPU cores will be used.

Raises:

TypeError – If strict_typing is enabled and the data type of pointings_det does not match the expected precision.

Notes

If no HWP is present, a standard 4π convolution is performed.
If HWP is present, the Mueller matrix is used to properly handle polarization.
Memory Management: To ensure maximum efficiency, the sky and beam coefficients are passed to the underlying engine using copy=False during type casting. If the input coefficients’ precision (e.g., complex128) differs from the requested convolution_params.single_precision, a conversion will occur. Users should ensure types are coherent to avoid even temporary memory overhead.
The function modifies tod_det in place by adding the convolved signal.

litebird_sim.beam_synthesis.gauss_beam_to_alm(lmax: int, mmax: int, fwhm_rad: float, ellipticity: float = 1.0, psi_ell_rad: float = 0.0, psi_pol_rad: float = 0.0, cross_polar_leakage: float = 0.0) → SphericalHarmonics#

Compute spherical harmonics coefficients a_ℓm representing a Gaussian beam.

The implementation is vectorized for performance and based on the physics described in Planck LevelS: https://github.com/zonca/planck-levelS/blob/master/Beam/gaussbeampol_main.f90

Parameters:

lmax (int) – Maximum spherical harmonic degree ℓ_max.
mmax (int) – Maximum spherical harmonic order m_max.
fwhm_rad (float) – Full width at half maximum (FWHM) of the beam in radians. Defined as fwhm = sqrt(fwhm_max*fwhm_min).
ellipticity (float, optional, default=1.0) – Beam ellipticity. Defined as fwhm_max/fwhm_min. Default is 1 (circular beam).
psi_ell_rad (float, optional, default=0.0) – Orientation of the beam’s major axis wrt the x-axis (radians).
psi_pol_rad (float, optional, default=0.0) – Polarization angle of the beam wrt the x-axis. If None, only the intensity (I) is computed (nstokes=1).
cross_polar_leakage (float, optional, default=0.0) – Cross-polar leakage factor (pure number).

Returns:

The spherical harmonics coefficients representing the beam.

Return type:

SphericalHarmonics

litebird_sim.beam_synthesis.gauss_bl(lmax: int, fwhm_rad: float, pol: bool = True) → ndarray#

Compute the Gaussian beam transfer function b_l analytically (Pure NumPy).

This implementation computes the beam window function including the polarization correction factors (Challinor et al 2000, astro-ph/0008228).

It returns components for T, E, B (excluding Stokes V).

Parameters:

lmax (int) – Maximum spherical harmonic degree ℓ_max.
fwhm_rad (float) – Full width at half maximum (FWHM) of the beam in radians.
pol (bool, optional) –
If True, returns an array of shape (3, lmax+1) containing:
- Index 0: Temperature beam b_l^T
- Index 1: E-mode polarization beam b_l^E
- Index 2: B-mode polarization beam b_l^B
If False, returns a 1D array of shape (lmax+1,) (Temperature only).

Returns:

Array containing the beam transfer function b_l.

Return type:

np.ndarray

litebird_sim.beam_synthesis.generate_gauss_beam_alms(observation: Observation, lmax: int, mmax: int | None = None, channels: FreqChannelInfo | list[FreqChannelInfo] | None = None, store_in_observation: bool = False) → dict[str, SphericalHarmonics]#

Generate Gaussian beam spherical harmonics coefficients for each detector.

This function computes the blms for a 2D Gaussian beam, accounting for detector-specific parameters such as beam width (FWHM), ellipticity, and polarization orientation.

Parameters:

observation (Observation) – Observation object containing detector parameters (names, fwhm, etc.).
lmax (int) – Maximum spherical harmonic degree ℓ_max.
mmax (int, optional) – Maximum spherical harmonic order m_max. If None, it defaults to lmax.
channels (FreqChannelInfo or list of FreqChannelInfo, optional) – Frequency channels to use. If None, uses detectors from the observation.
store_in_observation (bool, default=False) – If True, the computed blms will be stored in observation.blms.

Returns:

Dictionary mapping detector names (str) to SphericalHarmonics objects.

Return type:

dict

This file contains a port of the Grasp2alm repository <https://github.com/yusuke-takase/grasp2alm>, created by Yusuke Takase. We decided to incorporate it into LBS so that we can make the API more compatible with the framework, while preserving the original source code repository.

If you are patching a bug in this file, be sure to check if the same bug was present in Yusuke’s repository; if it is, please upload your patches there, so that both codes can benefit from your work.

class litebird_sim.grasp2alm.BeamCut(file_obj: TextIO)#

Class to hold the data from a beam cut file of GRASP.

This class only supports polar spherical cuts in the far field region.

Parameters:

vini (float) – Initial value.
vinc (float) – Increment.
vnum (int) – Number of values in cut.
c (numpy.ndarray) – Constant.
icomp (int) – Polarization control parameter.
icut (int) – Control parameter of cut.
ncomp (int) – Number of field components.
ncut (int) – Number of cuts.
amp (numpy.ndarray | None) – Amplitude, with a shape (2, num_of_phi_cuts, n_theta). The two fields are the complex amplitudes of the E_co and E_cx components of the electric field

to_polar(copol_axis: str = 'x') → BeamStokesPolar#

Converts beam in “cut” format to Stokes parameters on a polar grid. Assumes that cuts are evenly spaced in theta. The value of copol_axis specifies the alignment of the co-polar basis (‘x’ or ‘y’) of the input GRASP file.

Parameters:: copol_axis (str) – The axis of co-polarization. Must be either ‘x’ or ‘y’.
Returns:: The beam in polar coordinates.
Return type:: BeamStokesPolar
Raises:: ValueError – If the beam is not in the expected format.

class litebird_sim.grasp2alm.BeamGrid(file_obj: TextIO)#

Class to hold the data loaded from a TICRA GRASP beam grid file

This class only supports polar spherical grids in the far field region.

Parameters:

nset (int) – Number of field sets or beams (this class only supports one field)
klimit (int) – Specification of limits in a 2D grid.
field_component_type (SphericalFarFieldDecomposition) – Control parameter of field components. Only types 3 (copolar-crosspolar) and 9 (total power and \(\sqrt(\text{RHC}/\text{LHC}\) are supported)
num_of_components (int) – Number of field components (only ncomp==2 is supported)
grid_type (SphericalFarFieldGrid) – Control parameter of field grid type (only igrid==7 is supported)
ix (int) – Centre of set or beam No. \(i\).
iy (int) – Centre of set or beam No. \(i\).
xs (float) – Start x-coordinate of the grid.
ys (float) – Start y-coordinate of the grid.
xe (float) – End x-coordinate of the grid.
ye (float) – End y-coordinate of the grid.
nx (int) – Number of steps for the x-coordinate
ny (int) – Number of steps for the y-coordinate
frequency (float) – Value of the frequency.
frequency_unit (str) – Measurement unit for frequency.
amp (numpy.ndarray | None) – 2D array of complex amplitudes [\(\theta\), \(\phi\)].

to_polar(copol_axis='x') → BeamStokesPolar#

Converts beam in polar grid format into Stokes parameters on a polar grid.

The value of copol specifies the alignment of the co-polar basis (‘x’ or ‘y’) of the input GRASP file.

Parameters:: copol_axis (str) – The copolarization axis. Must be ‘x’ or ‘y’.
Returns:: The beam grid in polar coordinates.
Return type:: BeamStokesPolar
Raises:: ValueError – If the beam is not in the supported GRASP grid format.

class litebird_sim.grasp2alm.BeamHealpixMap(healpix_map)#

Represents a beam map.

Fields:: nside (int): The resolution parameter of the map. map (numpy.ndarray): The beam map data, in RING order.

to_alm(lmax: int, mmax: int, epsilon=1e-08, max_num_of_iterations=20) → ndarray#

Converts the beam map to spherical harmonic coefficients.

Parameters:

lmax (int) – Maximum l value for the spherical harmonic expansion.
mmax (int) – Maximum m value for the spherical harmonic expansion.
epsilon (float) – Precision of the result
max_num_of_iterations (int) – Maximum number of iterations

Returns:

The spherical harmonic coefficients, as a (3, N) array.

Return type:

numpy.ndarray

Raises:

AssertionError – If lmax is greater than 3*``nside``-1 or if mmax is greater than lmax.

class litebird_sim.grasp2alm.BeamStokesPolar(theta_phi_values_rad: ndarray[tuple[Any, ...], dtype[_ScalarT]], polar_basis_flag: bool = False)#

Representation of a beam as a set of Stokes parameters sampled over a sphere using θφ coordinates

This class stores the four Stokes parameters (\(I\), \(Q\), \(U\), and \(V\)) for a beam on a possibly irregular spherical grid (\(\theta\), \(\varphi\)). This type is used as input to the spherical harmonic transform. Internally, this package uses the third Ludwig’s definition for the polarization basis with the co-polar direction (positive \(Q\)) aligned with the y-axis.

Parameters:

theta_phi_values_rad (npt.NDArray) – A 2D matrix of shape \((N, 2)\) containing the values for \(\theta\) (in theta_phi_values_rad[:, 0]) and for \(\varphi\) (in theta_phi_values_rad[:, 1].
stokes (numpy.ndarray) – Array of shape \((4, N)\) containing the four Stokes parameters (\(I\), \(Q\), \(U\), \(V\), in this order).
polar_basis_flag (bool) – If True, the \(Q\) and \(U\) parameters have been rotated so that they are expressed in a polar basis. This is handy if you are going to compute the harmonic coefficients on them. By default, this is False, because typically BeamPolar instances are created out of TICRA GRASP files, which use Ludwig’s third definition of the polarization.

convert_to_polar_basis() → BeamStokesPolar#

Rotates \(Q\) and \(U\) Stokes parameters from the co-cross basis to the polar basis.

The \(Q\) and \(U\) Stokes parameters are usually represented in the co-cross basis, where the co-polar direction is aligned with the y-axis (consistent with Ludwig 3 convention). For the purposes of extracting the spherical harmonic coefficients, it is more useful to represent them in the polar basis. Unlike the original LevelS’s method, this function operates on a copy of the input BeamPolar, so self is not changed.

Returns:: A new instance of BeamPolar with the rotated stokes beam.
Return type:: BeamPolar

find_best_nside_for_resolution(resol_rad: float | None = None) → int#

Given a typical separation between two samples, estimate a good value for NSIDE

This function can be used to estimate which value of NSIDE to use for building a Healpix map with the representation of the beam, given the minimal separation between two samples in the GRASP file.

Note that this function should not be called for GRASP grids sampled using the (θ, φ) grid scheme, as in this case the spacing between points varies a lot between different co-latitude rings.

If the parameter resol_rad is not provided, the code will call the method minimum_angle_between_samples_rad().

minimum_angle_between_samples_rad() → float#

Return the minimum angle between samples in the GRASP file

This information is useful if you plan to project the samples over a Healpix map, as it provides a starting point for figuring out NSIDE.

See find_best_nside_for_resolution().

Returns:: The minimum separation between samples, in radians

to_map(nside: int, nstokes: int = 3, unseen_pixel_value: float = 0.0, nthreads: int | None = None) → BeamHealpixMap#

Convert the BeamPolar to a BeamMap.

Parameters:

nside (int) – The nside parameter for the HEALPix map. If you are sampling the sphere using a roughly uniform coverage (which means you are not using the IGRID=7 theta/phi sampling!), you can estimate this value using minimum_angle_between_samples_rad().
nstokes (int) – Number of Stokes parameters to project. The default is 3, which means that three maps are produced: I, Q, and U.
unseen_pixel_value (float) – Value to fill outside the valid theta range.
nthreads (int) – Number of threads to use in ducc’s Healpix methods. If None, the function reads from the OMP_NUM_THREADS environment variable.

Returns:

A new instance of BeamMap representing the beam map.

Return type:

BeamMap

class litebird_sim.grasp2alm.SphericalFarFieldDecomposition(value)#

class litebird_sim.grasp2alm.SphericalFarFieldGrid(value)#

litebird_sim.grasp2alm.beam_mapmaker(pixidx: ndarray[tuple[Any, ...], dtype[_ScalarT]], values: ndarray[tuple[Any, ...], dtype[_ScalarT]], output_map: ndarray[tuple[Any, ...], dtype[_ScalarT]], hit_map: ndarray[tuple[Any, ...], dtype[_ScalarT]]) → None#: A simple binning map-maker for Stokes parameters

litebird_sim.grasp2alm.ticra_cut_to_alm(**kwargs) → SphericalHarmonics#

Convert a GRASP .cut file to a spherical harmonic coefficients of beam map.

Parameters:

file_obj (file) – File object to read
nside (int) – Resolution parameter for the output beam map.
lmax (int, optional) – The desired lmax parameters for the analysis.
mmax (int, optional) – The desired mmax parameters for the analysis.
unseen_pixel_value (float) – Value to assign to pixels outside the valid theta range.
interp_method (str) – Interpolation method to use. Default is ‘linear’. Supported are ‘linear’, ‘nearest’, ‘slinear’, ‘cubic’, ‘quintic’ and ‘pchip’.
copol_axis (str, optional) – Axis of the co-polarization component. Defaults to ‘x’.
lmax – Maximum value for ℓ
mmax – Maximum value for m
epsilon (float, optional) – Target precision
max_num_of_iterations (int, optional) – Maximum number of iterations

Returns:

Spherical harmonic coefficients of the beam map.

Return type:

alm (numpy.ndarray)

Raises:

ValueError – If the file format is unknown.

litebird_sim.grasp2alm.ticra_grid_to_alm(**kwargs) → SphericalHarmonics#

Convert a GRASP .grd file to a spherical harmonic coefficients of beam map.

Parameters:

file_obj (file) – File object to read
nside (int) – Resolution parameter for the output beam map.
lmax (int, optional) – The desired lmax parameters for the analysis.
mmax (int, optional) – The desired mmax parameters for the analysis.
outOftheta_val (float) – Value to assign to pixels outside the valid theta range.
interp_method (str) – Interpolation method to use. Default is ‘linear’. Supported are ‘linear’, ‘nearest’, ‘slinear’, ‘cubic’, ‘quintic’ and ‘pchip’.
copol_axis (str, optional) – Axis of the co-polarization component. Defaults to ‘x’.
lmax – Maximum value for ℓ
mmax – Maximum value for m
epsilon (float, optional) – Target precision
max_num_of_iterations (int, optional) – Maximum number of iterations

Returns:

Spherical harmonic coefficients of the beam map.

Return type:

alm (numpy.ndarray)

Raises:

ValueError – If the file format is unknown.