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 A and S the analysis (map2alm) and synthesis (alm2map) operators and a, m and w, the $a_{\ell m}$, map and pixel mask vectors, the Jacobi iterative process reads
$\displaystyle \textbf{a}^{(n)} = \textbf{a}^{(n-1)} + \textbf{A}. \left( \textbf{w}.\textbf{m}- \textbf{S}.\textbf{a}^{(n-1)} \right),$     (10)

with
a(0) = A.w.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 Nside 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 Nside 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_iterativeinto a FITS file
map2alm_spin
analyze spin weighted maps.

Version 3.50, 2018-12-10