Map-making

Harlequin includes a routine that applies the destriping algorithm to time-ordered data, producing a maximum-likelihood map of the intensity and polarization signal.

The destriping algorithm used in the code is described in the paper Destiping CMB temperature and polarization maps, Kurki-Suonio et al., A/A 506, 1511-1539 (2009), and the source code closely follows the terminology introduced in that paper.

The destriping algorithm is effective in removing 1/f noise originated within the detectors of an instrument, provided that the noise is uncorrelated among detectors. It requires the user to specify a baseline, i.e., the maximum time span in which the noise can be assumed to be uncorrelated (i.e., white).

Since the destriper is effective only when much data is available, it is often the case that the input data to be fed to the algorithm is larger than the available memory on a machine. In this case, the destriper can take advantage of a distributed-memory environment and of the MPI libraries, if they have been loaded before Harlequin, like in the following example:

import MPI
import Harlequin  # Ok, use MPI whenever possible

# ...

You can check if MPI is being used with the function use_mpi.

Destriping is based on a datatype, DestripingData, and on the function destripe!. To use the destriper, you create an object of type DestripingData and then call destripe! on it. The TOD must be split in a set of observations, using the datatype Observation; we will see how to use them in the section Splitting a TOD into observations.

Splitting a TOD into observations

Harlequin.Observation — Type.

mutable struct Observation{T <: Real}

An observation is a sequence of time-ordered data (TOD) that has been acquired by some detector. It implements the following fields:

time: an array of N elements, representing the time of the sample. Being defined as an AbstractArray, a range (e.g., 0.0:0.1:3600.0) can be used to save memory
pixidx: an array of N elements, each being a pixel index in a map. In no way the Healpix pixelization scheme is enforced; in fact, the map must not even be spherical.
psi_angle: an array of N elements, each containing the orientation of the polarization angle with respect to some fixed reference system (e.g., the Ecliptic plane). The angles must be in radians.
tod: an array of N elements, containing the actual data measured by the instrument and already calibrated to some physical units (e.g., K_CMB)
sigma_squared: an array of N elements, containing the squared RMS of each sample. To save memory, you should usually use a RunLengthArray here.
flagged: a Boolean array of N elements, telling whether the sample should be discarded (true) or not (false) during the map-making process.
name: a string representing the detector. This is used only for debugging purposes.

To create an observation, the constructor takes keyword arguments instead of parameters; in this way, the code should be more readable. Also, sensible defaults will be provided for the missing fields:

# We do not specify "sigma_squared" nor "flagged", so they will be
# initialized to a vector of ones and a vector of `false`
obs = Observation{Float64}(
    time = 0.0:0.1:3600.0,
    pixidx = 1:3601,
    psi_angle = zeros(3601),
    tod = randn(3601),
)

Map-making

Splitting a TOD into observations

High-level functions

Low-level functions