.. _h-maps:

:math:`h`-maps
======================

.. contents:: Table of Contents
   :depth: 2
   :local:

Overview
--------

The ``h_maps`` module provides tools to generate **complex harmonic maps**
:math:`h_{n,m}` from LiteBIRD observations. These maps allow to perform a spin-based mapmaking, they encode the scanning
strategy information. The spin-based mapmaking can be performed with the `SMARTIES <https://github.com/simonsobs/smarties/tree/main>`_ package, it has an interface with LBS :math:`h` maps through the functions :func:`read_lbs_h_maps` and :func:`build_mask_from_lbs_h_maps`.

This implementation follows the formalism introduced in McCallum et al (2021)
(`arXiv:2109.05038 <https://arxiv.org/abs/2109.05038>`_), extended to the HWP case in Takase et al (2024) (`arXiv:2408.03040 <https://arxiv.org/abs/2408.03040>`_).

Mathematical formalism
----------------------
For a scanning strategy with HWP modulation, the harmonic map of spin
order :math:`n,m` at pixel :math:`p` is defined as:

.. math::
   :label: eq-h-maps

   h_{n,m}(p) = \frac{1}{N_p} \sum_{t \in p} e^{i n \psi_t + i m \chi_t}

where:

- :math:`N_p` is the number of observations (hits) in pixel :math:`p`
- :math:`\psi_t` is the parralactic angle of the detector at time sample :math:`t`
- :math:`\chi_t` is the HWP rotation angle at time :math:`t`
- The sum runs over all time samples :math:`t` that fall within pixel :math:`p`

These maps capture how the detector orientation is distributed across the sky
during the scanning strategy. Setting :math:`m=0` is equivalent to the definition of h maps in McCallum et al (2021) without HWP modulation.


The maps can be generated by using the :func:`make_h_maps` function with a list of observations, or by using the observations of a simulation through the method :func:`.make_h_maps` of :class:`simulation`

Example
-----------

.. code-block:: python

   import litebird_sim as lbs
   import numpy as np
   from astropy.time import Time

   # load an IMO
   imo = lbs.Imo()
   # Create a simple simulation
   sim = lbs.Simulation(
    base_path="./output_sim",
    start_time=0,
    duration_s=24 * 3600.0,
    random_seed=12345,
    imo=imo,
   )

   # ... (set up detectors, scanning strategy, etc.) ...

   #Create an observation

   sim.create_observations(
    detectors=dets,
    n_blocks_det=1,
    n_blocks_time=1,
    tod_dtype=precision,
    split_list_over_processes=False,
    allocate_tod=False,    #We only need the pointing information to compute h maps, so we can save memory by not allocating the full TOD
   )
   sim.prepare_pointings()

   result = lbs.make_h_maps(
       observations=sim.observations,
       nside=64,
       n_m_couples=np.array([[0, 0], [2, 0], [4, 0]]),
       output_coordinate_system=lbs.CoordinateSystem.Galactic,
       save_to_file=True,
       output_directory="./h_n_maps",
   )

The ``n_m_couples`` parameter
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The spin orders to compute are specified as a 2D array of ``(n, m)`` pairs::

   # Compute h_{0,0}, h_{2,0} and h_{4,0}
   n_m_couples = np.array([
       [0, 0],
       [2, 0],
       [4, 0],
   ])

Note: for the couple (0,0) the resulting map is not the one expected from the definition of eq :eq:`eq-h-maps`. This definition yield :math:`h_{0,0} = \frac{1}{N_p} \sum_{t \in p} e^{i 0 \psi_t + i 0 \chi_t} = 1` for all observed pixels, which is not very useful. Instead the :math:`h_{0,0}` is an hitmap, i.e. :math:`h_{0,0}(p) = N_p`.

Output: ``HMapsResult``
-----------------------

:func:`make_h_maps` returns a :class:`HMapsResult` object, which contains:

- ``h_maps``: a dictionary indexed by detector name, then by ``(n, m)``
  tuple, each entry being a :class:`.h_map_Re_and_Im` object.
- ``coordinate_system``: the output coordinate system
  (a :class:`.CoordinateSystem` object).
- ``detector_split``: the detector split used.
- ``time_split``: the time split used.
- ``duration_s``: the duration of the observation in seconds.
- ``sampling_rate_Hz``: the sampling rate of the observation in Hz.

Accessing individual maps
~~~~~~~~~~~~~~~~~~~~~~~~~

Each entry in ``h_maps`` is a :class:`.h_map_Re_and_Im` object with:

- ``real``: the real part of :math:`h_{n,m}` (NumPy array of shape ``(Npix,)``)
- ``imag``: the imaginary part (NumPy array of shape ``(Npix,)``)
- ``norm``: property returning :math:`\sqrt{\mathrm{Re}^2 + \mathrm{Im}^2}`,
  with unobserved pixels set to :attr:`.UNSEEN_PIXEL_VALUE`
- ``n``, ``m``: the spin orders
- ``det_info``: the detector name

Example::

   h_2_0 = result.h_maps["detector_name"][2, 0]
   print(h_2_0.real)   # real part
   print(h_2_0.imag)   # imaginary part
   print(h_2_0.norm)   # amplitude

Saving and loading maps
-----------------------

Maps are saved in **HDF5** format, one file per detector:

.. code-block:: text

   output_directory/
   └── h_maps_det_{detector_name}.h5
       ├── attrs: coordinate_system, det, detector_split,
       │          time_split, duration, sampling_rate
       ├── 0,0/
       │   ├── Re
       │   └── Im
       ├── 2,0/
       │   ├── Re
       │   └── Im
       └── ...

To reload maps from disk, use :func:`.load_h_maps_from_file`::

   from litebird_sim.mapmaking.h_maps import load_h_maps_from_file

   result = load_h_maps_from_file("./h_n_maps/h_maps_det_mydetector.h5")

Detector and time splits
------------------------

The ``detector_split`` and ``time_split`` parameters follow the same
conventions as the rest of the litebird_sim mapmaking framework (see
:ref:`mapmaking`). Use ``"full"`` (default) to include all detectors
and all time samples.

MPI support
-----------
The current implementation of :func:`make_h_maps` only allows to distribute observation by detector, i.e. each MPI process computes the h maps for a subset of detectors, but using all time samples.
Example:

.. code-block:: python

   sim.create_observations(
      detectors=dets,
      n_blocks_det=2, # The detectors can be split over several MPI tasks
      n_blocks_time=1,  # But the time samples cannot be split, all are used for each detector
      tod_dtype=precision,
      split_list_over_processes=False,
      allocate_tod=False,
   )
   sim.prepare_pointings()

   result = lbs.make_h_maps(
         observations=sim.observations,
         nside=64,
         n_m_couples=np.array([[0, 0], [2, 0], [4, 0]]),
         output_coordinate_system=lbs.CoordinateSystem.Galactic,
         save_to_file=True,
         output_directory="./h_n_maps",
      )


API reference
-------------

.. autoclass:: litebird_sim.mapmaking.h_maps.h_map_Re_and_Im
   :members: norm
   :undoc-members:

.. autoclass:: litebird_sim.mapmaking.h_maps.HMapsResult
   :undoc-members:

.. autofunction:: litebird_sim.mapmaking.h_maps.make_h_maps

.. autofunction:: litebird_sim.mapmaking.h_maps.save_h_maps

.. autofunction:: litebird_sim.mapmaking.h_maps.load_h_maps_from_file