(Sept 12) functions for read/write of lensing potential and deflection (spin) maps, alm rotation, etc.
(Dec 11) functions for transforming arrays of maps and alms, packed alm and array formats made public, HealpixMap_MakeSeparableNGfromG function for simulating local non-Gaussianity. Updates for Healpix version linking.
(April 11) fixed bug in the lensing routines, giving (small near phi=0) error in lensed maps
(January 11) Various compatibility improvements, gfortran changes and bug fixes; distribution now includes all updated WeightMixer power spectrum code and functions.
(October 10) Support for new CAMB unlensed power spectrum format and CLEψ (see arXiv:1101.2234)
(December 07) Minor fixes, lower memory usage (MPI 2), new routines for unlensed spherical transforms or arrays of maps and calculating all cross power spectra
(November 07) Major update including new cubic interpolation simulation method, speedups.
The code requires HEALPix 2.x and allows you to:
Download the code
A version-tracked (git) copy of the code is available on a shared repository site; after downloading please email if you would instead like to have access to the repository and hence latest builds.
There is no HEALPix-quality documentation, though details of the theory are given in the appendix of my original paper and the splining scheme in Appendix E.4 of arXiv:0801.0554. A discussion of the lensing potential-E correlation is in Lewis, Challinor & Hanson (arXiv:1101.2234). The download includes a sample program SimLens.f90 which should give you a good idea of how things work using the HealpixObj wrapper routines. If you prefer lower-level HEALPix-style routines you can also use spin_alm_tools.f90 directly. The code can be run on single CPUs without MPI, though for high resolution fast lensed map making you will likely need more than 2GB of memory. Rather than wasting lots of time documenting things now (which probably no one will read), please ask any questions on CosmoCoffee and I will respond there so everyone can read it..
Note for intel compilers: you may need to compile with -heap-arrays to avoid getting segmentation faults using some routines with large lmax.
To run the demo SimLens program you will need to edit the Makefile and job file, and also change params.ini to have the correct path to your HEALPix data files. Then run with "qsub job" or "./simlens params.ini". The program will read in some unlensed Cls from a text file, simulate unlensed a_lms, compute the gradient of the potential, then use this to construct a lensed map, outputting the text lensed C_ls in that realization.
The way MPI is implemented is that you call a HealpixInit(H,..) function with an nside and lmax you want to use before anything else. All but the main thread stay in this routine until HealpixFree(H) is called by the main thread. The main thread calls transforms like HealpixMap2Alm(H,....), which are automatically farmed out over the other threads. This is demonstrated in the SimLens demo program. There are numerous ways things could be improved.
There are almost certainly bugs, please let me know if you find out what they are, as well as any improvements. The implementation is not complete, in that there are no fits routines for storing spin fields or lensing maps, etc, etc. Using interpolation for the pixel re-mapping might well save a lot of time and memory, and probably worth implementing if you want to do lots of simulations.
For an accurate calculation of the lensing potential, unlensed and lensed CMB power spectra see CAMB. For parameter estimation including CMB lensing see CosmoMC.
The cubic interpolation method works by simulating an equicylindricallly-pixellized high-resolution map, then doing cubic interpolation (ref) on this grid for the deflection angles in each Healpix pixel. The high-res pixels are equally spaced in θ, but the number of pixels in φ depends on the number of MPI threads (for each MPI slice it chooses the number of pixels in φ so that the maximum angular spacing is similar to or smaller than the θ-spacing, and the number of φ-pixels is of the form 2n3m so that the FFT is fast). To this extent the result is not the invariant under changing the number of threads. Using a simpler cubic interpolation scheme would speed up the code significantly, however it appears to make it much less stable - the result does then not converge as you increase interp_factor.
For a given lmax using nside=2048 and interp_factor=1.5 seems to be give power spectra accurate to about 0.1% at l<2000 (or nside=1024 and interp_factor=3; for accurate lensing basically needs an effective high-res nside of about 3000). Errors on BB power spectrum are larger. Note that all spectra are affected by choice of lmax; BB especially requires high lmax for accurate results (and an input power spectrum with non-linear corrections).
The expectation over simulations of the simulated lensed Cl should agree with the lensed power spectrum from CAMB (using accurate full-sky method of astro-ph/0502425). For this to be true either lmax needs to be much larger than the range of interest, so all relevant power is included, or lmax must be the same in both cases (and of course pixellization/interpolation errors must be negligible). For lmax=2250 the smoothed average fractional difference from 50 simulated spectra does indeed agree well (~0.1%) on all spectra (with usual caveats about BB).