# How to compute spherical harmonics coefficients using orthogonality?

I use the lunar gravity model GRGM900C to calculate the gravitational potential across a $$360{^\circ}*360{^\circ}$$ grid. This calculation follows the equation: $$V(\varphi,\lambda)=\sum_{\mathrm{n=0}}^\infty\sum_{m=0}^\mathrm{n}[\overline{\mathstrut a}_{nm}\overline{\mathstrut R}_{nm}(\varphi,\lambda)+\overline{\mathstrut b}_{nm}\overline{\mathstrut S}_{nm}(\varphi,\lambda)]^{\leftarrow}$$ where $$\overline a_{nm}$$,$$\overline b_{nm}$$ are the normalized spherical harmonics given in the model,and $$\overline{\mathstrut P}_{nm}$$ is the normalized associated Legendre polynomial:

$$\overline{P}_{nm}(\sin\varphi)=\sqrt{2(2n+1)\frac{(n-m)!}{(n+m)!}}P_{nm}(\sin\varphi)$$ $$\\\begin{cases}\overline{R}_{nm}(\varphi,\lambda)=\overline{P}_{nm}(\mathrm{sin}\varphi)\cdot\mathrm{cos}m\lambda \\\overline{S}_{nm}(\varphi,\lambda)=\overline{P}_{nm}(\mathrm{sin}\varphi)\cdot\mathrm{sin}m\lambda&\end{cases}$$

Next, based on orthogonality, I use the following equations to compute the spherical harmonics coefficients: $$\overline{a}_{\mathrm{nm}}=\frac1{4\pi}\iint V(\varphi,\lambda)\cdot\overline{R}_{\mathrm{nm}}(\varphi,\lambda)d\Omega\leftarrow \\\overline{b}_{\mathrm{nm}}=\frac1{4\pi}\iint V(\varphi,\lambda)\cdot\overline{S}_{\mathrm{nm}}(\varphi,\lambda)d\Omega^{\leftarrow}$$ where: $$d\Omega=\cos(\varphi)\mathrm{d}\varphi\mathrm{d}\lambda$$

Most of the harmonics coefficients match those of GRGM900C, but discrepancies are notable in J2 and J4. Specifically, when the order is 0 and 1, the spherical harmonics coefficients are not equal to 0.

here is my code:

degreemax=60;
[anm, bnm] = readCoefficients('GRGM900C_SHA.TAB', degreemax, degreemax);
GM=4902.799967088640;%km3/s2

V0=GM/R;
V=zeros(181,360);
V=V+V0;
for lat=-90:90 %latitude
for lon=0:359 %longitude
for n=1:degreemax
p=associatedlegendre(n,lat);
for m=0:n
rnm=p(n+1,m+1)*cosd(m*lon);
snm=p(n+1,m+1)*sind(m*lon);
term=anm(n+1,m+1)*rnm+bnm(n+1,m+1)*snm;
V(lat+91,lon+1)=V(lat+91,lon+1)+term;
end
end
end
end

%fit a 60*60  lunar gravity model by gravity potential at discrete points
maxn=60;
cnm=zeros(maxn+1,maxn+1);
snm=zeros(maxn+1,maxn+1);
R=1738;
k=1/(4*pi);

for n=0:maxn
for m=0:n
for lat=-90:90 %latitude
p=associatedlegendre(n,lat);
for lon=0:359 %longitude

Rnm=p(n+1,m+1)*cosd(m*lon);
Snm=p(n+1,m+1)*sind(m*lon);

end
end
end
end
cnm=cnm*k;
snm=snm*k;


• Welcome to Stack Exchange! What an excellent first question. There is more than on e convention for coefficients. You might (or might not) find something helpful in Space SE's Ceres gravity from spherical harmonics from Dawn, how to get the coefficients, definitions and potential? and How can I verify my reconstructed gravity field of Ceres from spherical harmonics?
– uhoh
May 12 at 4:29
• @uhoh Thank you so much for the references you shared; they are helpful. I reviewed the normalization and associated Legendre functions and discovered that the problem didn't lie with them, but rather with other issues in the code. Subsequently, I realized that the discrepancies weren't limited to just J2 and J4; they extended to all even-order terms. Upon closer examination, I pinpointed the mistake. When computing the gravitational potential on the spherical surface, I mistakenly included the zeroth-order term, GM/r, which should actually be omitted. May 13 at 17:00
• Excellent! It is always OK to post answers to your own question, the goal is always to generate answers for the benefit of both the question author and for future readers. (see uhoh's lemmas) It would be great if you rewrote that comment as a short answer post, wait a few days, and if there are no further answers go ahead and accept your own answer.
– uhoh
May 13 at 21:49

The equations are mostly correct, but the equation for calculating the gravitational potential on the lunar surface needs to be modified to:

$$V(\varphi,\lambda)=\sum_{\mathrm{n}=2}^\infty\sum_{m=0}^\mathrm{n}[\overline a_{nm} \overline R_{nm}(\varphi,\lambda)+\overline b_{nm} \overline S_{nm}(\varphi,\lambda)]^\leftarrow$$

Using this equation, you can generate a 360x360 grid map of the gravitational potential distribution.

Next, compute the spherical harmonics coefficients based on orthogonality.

To verify the results, use the fitted harmonics to draw the gravitational potential map again. It should match the one generated using the gravity model.

New contributor
Picaqiu is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.