anafast

This program performs harmonic analysis of the HEALPix maps up to a user specified maximum spherical harmonic order $\ell_{\mathrm{max}}$. The integrals are computed on the whole sphere, unless the user chooses a provided option to excise from the input map(s) a simple, constant latitude, symmetric cut, and/or apply an arbitray cut read from an external file. Scalar, or scalar and tensor, spherical harmonic coefficients are evaluated from the map(s) if the input provides, respectively, only the temperature, or temperature and polarisation maps. The total operation count scales as ${\cal {O}}(N_{\mathrm{pix}}^{1/2}\ell_{\mathrm{max}}^2 )$.
Anafast reads one (two) file(s) containing the map(s) and produces a file containing the temperature auto- (or cross-) power spectrum $C^{TT}_\ell$ and, if requested, also the polarisation power spectra $C^{EE}_\ell$, $C^{BB}_\ell$, $C^{TE}_\ell$, $C^{TB}_\ell$, $C^{EB}_\ell$ (as well as $C^{ET}_\ell$, $C^{BT}_\ell$, $C^{BE}_\ell$ if two maps are provided). The $a_{\ell m}$ coefficients computed during the execution also can be written to a (two) file(s) if requested.
Anafast executes an approximate, discrete point-set quadrature on a sphere sampled at the HEALPix pixel centers. Spherical harmonic transforms are computed using recurrence relations for Legendre polynomials on co-latitude, $\theta$, and Fast Fourier Transforms on longitude, $\phi$.
Anafast is provided with an option to use precomputed Legendre Polynomials; please note that since version 2.20 this will most likely reduce performance instead of increasing it.
Anafast permits two execution options which allow a significant improvement of accuracy of the approximate quadrature performed by this facility:
$\bullet$ Improved analyses using either the provided ring weights, which correct the quadrature on iso-latitude rings, or pixel-based weights which improve the quadrature on every pixels, and/or
$\bullet$ An iterative scheme using in succession several backward and forward harmonic transforms of the maps.

Location in HEALPix directory tree: src/f90/anafast/anafast.f90 

RECOMMENDATIONS FOR USERS

Execution of anafast requires a user to specify the maximum spherical harmonic order $\ell_{\mathrm{max}}$ up to which the harmonic decomposition of the input maps will be performed. Since there are no formal limits on parameter $\ell_{\mathrm{max}}$ enforced by anafast, the user should make his/her choices judiciously. Hereafter it is convenient to specify $\ell_{\mathrm{max}}$ in terms of the HEALPix map resolution parameter N_side (called nsmax in some other contexts).

If the function to be analysed is strictly band-width limited, or nearly band-width limited (as in the case of a Gaussian beam smoothed signal discretized at a rate of a few pixels per beam area), it is sufficient to run anafast with $\ell_{\mathrm{max}}\approx 2\cdot N_{\mathrm{side}}$, with a very good $C_\ell$ error performance already in the raw (i.e. uncorrected quadrature) harmonic transform mode. If quadrature corrections are still desired in this case, it should be sufficient to use, at no extra cost in execution time, the ring-weighted quadrature scheme. This is the recommended mode of operation of anafast for essentially error and worry free typical applications, e.g. CPU-intensive Monte Carlo studies.

A new set of pixel-based quadrature weights was introduced in HEALPix 3.40. Pre-computed to inforce a (near) ideal integration of the spherical harmonics $Y_{\ell m}$ on the pixelized sphere (ie $\frac{4\pi}{N_{\mathrm{pix}}} \sum_p w(p) Y_{\ell m}(p) = \sqrt{4 \pi} \delta_{\ell 0}\delta_{m 0}$) for $\vert m\vert \le \ell \le 3 N_{\mathrm{side}}$, they can be used to insure that the $a_{\ell m}$ and $C_\ell$ computed by anafast are perfectly accurate (almost to machine precision) without the need for iterations, but only for band-width limited input signal with $\ell_{\mathrm{max}}\le 1.5 N_{\mathrm{side}}$.

If more aggressive attempts are undertaken to extract from a map the spectral coefficients at $\ell > 2\cdot N_{\mathrm{side}}$ (for example, as in a possible case of an attempt to analyse an existing map, which was irreversibly binned at a suboptimal resolution) the following should be kept in mind:

$\bullet$ Spherical harmonics discretized using HEALPix (either sampled at pixel centers, or avaraged over pixel areas) form a linearly independent system up to $\ell_{\mathrm{max}}= 3 \cdot N_{\mathrm{side}}-1$. Hence, the functions which are strictly band-width limited to $\ell_{\mathrm{max}}= 3 \cdot N_{\mathrm{side}}-1$ can be fully spectrally resolved with anafast, albeit with integration errors in the uncorrected quadrature mode, which grow up to $\delta C_\ell \propto \epsilon \cdot C_\ell$, with $\epsilon <0.1$, at the highest values of $\ell$. These integration errors can be efficiently reduced using anafast in the iterative mode. Although this $\ell_{\mathrm{max}}$ range — $2 \cdot N_{\mathrm{side}}< \ell_{\mathrm{max}}< 3 \cdot N_{\mathrm{side}}- 1$ — is easily manageable with anafast used on strictly band-width limited functions, it should be used with caution in basic and automated applications, e.g. Monte Carlo simulations.

$\bullet$ As with any discrete Fourier transform, anafast application to functions which are not band-width limited results with aliasing of power, which can not be remedied. If the particular case of interest may result in such a band-width violation (i.e. there is significant power in the function at $\ell > 3 \cdot N_{\mathrm{side}}-1$), the function should be smoothed before the application of anafast, or discretized and then analysed, on a refined HEALPix grid (with larger N_side).

$\bullet$ REMEMBER: A peculiar property of the sphere, which usually surprises those whose intuition is built on experience with FFTs on a segment, or on a Euclidean

multidimensional domain, is the lack of a regular and uniform point-set at arbitrary resolution, and the resulting non-commutativity of the forward and backward discrete Fourier transforms on nearly-uniform point-sets, e.g. HEALPix. Hence, as in any case of attempting an extreme application of an off-the-shelf software, use caution and understand your problem well before executing anafast under such circumstances!

FORMAT

% anafast [options] [parameter_file]

COMMAND LINE OPTIONS

-d

--double

double precision mode (see Notes on double/single precision modes)

-s

--single

single precision mode (default)

QUALIFIERS

infile =

Defines the input map file. (default= map.fits) If not blank, the filename should never be put between quotes even if it contains symbols such as &, [, ], ?, = which should be typed literally (ie, unprotected). For instance infile = http://site/action?file.fits[2][col FLUX] is just fine.

infile2 =

Defines the 2nd input map file, to be cross-correlated with the first one. The 2 maps should match in resolution ( N_side) and coordinate. (default= `', only the auto-correlation of the first map will be computed)

outfile =

Defines the output file with the power spectrum. If only one input map is provided, outfile will contains its auto-spectra, if 2 maps are provided, outfile will contain their cross-spectra. Note in particular that in the latter case, the C^{T x E}_l power spectrum will be build from the T field of the 1st (possibly polarized) map, and the E field of the second polarized map. (default= cl_out.fits)

simul_type =

Defines which map(s) to analyse, 1=temperature only, 2=temperature AND polarisation. (default= 1)

nlmax =

Defines the maximum $\ell$ value to be used. See the Recommendations for Users. (default= 64)

maskfile =

Defines the FITS file containing the pixel mask(s) or weighting scheme(s) by which the map(s) read from infile will be multiplied before analysis. If the file contains several fields, the first one in which at least one pixel is non-zero will be used. This option can be used to, for instance, apply the WMAP Kp intensity mask to the data (see https://lambda.gsfc.nasa.gov), but it will not handle the WMAP composite mask correctly. Can be used in conjonction with theta_cut_deg. Masked or weighted pixels will be correctly accounted when performing the monopole and dipole regression.
Note: The mask's resolution ( N_side) and ordering can be different from the input map(s) one's, and the mask will be pro/down-graded and reordered to match the map. On the other hand, the mask and maps coordinates will always be presumed to match (ie, no attempt of rotation of the mask will be made). (default= `': no mask, all valid pixels are used)

theta_cut_deg =

Defines the latitude (in degrees) of an optional, straight symmetric cut around the equator. Pixels located within that cut (|b|<theta_cut_deg) are ignored. (default= 0°: no cut)

iter_order =

Defines the maximum order of quadrature iteration to be used. (default=0, no iteration). For details, see the map2alm_iterative routine described in the ”Fortran Subroutines” document. iter_order>0 can not be used together with won=2.

outfile_alms =

Defines the name of the file containing the $a_{\ell m}$ coefficients of the first map which can be written optionally. (default= no entry — $a_{\ell m}$s are not written to a file)

outfile_alms2 =

Defines the name of the file containing the $a_{\ell m}$ coefficients of the second map, if any, which can be written optionally. (default= no entry — $a_{\ell m}$s are not written to a file)

plmfile =

Defines the name for an input file containing precomputed Legendre polynomials $P_{\ell m}$. (default= no entry — anafast executes the recursive evaluation of $P_{\ell m}$s)

w8file =

Defines name for an input file containing ring-based or pixel-based weights (depending on the value of won) in the improved quadrature mode (default= no entry — see ”Default file names and directories” in introduction)

w8filedir =

Gives the directory where the ring weight files are to be found (default= no entry — anafast searches the default directories, see ”Default file names and directories” in introduction)

won =

Set this to 1 if ring-based quadrature weight files are to be used, or to 2 to use pixel-based weight files instead; otherwise set it to 0. (default= 0). won=2 can not be used together with iter_order>0.

regression =

Sets the degree of the regression made on the input map before doing the power spectrum analysis. The regression is a minimal variance fit (assuming a uniform noise) made on valid (unflagged and unmasked) pixels, out of the symmetric cut (if any). In case of cut sky analysis, such a regression reduces the monopole and dipole leakage to higher $\ell$'s.
0 : no regression, does the $a_{\ell m}$ analysis on the raw map
1 : removes the best fit monopole first
2 : removes the best fit monopole and dipole first
default = 0.

DESCRIPTION

Anafast reads one or two binary FITS-files containing a HEALPix map. These files can each contain a temperature map or both temperature and polarisation (Q,U) maps. Anafast analyses the map(s) and makes an output ascii-FITS file containing the angular auto or cross power spectra $C^{TT}_\ell$s (and $C^{EE}_\ell$, $C^{BB}_\ell$, $C^{TE}_\ell$, $C^{TB}_\ell$ and $C^{EB}_\ell$ if specified, as well as $C^{ET}_\ell$, $C^{BT}_\ell$ and $C^{BE}_\ell$ if two maps are provided). Here $C^{TE}_\ell$ is meant as the power spectrum built from the T field of the first (polarized) map, and the E field of the second polarized map, while it is the other way around for $C^{ET}_\ell$. Anafast produces $C_\ell$s up to a specified maximum $\ell$-value (see Recommendations for Users). If requested, the computed $a_{\ell m}$ coefficients can be written to a FITS file. This file can be used in the constrained realisation mode of synfast.

Anafast permits two execution modes that allow to improve the quadrature accuracy: (1) the ring weight corrected quadrature, and (2) the iterative scheme. Using the ring weights does not increase the execution time. The precomputed ring weights to be used for each HEALPix resolution N_side are provided in the $HEALPIX/data directory. The more sophisticated iterative scheme increases the accuracy more effectively than the weighted ring scheme, but its disadvantage is that the time for the analysis increases, 1 iteration takes 3 times as long, 2 iterations 5 times as long on so forth, since each order of iteration requires one more forward and backward transform.

The spherical harmonics evaluation uses a recurrence on associated Legendre polynomials $P_{\ell m}(\theta)$. This recurrence consumed most of the CPU time used by anafast up to version 2.15. We have therefore included an option to load precomputed values for the $P_{\ell m}(\theta)$ from a file generated by the HEALPix facility plmgen. Since the introduction of accelerated spherical harmonic transforms in HEALPix v2.20, this feature is obsolete and should no longer be used.

When dealing with polarized signal maps, the anafast behavior will depend on the value of the POLCCONV FITS keyword (see note on POLCCONV in The HEALPix Primer)

In version 3.50 a bug affecting previous versions of anafast has been fixed. (It occured previously when iter_order>0 was used in conjonction with a maskfile and/or a restrictive theta_cut_deg, see map2alm_iterative for details). The result was correct when the mask (if any) was applied to the map prior to the anafast calling, or when no iteration was requested.

DATASETS

The following datasets are involved in the anafast processing.

Dataset	Description
data/weight_ring_n0xxxx.fits	Files containing ring weights for the anafast improved quadrature mode.

SUPPORT

This section lists those routines and facilities (including those external to the HEALPix distribution) which can assist in the utilisation of anafast.

synfast

This HEALPix facility can generate a map for analysis by anafast.

alteralm

This HEALPix Fortran facility can be used to modify the $a_{\ell m}$ extracted by anafast before they are used as constraints on a synfast run.

plmgen

This HEALPix facility can be used to generate precomputed Legendre polynomials.

EXAMPLE # 1:

anafast

Anafast runs in interactive mode — self-explanatory.

EXAMPLE # 2:

anafast filename

When 'filename' is present, anafast enters the non-interactive mode and parses its inputs from the file 'filename'. This has the following structure: the first entry is a qualifier which announces to the parser which input immediately follows. If this input is omitted in the input file, the parser assumes the default value. If the equality sign is omitted, then the parser ignores the entry. In this way comments may also be included in the file. In this example, the file contains the following qualifiers:
simul_type= 1
nlmax= 64
theta_cut_deg= 0
iter_order= 0
infile= map.fits
outfile= cl_out.fits
regression= 0
Anafast reads the map from map.fits, makes an analysis and produces C^T_ls up to l=64. This powerspectrum is saved in the file cl_out.fits. No galactic cut is excised and no iterations are performed. As regressionis set to 0 (its default value) the map is analyzed as is, without prior best fit removal of the monopole nor the dipole.

Since
infile2
outfile_alms
outfile_alms2
w8file
w8filedir
plmfile
maskfile
were omitted, they take their default values (empty strings). This means that no file for precomputed Legendre polynomials is read, no second map is read, no mask is applied, and anafast does not save the $a_{\ell m}$ values from the analysis.

Also since
won
is not given, it takes its default value 0, which means that quadrature weights are not used.

RELEASE NOTES

■ Initial release (HEALPix 0.90)

■ Optional non-interactive operation. Proper FITS file support. Improved reccurence algorithm for $P_{\ell m}(\theta)$ which can compute to higher $\ell$ values. New functionality: arbitrary order of iterations, precomputed $P_{\ell m}$, dumping of $a_{\ell m}$. (HEALPix 1.00)

■ New functionality: possibility of removing the best fit monopole and dipole. New Parser. Can be linked to FFTW (HEALPix 1.20)

■ New functionality: addition of maskfile (HEALPix 2.0)

■ Bug correction: correct interaction of iterative scheme with masked pixels (HEALPix 2.01)

■ New functionality: cross-correlation of 2 maps; Correction of this documentation: the code expects maskfile and not mask_file (HEALPix 2.1)

■ Bug correction: now correctly supports mask pro/down-grading

■ Support for pixel-based quadrature weights when won=2 (HEALPix 3.40)

■ Correction of a bug in map2alm_iterative (HEALPix 3.50)

MESSAGES

This section describes error messages generated by anafast

Message	Severity	Text
can not allocate memory for array xxx	Fatal	You do not have sufficient system resources to run this facility at the map resolution you required. Try a lower map resolution.