1
$\begingroup$

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;
 
$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – uhoh
    May 12 at 4:29
  • 1
    $\begingroup$ @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. $\endgroup$
    – Picaqiu
    May 13 at 17:00
  • $\begingroup$ 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. $\endgroup$
    – uhoh
    May 13 at 21:49

1 Answer 1

1
$\begingroup$

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

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.
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – uhoh
    May 16 at 23:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.