Numerical computation of the Spherical harmonics expansion

Question

I'm trying to understand the Spherical harmonics expansion in order to solve a more complex problem but the result I'm expecting from a very simple calculation is not correct. I have no clue why this is happening.

A bit of theory: It is well known that a function on the surface of a sphere ( $\mathbb{R}^{3}$ ) can be defined as an infinite sum of some constant coefficients $f_{m}^{l}$ and the spherical harmonics $Y_{m}^{l}(\Theta,\varphi)$ :

$f(\Theta,\varphi)=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}f_{m}^{l}Y_{m}^{l}(\Theta,\varphi)$

The spherical harmonics are defined as :

$Y_{m}^{l}(\Theta,\varphi)=\sqrt{\frac{(2l+1)\cdot(l-m)!}{4\pi\cdot(l+m)!}}P_{m}^{l}(cos\Theta)e^{im\varphi}$

where $P_{m}^{l}$ are the associated Legendre polynomials.

An finally, the constant coefficients $f_{m}^{l}$ can be calculated (similarly to the Fourier transform) as follow:

$f_{m}^{l}=\oint_{\Omega}^{}f(\Theta,\varphi)Y_{m}^{l}(\Theta,\varphi)d\Omega$

The problem: Let's assume we have a sphere centered in $(0,0,0)$ where the function on the surface is equal to $1$ for all points $(\Theta,\varphi)$ . We want to calculate the constant coefficients and then calculate back the surface function by approximation. Since $f(\Theta,\varphi)=1$ the calculation of the constant coefficients reduces to :

$f_{m}^{l}=\sum_{\Omega}Y_{m}^{l}(\Theta,\varphi)$

which numerically (in Python) can be approximated using:

def Ylm(l,m,theta,phi):
    return scipy.special.sph_harm(m,l,theta,phi)

def flm(l,m):
    phi, theta = np.mgrid[0:pi:101j, 0:2*pi:101j]
    return Ylm(l,m,theta,phi).sum()

Then, by computing a band limited sum over $l$ I'm expecting to see $f'(\Theta,\varphi)\rightarrow1$ when $l\rightarrow\infty$ for any given point $(\Theta,\varphi)$ .

$f'(\Theta,\varphi)=\sum_{l=0}^{L}\sum_{m=-l}^{l}f_{m}^{l}Y_{m}^{l}(\Theta,\varphi)$

L = 20
f = 0
theta0, phi0 = 0.0, 0.0
for l in xrange(0,L+1):
    for m in xrange(-l,l+1):
        f += flm(l,m)*Ylm(l,m,theta0,phi0)
print f

but for $L=20$ it gives me $(12860.5407362+0j)$ and not $\sim(1+0j)$ . For $L=40$ it gives me $(28156.0642429+0j)$

I know it seems more a Mathematics problem but the formulas should be correct. The problem seems being on my computation. It could be a really stupid mistake but I cannot spot it. Any suggestion?

Thanks

Numerical computation of the Spherical harmonics expansion

Answers (1)

Related Questions