GAPCYNSKI1969_LunarGravityModel.m

%rendering of gravity model, Gapcynski et al., 1969
%coefficients taken from https://arc.aiaa.org/doi/10.2514/3.5479

clear;
clc;
close all;

resolution =75;

GM = 4902.58;
a = 1738.09;

r = a+0.00001;

lambda = deg2rad(linspace(-180,180,resolution));
phi = deg2rad(linspace(-90,90,resolution));

[LAMBDA, PHI] = meshgrid(lambda, phi);

degree = 7;

C = zeros(degree,degree+1);
S = zeros(degree,degree+1);

%
C(2,1) = -2.192E-4;       %C20
C(2,2) = -0.9019E-5;      %C21
C(2,3) = 0.3293E-4;       %C22
%
C(3,1) = -0.1953E-4;      %C30
C(3,2) = 0.2718E-4;       %C31
C(3,3) = 0.6932E-5;       %C32
C(3,4) = 0.2583E-5;       %C33
%
C(4,1) = -0.8989E-5;      %C40
C(4,2) = 0.8125E-5;       %C41
C(4,3) = -0.2822E-5;      %C42
C(4,4) = -0.1828E-6;      %C43
C(4,5) = 0.6060E-7;       %C44
%
C(5,1) = 0.3220E-5;       %C50
C(5,2) = -0.3245E-5;      %C51
C(5,3) = 0.8523E-6;       %C52
C(5,4) = -0.4586E-6;      %C53
C(5,5) = -0.2826E-7;      %C54
C(5,6) = -0.1936E-9;      %C55
%
C(6,1) = 0.1506E-4;       %C60
C(6,2) = 0.1551E-4;       %C61
C(6,3) = -0.1652E-5;      %C62
C(6,4) = -0.2460E-6;      %C63
C(6,5) = -0.5457E-7;      %C64
C(6,6) = 0.2515E-8;       %C65
C(6,7) = -0.3282E-10;     %C66
%
C(7,1) = 0.5640E-4;       %C70
C(7,2) = 0.5677E-5;       %C71
C(7,3) = -0.5746E-6;      %C72
C(7,4) = -0.6442E-7;      %C73
C(7,5) = -0.1880E-7;      %C74
C(7,6) = -0.5815E-9;      %C75
C(7,7) = -0.5716E-9;      %C76
C(7,8) = -0.1045E-9;      %C77
%
S(2,2) = -0.5564E-5;      %S21
S(2,3) = -0.3910E-5;      %S22
%
S(3,2) = 0.4071E-5;       %S31
S(3,3) = 0.4413E-5;       %S32
S(3,4) = 0.6349E-6;       %S33
%
S(4,2) = 0.1862E-5;       %S41
S(4,3) = -0.1717E-5;      %S42
S(4,4) = -0.1571E-5;      %S43
S(4,5) = -0.2726E-6;      %S44
%
S(5,2) = -0.1005E-4;      %S51
S(5,3) = -0.6792E-6;      %S52
S(5,4) = 0.8224E-6;       %S53
S(5,5) = 0.1271E-7;       %S54
S(5,6) = 0.3457E-9;       %S55
%
S(6,2) = -0.5958E-5;      %S61
S(6,3) = -0.6825E-6;      %S62
S(6,4) = -0.1144E-6;      %S63
S(6,5) = -0.1924E-7;      %S64
S(6,6) = -0.2040E-7;      %S65
S(6,7) = -0.9527E-9;      %S66
%
S(7,2) = -0.7902E-5;      %S71
S(7,3) = -0.1048E-5;      %S72
S(7,4) = 0.3588E-7;       %S73
S(7,5) = 0.6937E-8;       %S74
S(7,6) = 0.4537E-8;       %S75
S(7,7) = 0.8406E-9;       %S76
S(7,8) = 0.1266E-9;       %S77

Pnm = 0;

U = zeros(length(phi),length(lambda));

for ii = 1:length(phi)
  for jj = 1:length(lambda)
    for n = 2:degree
      Pnm = legendre(n,sin(phi(ii)));
      for m = 0:(n-1)
        U(ii,jj) += ((a/r)^n)*Pnm(m+1)*...
        ((C(n,m+1)*cos(m*(lambda(jj)+pi)))+(S(n,m+1)*sin(m*(lambda(jj)+pi))));
      endfor
    endfor
    U(ii,jj) = (GM/r)*(1+U(ii,jj));
  endfor
endfor

[X,Y,Z] = sph2cart(LAMBDA,PHI,U);
surf(X,Y,Z,sqrt(X.^2+Y.^2+Z.^2));
colormap jet;
shading interp;
colorbar;
set(gca,'DataAspectRatio',[1 1 1])

figure(2)
[Cplot,h_plot] = contourf(rad2deg(LAMBDA),rad2deg(PHI),U,15);
colormap jet;
colorbar;
clabel(Cplot,h_plot);
grid on;
grid minor;
pbaspect([2 1 1]);
xlabel("Lunar Longitude");
ylabel("Lunar Latitude");
xticks([-180, -150, -120, -90, -60, -30, 0, 30, 60, 90, 120, 150, 180]);
yticks([-90, -60, -30, 0, 30, 60, 90]);