Is there any way to integrate spherical harmonics over the surface of a unit sphere in python?

Question

I am trying to integrate this in Python

f1 = lambda phi, theta: ( np.conj(sph_harm(m1, l1, theta, phi)) \
                           * sph_harm(m2, l2, theta, phi)     #theta--> 0-2pi, phi-->0-pi
                  * np.sin(phi) / np.power(1+epsilon*np.cos(n*phi), l2))

coeff1 = (l1*l2*(l2+1)) /(2*l1*(l1+1))
result = coeff1 * term1
print(result)

where m1, l1, m2, l2 have certain values. For example m1,m2=-2, l1,l2=2.
I have also used Mathematica to check the values I get from python. They don't match for different values. For example, for the above mentioned values of m1,l1,m2,l2, from python I obtained 0.01673 while Mathematica gives me 1.001433. For other higher values of l and ms the difference is monumental sometimes.

Then, I have tried Gaussian quadrature method to numerically calculate my result. Though it matches for some values of l and m with the result from Mathematica, It does not hold for all values.

The weird thing is I found the problem occurring only when I consider np.sin(phi) / np.power(1+epsilon*np.cos(n*phi), l2)) in the integration. If I only try to integrate the spherical harmonics, it does not show any difference. Do you have any idea why this is happening and how to fix it?

Answer 1

Note in many cases on potential theory, I was better of with something like

import numpy as np
from scipy.integrate import simps

xl = np.linspace( a, b, A ) 
yl = np.linspace( c, d, C )
xx, yy = np.meshgrid( xl, yl )
zz = xx.copy()
for i in range( len( yy ) ):
    for j in range( len( xx[0] ) ):
        x0 = xx[ i, j ]
        y0 = yy[ i, j ]
        zz[ i, j ] = f( x0, y0 )
result = simps( simps( zz, xl ), yl )

rather than using dblquad.