Aerospace Toolbox 2.4

Comparing Zonal Harmonic Gravity Model to Other Gravity Models

Examine the zonal harmonic, spherical and 1984 World Geodetic System (WGS84) gravity models for latitudes from +/- 90 degrees at the surface of the Earth.

Contents

Determine Earth-Centered Earth-Fixed (ECEF) Position

Since ECEF coordinate system is geocentric,you can use spherical equations to determine the position from the latitude, longitude and geocentric radius.

Calculate the geocentric radii in meters at an array of latitudes from +/- 90 degrees using geocradius.

lat = -90:90;
r = geocradius( lat );
rlat = convang( lat, 'deg', 'rad');
z = r.*sin(rlat);

Because longitude has no effect for zonal harmonic gravity model, assume that the y position is zero.

x = r.*cos(rlat);
y = zeros(size(x));

Compute Zonal Harmonic Gravity for Earth

Use gravityzonal to calculate array of zonal harmonic gravity in ECEF coordinates for array of ECEF positions in meters per seconds squared.

[gx_zonal gy_zonal gz_zonal] = gravityzonal([x' y' z']);

Calculate WGS84 Gravity

Use gravitywgs84 to compute WGS84 gravity in down-axis and north-axis at the Earth's surface, an array of geodetic latitudes in degrees and 0 degrees longitude using the exact method with atmosphere, no centrifugal effects, and no precessing.

lat_gd = geoc2geod(lat, r);
long_gd = zeros(size(lat));
[gd_wgs84 gn_wgs84] = gravitywgs84( zeros(size(lat)), lat_gd, long_gd, 'Exac
t', [false true false 0] );

Determine Gravity for Spherical Earth

Compute array of spherical gravity for array of geocentric radii in meters per second squared using Earth's gravitational parameter in meters cubed per second squared.

GM  = 3.986004415e14;
gd_sphere = -GM./(r.*r);

Comparison Plots for Different Gravity Models

To compare the gravity models, their outputs must be in the same coordinate system. You can transform zonal gravity from ECEF coordinates to NED coordinates by using the Direction Cosine Matrix from dcmecef2ned.

dcm_ef = dcmecef2ned( lat_gd, long_gd );
gxyz_zonal = reshape([gx_zonal gy_zonal gz_zonal]', [3 1 181]);
for m = length(lat):-1:1
    gned_zonal(m,:) = (dcm_ef(:,:,m)*gxyz_zonal(:,:,m))';
end

figure(1)
plot(lat_gd, gned_zonal(:,3), lat, -gd_sphere, ':', lat_gd, gd_wgs84, '-.')
legend('Zonal Harmonic', 'Spherical', 'WGS84','Location','North')
xlabel('Geodetic Latitude (degrees)')
ylabel('Down gravity (m/s^2)')
grid

Figure 1: Gravity in the Down-axis in meters per second squared

figure(2)
plot( lat_gd, gned_zonal(:,1), lat_gd, gn_wgs84, 'r-.' )
hold on
plot([lat(1) lat(end)], [0 0], 'g:', 'LineWidth', 1.5 )
legend('Zonal Harmonic', 'WGS84', 'Spherical','Location','Best')
xlabel('Geodetic Latitude (degrees)')
ylabel('North gravity (m/s^2)')
grid
hold off

Figure 2: Gravity in the North-axis in meters per second squared

Calculate total gravity for WGS84 and from zonal gravity vector in meters per second squared.

gtotal_wgs84 = gravitywgs84( zeros(size(lat)), lat_gd, zeros(size(lat)), 'Ex
act', [false true false 0] );
gtotal_zonal = sqrt(sum([gx_zonal.^2 gy_zonal.^2 gz_zonal.^2],2));

figure(3)
plot(lat, gtotal_zonal, lat_gd, gtotal_wgs84)
legend('Zonal Harmonic', 'WGS84','Location','North')
xlabel('Geodetic Latitude (degrees)')
ylabel('Total gravity (m/s^2)')
grid

Figure 3: Total gravity in meters per second squared