<p>The equations are mostly correct, but the equation for calculating the gravitational potential on the lunar surface needs to be modified to:</p>
<p><span class="math-container">$$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 $$</span></p>
<p>Using this equation, you can generate a 360x360 grid map of the gravitational potential distribution.
<a href="https://i.sstatic.net/3KJgax4l.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/3KJgax4l.png" alt="enter image description here" /></a></p>
<p>Next, compute the spherical harmonics coefficients based on orthogonality.</p>
<p>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.</p>

CC BY-SA 4.0
2024-05-16T23:23:30.950Z
Feel free to accept your own answer, the only criteria is that it answers your question or solves your problem, and acceptance flags future readers doing searches that an accepted answer exists.
uhoh
CC BY-SA 4.0
2024-05-14T11:55:38.123Z
2024-05-14T11:55:38.123Z
Picaqiu
<p>I use the lunar gravity model GRGM900C to calculate the gravitational potential across a <span class="math-container">$360{^\circ}*360{^\circ}$</span> grid. This calculation follows the equation:
<span class="math-container">$$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}
$$</span>
where  <span class="math-container">$\overline a_{nm}$</span>,<span class="math-container">$\overline b_{nm}$</span> are the normalized spherical harmonics given in the model,and <span class="math-container">$\overline{\mathstrut P}_{nm}$</span>   is the normalized associated Legendre polynomial:</p>
<p><span class="math-container">$$\overline{P}_{nm}(\sin\varphi)=\sqrt{2(2n+1)\frac{(n-m)!}{(n+m)!}}P_{nm}(\sin\varphi)$$</span>
<span class="math-container">$$
\\\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&amp;\end{cases}
$$</span></p>
<p>Next, based on orthogonality, I use the following equations to compute the spherical harmonics coefficients:
<span class="math-container">$$\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}$$</span>
where:
<span class="math-container">$$d\Omega=\cos(\varphi)\mathrm{d}\varphi\mathrm{d}\lambda $$</span></p>
<p>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.</p>
<p>here is my code:</p>
<pre><code>degreemax=60;
[anm, bnm] = readCoefficients('GRGM900C_SHA.TAB', degreemax, degreemax);
R=1738;%reference radius of the moon,km
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);
deg2rad=pi/180;
    
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);
               
                cnm(n+1,m+1)=cnm(n+1,m+1)+V(lat+91,lon+1)*Rnm*cosd(lat)*deg2rad*deg2rad;
                snm(n+1,m+1)=snm(n+1,m+1)+V(lat+91,lon+1)*Snm*cosd(lat)*deg2rad*deg2rad;
            end
        end
    end
end
cnm=cnm*k;
snm=snm*k;
 
</code></pre>

CC BY-SA 4.0
2024-05-12T04:29:39.647Z
*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?](https://space.stackexchange.com/q/27689/12102) and [How can I verify my reconstructed gravity field of Ceres from spherical harmonics?](https://space.stackexchange.com/q/27779/12102)
uhoh
CC BY-SA 4.0
2024-05-13T17:00:48.130Z
@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.
Picaqiu
CC BY-SA 4.0
2024-05-13T21:49:42.510Z
*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](https://meta.stackexchange.com/q/377004/303080)) 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
<geophysics><gravity><geodesy>
How to compute spherical harmonics coefficients using orthogonality?
Picaqiu
# How to compute spherical harmonics coefficients using orthogonality?

(Asked by: username_0 on 2024-05-11)

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);
R=1738;%reference radius of the moon,km
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);
deg2rad=pi/180;
    
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);
               
                cnm(n+1,m+1)=cnm(n+1,m+1)+V(lat+91,lon+1)*Rnm*cosd(lat)*deg2rad*deg2rad;
                snm(n+1,m+1)=snm(n+1,m+1)+V(lat+91,lon+1)*Snm*cosd(lat)*deg2rad*deg2rad;
            end
        end
    end
end
cnm=cnm*k;
snm=snm*k;
 

```

**Tags:** geophysics, gravity, geodesy

> **Comments:**
> **username_1**: *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?](https://space.stackexchange.com/q/27689/12102) and [How can I verify my reconstructed gravity field of Ceres from spherical harmonics?](https://space.stackexchange.com/q/27779/12102)  
> **username_0**: @username_1 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.  
> **username_1**: *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 [username_1's lemmas](https://meta.stackexchange.com/q/377004/303080)) 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.  

---

## 1 answers

### Answer by username_0 on 2024-05-14. Score: 1 | Accepted answer

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.
[![enter image description here](https://i.sstatic.net/3KJgax4l.png)](https://i.sstatic.net/3KJgax4l.png)

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.

> **Comments:**
> **username_1**: Feel free to accept your own answer, the only criteria is that it answers your question or solves your problem, and acceptance flags future readers doing searches that an accepted answer exists.  

---

