.. _beamconvolution: Convolve Alms with a Beam to fill a TOD ======================================= .. contents:: Table of Contents :depth: 2 :local: The framework provides the function :func:`.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 :cite:`2010:prezeau:conviqt` and :cite:`2019:duivenvoorden:beamconv`. To populate an existing TOD with a signal, use the function :func:`.add_convolved_sky_to_observations` as demonstrated in the following example: .. testcode:: 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}") .. testoutput:: 113.796 116.186 118.567 120.935 123.300 125.655 128.000 130.335 132.661 134.975 Input Data Format ----------------- The input sky alms and beam alms must be encapsulated in instances of the :class:`SphericalHarmonics` class. These can be stored in a dictionary, using either the channel or detector name as a keyword, or passed directly to :func:`.add_convolved_sky_to_observations`. The routines in :ref:`input_sky` already provide correctly formatted inputs. See below for a desctription of the dataclass :class:`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 :ref:`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 :class:`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 :class:`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 :class:`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 :ref:`maps_and_harmonics`. Example usage: .. testcode:: 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}") .. testoutput:: -0.20471 0.47894 -0.51944 -0.55573 -1.29622 0.27499 0.22891 0.47699 3.24894 Elliptical Gaussian Beam Spherical Harmonics -------------------------------------------- The framework provides a function :func:`.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 :class:`SphericalHarmonics` object with the intensity and (if requested) polarized components of the beam in harmonic space. The function :func:`.generate_gauss_beam_alms` provides a convenient way to compute Gaussian beam harmonics for all detectors in an :class:`.Observation`. It wraps around :func:`.gauss_beam_to_alm` and automatically pulls relevant beam parameters from the Observation object. The parameters are: - `observation`: An :class:`.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 :func:`.generate_gauss_beam_alms` returns a dictionary mapping each detector name to its corresponding :class:`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 :func:`.ticra_grid_to_alm` and :func:`.ticra_cut_to_alm`, which convert a grid/cut file into a :class:`.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: 1. The GRASP file is loaded in memory, using either a :class:`.BeamCut` or a :class:`.BeamGrid` instance. 2. The electric field E is converted into Stokes I/Q/U parameters and saved in a :class:`.BeamStokesPolar` object, which uses (θ, φ) coordinates to store the directions on the 4π sphere where the Stokes parameters have been calculated. 3. The :class:`.BeamStokesPolar` object is converted into a Healpix map via the class :class:`.BeamHealpixMap`, which employs a simple binning procedure. 4. The spherical harmonic coefficients are computed from the Healpix map using Ducc. Methods of the Simulation class ------------------------------- The :class:`.Simulation` class offers various functions to streamline convolution: - :func:`.Simulation.get_gauss_beam_alms`: Generates Gaussian beam alms for all detectors using :class:`.DetectorInfo`. - :func:`.Simulation.get_sky`: Produces sky alms based on an instance of :class:`.SkyGenerationParams`. - :func:`.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: .. testcode:: 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 ------------- .. automodule:: litebird_sim.beam_convolution :members: :undoc-members: :show-inheritance: .. automodule:: litebird_sim.beam_synthesis :members: .. automodule:: litebird_sim.grasp2alm :members: