map2alm_iterative

This routine covers and extends the functionalities of map2alm: it analyzes a (polarised) HEALPix RING ordered map and returns its $a_{\ell m}$ coefficients for temperature (and polarisation) up to a specified multipole, and use precomputed harmonics if those are provided, but it also can also perform an iterative (Jacobi) determination of the $a_{\ell m}$, and apply a pixel mask if one is provided.
Denoting $\textbf{A}$ and $\textbf{S}$ the analysis (map2alm) and synthesis (alm2map) operators and $\textbf{a}, \textbf{m}$ and $\textbf{w}$, the $a_{\ell m}$, map and pixel mask vectors, the Jacobi iterative process reads

$$\textbf{a}^{(n)} = \textbf{a}^{(n-1)} + \textbf{A}. \left( \textbf{w}.\textbf{m}- \textbf{S}.\textbf{a}^{(n-1)} \right),$$ (10)

with

$$\textbf{a}^{(0)} = \textbf{A}.\textbf{w}.\textbf{m}.$$ (11)

During the processing, the standard deviation of the input map $\left(\textbf{w}.\textbf{m}\right)$ and the current residual map $\left(\textbf{w}.\textbf{m}- \textbf{S}.\textbf{a}^{(n-1)}\right)$ is printed out, with the latter expected to get smaller and smaller as $n$ increases.
The standard deviation of map $x$ has the usual definition $\sigma \equiv \sqrt{\sum_{p=1}^{N}\frac{(x(p)-\bar{x})^2}{N-1}}$, where the mean is $\bar{x} \equiv \sum_{p=1}^{N} \frac{x(p)}{N}$, and the index $p$ runs over all pixels.
In version 3.50 a bug affecting previous versions of map2alm_iterative has been fixed. (It occured when iter_order$> 0$ was used in conjonction with a mask and/or a restrictive zbounds, with a magnitude that depended on each of those factors and was larger for non-boolean masks (ie, $\textbf{w}^2 \ne \textbf{w}$). To assess the impact of this bug on previous results, the old implementation remains available in map2alm_iterative_old). The result was correct when the mask (if any) was applied to the map prior to the map2alm_iterative calling, or when no iteration was requested.

Location in HEALPix directory tree: src/f90/mod/alm_tools.F90 

FORMAT

call map2alm_iterative( nsmax, nlmax, nmmax, iter_order, map_TQU, alm_TGC[, zbounds, w8ring_TQU, plm, mask] )

ARGUMENTS

name & dimensionality	kind	in/out	description
nsmax	I4B	IN	the $N_{\mathrm{side}}$ value of the map to analyse.
nlmax	I4B	IN	the maximum $\ell$ value ( $\ell_{\mathrm{max}}$) for the analysis.
nmmax	I4B	IN	the maximum $m$ value for the analysis.
iter_order	I4B	IN	the order of Jacobi iteration. Increasing that order improves the accuracy of the final $a_{\ell m}$ but increases the computation time $T_{\mathrm{CPU}} \propto 1 + 2 \times$iter_order. iter_order $=0$ is a straight analysis, while iter_order $=3$ is usually a good compromise.
map_TQU(0:12*nsmax**2-1, 1:p)	SP/ DP	INOUT	input map. $p$ is 1 or 3 depending if temperature (T) only or temperature and polarisation (T, Q, U) are to be analysed. It will be altered on output if a mask is provided and/or if iter_order$> 0$ and zbounds is provided.
alm_TGC(1:p, 0:nlmax, 0:nmmax)	SPC/ DPC	OUT	The $a_{\ell m}$ values output from the analysis. $p$ is 1 or 3 depending on whether polarisation is included or not. In the former case, the first index is (1,2,3) corresponding to (T,E,B).
zbounds(1:2), OPTIONAL	DP	IN	section of the map on which to perform the $a_{\ell m}$ analysis, expressed in terms of $z=\sin(\mathrm{latitude}) = \cos(\theta).$ If zbounds(1)$<$zbounds(2), it is performed on the strip zbounds(1)$<z<$zbounds(2); if not, it is performed outside the strip zbounds(2)$\le z \le$zbounds(1). If absent, the whole map is processed.

w8ring_TQU(1:2*nsmax,1:p), OPTIONAL	DP	IN	ring weights for quadrature corrections. p is 1 for a temperature analysis and 3 for (T,Q,U). If absent, the ring weights are all set to 1.
plm(0:,1:p), OPTIONAL	DP	IN	If this optional matrix is passed, precomputed scalar (and tensor) $P_{\ell m}(\theta)$ are used instead of recursion.
mask(0:12*nsmax**2-1,1:q), OPTIONAL	SP/ DP	IN	pixel mask, assumed to have the same resolution (and RING ordering) as the map. The map map_TQU is multiplied by that mask before being analyzed, and will therefore be altered on output. $q$ should be in $\{1,2,3\}$. If $p=q=3$, then each of the 3 masks is applied to the respective map. If $p=3$ and $q=2$, the first mask is applied to the first map, and the second mask to the second (Q) and third (U) map. If $p=3$ and $q=1$, the same mask is applied to the 3 maps. Note: the output $a_{\ell m}$ are computed directly on the masked map, and are not corrected for the loss of power, correlation or leakage created by the mask.

EXAMPLE:

use healpix_types
use alm_tools
use pix_tools
integer(i4b) :: nside, lmax, npix, iter
real(sp), allocatable, dimension(:,:) :: map
real(sp), allocatable, dimension(:) :: mask
complex(spc), allocatable, dimension(:,:,:) :: alm

nside = 256
lmax = 512
iter = 2
npix = nside2npix(nside)
allocate(map(0:npix-1,1:3))
allocate(mask(0:npix-1))
mask(0:) = 0. ! set unvalid pixels to 0
mask(0:10000-1) = 1. ! valid pixels
allocate(alm(1:3, 0:lmax, 0:lmax)
call map2alm_iterative(nside, lmax, lmax, iter, map, alm, mask=mask)

Analyses temperature and polarisation signals in the first 10000 pixels of map (as determined by mask). The map has an $N_{\mathrm{side}}$ of 256, and the analysis is supposed to be performed up to 512 in $\ell$ and $m$. The resulting $a_{\ell m}$ coefficients for temperature and polarisation are returned in alm. Uniform weights are assumed. In order to improve the $a_{\ell m}$ accuracy, 2 Jacobi iterations are performed.

MODULES & ROUTINES

This section lists the modules and routines used by map2alm_iterative.

map2alm

Performs the alm analysis

alm2map

Performs the map synthesis

misc_util

module, containing:

assert_alloc

routine to print error message when an array is not properly allocated

RELATED ROUTINES

This section lists the routines related to map2alm_iterative

anafast

executable using map2alm_iterative to analyse maps.

alm2map

routine performing the inverse transform of map2alm_iterative.

alm2map_spin

synthesize spin weighted maps.

dump_alms

write $a_{\ell m}$ coefficients computed by map2alm_iterative into a FITS file

map2alm_spin

analyze spin weighted maps.