$\newcommand{\d}{\mathrm d}$First time using StackExchange, so please excuse any mistakes:
My goal is to express functions on the sphere as a sum of zonal harmonics.
I've spent a while with the simpler case, where we have a function of one variable and we want to express it with a sum of Legendre polynomials, and have gotten that to work: For example, if we have the function $f(x) = sin(\pi x)$, we can then use the fact that $$ f(x) = \sum_{l=0}^\infty a_l P_l(x)$$ where $P_l(x)$ is the Legendre polynomial of order $l$ to first find coefficients $$ a_l = \frac{2l+1}{2} \int_{-1}^1 f(x) P_l(x) dx $$
Numerically, my approach to build approximation function $g(x)$ is to first find my coefficients with the equation $$a_l = \frac{2l+1}{2} \int_{0}^{\pi} f(\cos(x))P_l(\cos(x))\sin(x) \d x $$ and then to express $$ g(x) = \sum_{l=0}^\infty a_l P_l(x)$$
As far as I can tell, this has been working fine. My problem is when I now have equations of the form $f(\phi, \theta)$. I'm trying to use the same approach as above, but the details are eluding me: such an attempt as the following yields horribly wrong results: $$a_l = \frac{2l+1}{2} \int_{0}^\pi \int_{0}^{\pi}f(\phi, \theta)P_l(\cos\theta)P_l(\cos\phi)\sin(\theta) \sin(\phi)\d\theta \d\phi$$ $$ g(\phi, \theta) = \sum_{l=0}^\infty a_l P_l(\cos \theta) P_l(\cos \phi)$$
Am I simply making mistakes in my computation, or is there something wrong in this approach for spherical functions? Any help would be much appreciated. Thanks!