Interface with ducc0 noise generation and documentation improvements

CHANGELOG.md

@@ -1,5 +1,9 @@
 # HEAD
 
+-   Ducc0 1/f noise generation module included [#490](https://github.com/litebird/litebird_sim/pull/490)
+
+-   Small optimizations in beam convolver [#492](https://github.com/litebird/litebird_sim/pull/492)
+
 -   Fixed a TypeError in Observation when allocate_tod=False in MPI jobs [#491](https://github.com/litebird/litebird_sim/pull/491)
 
 -   Save memory in pointing generation [#488](https://github.com/litebird/litebird_sim/pull/488)

docs/source/beam_convolution.rst

@@ -80,6 +80,7 @@ as demonstrated in the following example:
         values=hp.synalm(np.ones((4,lmax+1)),lmax=lmax, mmax=lmax,),
         lmax=lmax,
         mmax=lmax,
+        coordinates=lbs.CoordinateSystem.Galactic,
     )
 
     Convparams = lbs.BeamConvolutionParameters(
@@ -193,7 +194,8 @@ It supports:
 - 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::
@@ -402,9 +404,6 @@ API reference
     :undoc-members:
     :show-inheritance:
 
-.. automodule:: litebird_sim.maps_and_harmonics
-    :members:
-
 .. automodule:: litebird_sim.beam_synthesis
     :members:
 

docs/source/index.rst

@@ -19,7 +19,6 @@ This is the User’s Manual of the LiteBIRD Simulation Framework.
    part5.rst
    part6.rst
    appendix.rst
-   installation
 
 Indices and tables
 ==================

docs/source/map_scanning.rst

@@ -210,7 +210,7 @@ transparent:
   sim.prepare_pointings()
 
   values_tuple = (np.ones(12*nside*nside)*1e-4, np.ones(12*nside*nside)*1e-4, np.ones(12*nside*nside)*1e-4)
-  sky_signal = lbs.HealpixMap(values=values_tuple, nside=nside)
+  sky_signal = lbs.HealpixMap(values=values_tuple, nside=nside, coordinates=lbs.CoordinateSystem.Galactic)
 
   sim.fill_tods(sky_signal)
 

docs/source/noise.rst

@@ -8,7 +8,10 @@ function :func:`.add_noise_to_observations` and the low-level versions
 :func:`.add_noise`, :func:`.add_white_noise`, and
 :func:`.add_one_over_f_noise`.
 
-Here is a short example that shows how to add noise:
+White noise
+-----------
+
+Here is a short example that shows how to add white noise to timelines:
 
 .. testcode::
 
@@ -90,7 +93,42 @@ call the low level function directly:
    custom_sigma_uk = 1234
    lbs.noise.add_white_noise(obs[0].tod[0], custom_sigma_uk, random=sim.dets_random[0])
 
-We can also add 1/f noise using a very similar call to the above:
+
+1/f Models and Engines
+----------------------
+
+The framework supports also the generatio of 1/f noise. Here you can choose the computational **engine** 
+(how it is calculated) and the physical **model** (the shape of the power spectrum).
+
+Engines
+^^^^^^^
+
+The engine is selected via the ``engine`` parameter:
+
+1.  **"fft" (Default)**: Generates noise in the Fourier domain.
+2.  **"ducc"**: Uses time-domain infinite impulse response (IIR) filtering provided by the `ducc0` library. It only supports the "keshner" model.
+
+Models
+^^^^^^
+
+The physical shape of the Power Spectral Density (PSD) is selected via the ``model`` parameter:
+
+1.  **"toast" (Default)**:
+    The classic power-law ratio, also implemented in https://github.com/hpc4cmb/toast/blob/372fa7642bbe61a5f01d239e707c04b80ad4bf46/src/toast/tod/sim_noise.py#L74. The PSD is proportional to:
+
+    .. math::
+
+        P(f) \propto \frac{f^\alpha + f_{knee}^\alpha}{f^\alpha + f_{min}^\alpha}
+
+2.  **"keshner"**:
+    Corresponds to a sum of relaxation processes. This is the native model of the `ducc` engine. The PSD is proportional to:
+
+    .. math::
+
+        P(f) \propto \left( \frac{f^2 + f_{knee}^2}{f^2 + f_{min}^2} \right)^{\alpha/2}
+
+
+This call allows to add 1/f noise:
 
 .. testcode::
 

litebird_sim/beam_synthesis.py

@@ -286,6 +286,7 @@ def gauss_bl(lmax: int, fwhm_rad: float, pol: bool = True) -> np.ndarray:
           - 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

litebird_sim/hwp_diff_emiss.py

@@ -68,13 +68,13 @@ def add_2f_to_observations(
     By default, the TOD is added to ``Observation.tod``. If you want to add it to some
     other field of the :class:`.Observation` class, use `component`::
 
-    for cur_obs in sim.observations:
-        # Allocate a new TOD for the 2f alone
-        cur_obs.2f_tod = np.zeros_like(cur_obs.tod)
+        for cur_obs in sim.observations:
+            # Allocate a new TOD for the 2f alone
+            cur_obs.hwp_2f_tod = np.zeros_like(cur_obs.tod)
 
-        # Ask `add_2f_to_observations` to store the 2f
-        # in `observations.2f_tod`
-        add_2f_to_observations(sim.observations, component="2f_tod")
+            # Ask `add_2f_to_observations` to store the 2f
+            # in `observations.hwp_2f_tod`
+            add_2f_to_observations(sim.observations, component="hwp_2f_tod")
     """
     if isinstance(observations, Observation):
         obs_list = [observations]

litebird_sim/hwp_harmonics.py

@@ -486,7 +486,7 @@ def fill_tod(
             ``allocate_tod=False`` in :class:`.Observation`.
 
         maps : HealpixMap | SphericalHarmonics | dict[str, HealpixMap] |
-               dict[str, SphericalHarmonics]
+            dict[str, SphericalHarmonics]
             Sky model to be scanned. In dictionary form, the keys must match the
             entries of `input_names`.
 

litebird_sim/maps_and_harmonics.py

@@ -2319,8 +2319,7 @@ def synthesize_alm(
 
     Optimized behavior based on input correlations:
     - 'TT' only: Scalar generation (fastest).
-    - 'TT', 'EE', 'BB', 'TE': Block-diagonal generation. Computes (T, E)
-      as a correlated 2x2 pair, and B independently.
+    - 'TT', 'EE', 'BB', 'TE': Block-diagonal generation. Computes (T, E) as a correlated 2x2 pair, and B independently.
     - 'TB' or 'EB' present: Full 3x3 covariance generation.
 
     Parameters