I'm trying to understand how to calculate numerically the spherical harmonic expansion's coefficients considering a specific set of points on the surface of a sphere.
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)$$
and that 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$$
So my question is how to calculate these coefficients given a set of point on the sphere? In this thread Integrating over particular grids to obtain Spherical Harmonic coefficients they pointed out that the Gaussian quadrature is one solution but I am struggling to understand how to retrieve the coefficients $\omega$.
$$f_{m}^{l}=A\sum_{\Theta}^{}\sum_{\varphi}^{}f(\Theta,\varphi)Y_{m}^{l}(\Theta,\varphi)\omega(\Theta,\varphi)d\Omega$$
Does anyone have suggestions how to tackle the problem?