HEALpy - Getting spherical harmonics from as function of pixel, l and m

Question

I am in need to calculate an integral of spherical harmonics and their products of different orders (l,m) over a fraction of a sphere (let's say for simplicity over a half a sphere).

Healpix/Healpy package makes it very easy and fast to calculate spherical harmonics transform be replacing an integral of a function times spherical harmonics (SH) with a sum of function*SH values in a predefined pixel (see this link).

However, I do not see how to get explicit values of spherical harmonics from healpy directly (i.e. just a function something like Ylm(index,l,m)).

I tried so far to use a package sphtools which has spherical harmonics function, and translate it into healpix with a simple code:

import numpy as np
    
import pyshtools as pysh
import healpy as hp
    
def Ylm_sphtools(theta,phi,l,m,csphase=-1):
    Ylm_sphtools= np.vectorize(pysh.expand.spharm_lm)(l,m,theta,phi,normalization='ortho',kind='complex',                                     degrees=False,csphase=csphase)
    return Ylm_sphtools

def get_Ylm_healpy(nside,l,m):
    npix=hp.nside2npix(nside)
    theta,phi=hp.pix2ang(nside,np.arange(npix))
    return Ylm_sphtools(theta,phi,l,m)

>get_Ylm_healpy(512,10,5)

does what it needs: gives me the values of Y(l=10,m=5) for every pixel on a healpy map with NSIDE=512, but as one may have understood this is very slow (compared with the speed of healpy.anafast function and so on).

So my question is: how can I get spherical harmonics (real and imag.) pars over a healpix map? The integration from this point is simple: sum the values of your Ylm over a mask of interest and it's done.

Thanks.

question from:https://stackoverflow.com/questions/65844496/healpy-getting-spherical-harmonics-from-as-function-of-pixel-l-and-m

Answer 1

The answer from one of healpix's designers, Erik Hivon.

To get the spherical harmonics Ylm or order l,m in healpix format with nside, one does:

import healpy, numpy
alms  = numpy.zeros(  ((l+1)*(l+2) - (l-m)*(l-m+1))//2 , dtype=complex)
alms[-1] = amplitude
ylm = healpy.alm2map(alms, nside=nside, pol=False, sigma=0.0, lmax=l, mmax=m)

amplitude = 0.5 for the real part of the function, and amplitude = 0.5j for the imaginary.

The result is the difference of real and imaginary part (imaginary part should be multiplied by 1j).

Seems legit, I will test it and report the results soon here.

EDIT: Yes, this works. It reproduces well the map2alm command if one uses not healpy.map2alm, but integrals of spherical harmonics obtained above.