Numerical computation of the Spherical harmonics expansion
<p>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.</p>

<hr>

<p><strong>A bit of theory:</strong> It is well known that a function on the surface of a sphere (<img src="http://latex.codecogs.com/gif.latex?%5Cmathbb%7BR%7D%5E%7B3%7D" alt="\mathbb{R}^{3}">) can be defined as an infinite sum of some constant coefficients <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=f_%7Bm%7D%5E%7Bl%7D" alt=""> and the spherical harmonics <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=Y_%7Bm%7D%5E%7Bl%7D(%5CTheta,%5Cvarphi)" alt=""> : </p>

<p><img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=f(%5CTheta,%5Cvarphi)=%5Csum_%7Bl=0%7D%5E%7B%5Cinfty%7D%5Csum_%7Bm=-l%7D%5E%7Bl%7Df_%7Bm%7D%5E%7Bl%7DY_%7Bm%7D%5E%7Bl%7D(%5CTheta,%5Cvarphi)" alt=""></p>

<p>The spherical harmonics are defined as : </p>

<p><img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=Y_%7Bm%7D%5E%7Bl%7D(%5CTheta,%5Cvarphi)=%5Csqrt%7B%5Cfrac%7B(2l%2b1)%5Ccdot(l-m)!%7D%7B4%5Cpi%5Ccdot(l%2bm)!%7D%7DP_%7Bm%7D%5E%7Bl%7D(cos%5CTheta)e%5E%7Bim%5Cvarphi%7D" alt=""></p>

<p>where <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=P_%7Bm%7D%5E%7Bl%7D" alt=""> are the associated Legendre polynomials.</p>

<p>An finally, the constant coefficients <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=f_%7Bm%7D%5E%7Bl%7D" alt=""> can be calculated (similarly to the Fourier transform) as follow:</p>

<p><img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=f_%7Bm%7D%5E%7Bl%7D=%5Coint_%7B%5COmega%7D%5E%7B%7Df(%5CTheta,%5Cvarphi)Y_%7Bm%7D%5E%7Bl%7D(%5CTheta,%5Cvarphi)d%5COmega" alt=""></p>

<hr>

<p><strong>The problem:</strong> Let's assume we have a sphere centered in <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=(0,0,0)" alt=""> where the function on the surface is equal to <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=1" alt=""> for all points  <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=(%5CTheta,%5Cvarphi)" alt="">. We want to calculate the constant coefficients and then calculate back the surface function by approximation.  Since <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=f(%5CTheta,%5Cvarphi)=1" alt=""> the calculation of the constant coefficients reduces to :</p>

<p><img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=f_%7Bm%7D%5E%7Bl%7D=%5Csum_%7B%5COmega%7DY_%7Bm%7D%5E%7Bl%7D(%5CTheta,%5Cvarphi)" alt=""></p>

<p>which numerically (in Python) can be approximated using:</p>

<pre><code>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()
</code></pre>

<p>Then, by computing a band limited sum over <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=l" alt=""> I'm expecting to see <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=f&#39;(%5CTheta,%5Cvarphi)%5Crightarrow1" alt="">  when <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=l%5Crightarrow%5Cinfty" alt=""> for any given point <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=(%5CTheta,%5Cvarphi)" alt="">.</p>

<p><img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=f&#39;(%5CTheta,%5Cvarphi)=%5Csum_%7Bl=0%7D%5E%7BL%7D%5Csum_%7Bm=-l%7D%5E%7Bl%7Df_%7Bm%7D%5E%7Bl%7DY_%7Bm%7D%5E%7Bl%7D(%5CTheta,%5Cvarphi)" alt=""></p>

<pre><code>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
</code></pre>

<p>but for <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=L=20" alt=""> it gives me <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=(12860.5407362%2b0j)" alt=""> and not <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=%5Csim(1%2b0j)" alt="">. For <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=L=40" alt=""> it gives me <img src="https://chart.googleapis.com/chart?cht=tx&amp;chl=(28156.0642429%2b0j)" alt=""></p>

<p>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?</p>

<p>Thanks </p>

<p>The spherical harmonics are orthonormal with the inner product</p>

<pre><code>&lt;f|g&gt; = Integral( f(theta,phi)*g(theta,phi)*sin(theta)*dphi*dtheta)
</code></pre>

<p>So you should calulate the coefficients by</p>

<pre><code>  clm = Integral( Ylm( theta, phi) * sin(theta)*dphi*dtheta)
</code></pre>

python math