GeographicLib  1.47
Gravity models
Back to Geoid height. Forward to Magnetic models. Up to Contents.

GeographicLib can compute the earth's gravitational field with an earth gravity model using the GravityModel and GravityCircle classes and with the Gravity utility. These models expand the gravitational potential of the earth as sum of spherical harmonics. The models also specify a reference ellipsoid, relative to which geoid heights and gravity disturbances are measured. Underlying all these models is normal gravity which is the exact solution for an idealized rotating ellipsoidal body. This is implemented with the NormalGravity class and some notes on are provided in section Normal gravity

Go to

The supported models are

See

• W. A. Heiskanen and H. Moritz, Physical Geodesy (Freeman, San Francisco, 1967).

Acknowledgment: I would like to thank Mathieu Peyréga for sharing EGM_Geoid_CalculatorClass from his Geo library with me. His implementation was the first I could easily understand and he and I together worked through some of the issues with overflow and underflow the occur while performing high-degree spherical harmonic sums.

# Installing the gravity models

Available gravity models
name max
degree
size
(kB)
tar file Windows
installer
zip file
egm84
180
27
egm96
360
2100
egm2008
2190
76000
wgs84
20
1

The "size" column is the size of the uncompressed data.

For Linux and Unix systems, GeographicLib provides a shell script geographiclib-get-gravity (typically installed in /usr/local/sbin) which automates the process of downloading and installing the gravity models. For example

   geographiclib-get-gravity all  # to install egm84, egm96, egm2008, wgs84
geographiclib-get-gravity -h   # for help


This script should be run as a user with write access to the installation directory, which is typically /usr/local/share/GeographicLib (this can be overridden with the -p flag), and the data will then be placed in the "gravity" subdirectory.

Windows users should download and run the Windows installers. These will prompt for an installation directory with the default being

   C:/ProgramData/GeographicLib


(which you probably should not change) and the data is installed in the "gravity" sub-directory. (The second directory name is an alternate name that Windows 7 uses for the "Application Data" directory.)

Otherwise download either the tar.bz2 file or the zip file (they have the same contents). To unpack these, run, for example

   mkdir -p /usr/local/share/GeographicLib
tar xofjC egm96.tar.bz2 /usr/local/share/GeographicLib
tar xofjC egm2008.tar.bz2 /usr/local/share/GeographicLib
etc.


and, again, the data will be placed in the "gravity" subdirectory.

However you install the gravity models, all the datasets should be installed in the same directory. GravityModel and Gravity uses a compile time default to locate the datasets. This is

• /usr/local/share/GeographicLib/gravity, for non-Windows systems
• C:/ProgramData/GeographicLib/gravity, for Windows systems

consistent with the examples above. This may be overridden at run-time by defining the GEOGRAPHICLIB_GRAVITY_PATH or the GEOGRAPHICLIB_DATA environment variables; see GravityModel::DefaultGravityPath() for details. Finally, the path may be set using the optional second argument to the GravityModel constructor or with the "-d" flag to Gravity. Supplying the "-h" flag to Gravity reports the default path for gravity models for that utility. The "-v" flag causes Gravity to report the full path name of the data file it uses.

# The format of the gravity model files

The constructor for GravityModel reads a file called NAME.egm which specifies various properties for the gravity model. It then opens a binary file NAME.egm.cof to obtain the coefficients of the spherical harmonic sum.

The first line of the .egm file must consist of "EGMF-v" where EGMF stands for "Earth Gravity Model Format" and v is the version number of the format (currently "1").

The rest of the File is read a line at a time. A # character and everything after it are discarded. If the result is just white space it is discarded. The remaining lines are of the form "KEY WHITESPACE VALUE". In general, the KEY and the VALUE are case-sensitive.

GravityModel only pays attention to the following keywords

• keywords that affect the field calculation, namely:
• ReferenceRadius (required), the equatorial radius a for the reference ellipsoid meters.
• ModelMass (required), the mass constant GM for the model in meters3/seconds2.
• ReferenceMass (required), the mass constant GM for the reference ellipsoid in meters3/seconds2.
• AngularVelocity (required), the angular velocity ω for the model and the reference ellipsoid in rad seconds−1.
• Flattening, the flattening of the reference ellipsoid; this can be given a fraction, e.g., 1/298.257223563. One of Flattening and DynamicalFormFactor is required.
• DynamicalFormFactor, the dynamical form factor J2 for the reference ellipsoid. One of Flattening and DynamicalFormFactor is required.
• HeightOffset (default 0), the constant height offset (meters) added to obtain the geoid height.
• CorrectionMultiplier (default 1), the multiplier for the "zeta-to-N" correction terms for the geoid height to convert them to meters.
• Normalization (default full), the normalization used for the associated Legendre functions (full or schmidt).
• ID (required), 8 printable characters which serve as a signature for the .egm.cof file (they must appear as the first 8 bytes of this file).
The parameters ModelRadius, ModelMass, and AngularVelocity apply to the gravity model, while ReferenceRadius, ReferenceMass, AngularVelocity, and either Flattening or DynamicalFormFactor characterize the reference ellipsoid. AngularVelocity (because it can be measured precisely) is the same for the gravity model and the reference ellipsoid. ModelRadius is merely a scaling parameter for the gravity model and there's no requirement that it be close to the equatorial radius of the earth (although that's typically how it is chosen). ModelMass and ReferenceMass need not be the same and, indeed, they are slightly difference for egm2008. As a result the disturbing potential T has a 1/r dependence at large distances.
• keywords that store data that the user can query:
• Name, the name of the model.
• Description, a more descriptive name of the model.
• ReleaseDate, when the model was created.
• keywords that are examined to verify that their values are valid:
• ByteOrder (default little), the order of bytes in the .egm.cof file. Only little endian is supported at present.

Other keywords are ignored.

The coefficient file NAME.egm.cof is a binary file in little endian order. The first 8 bytes of this file must match the ID given in NAME.egm. This is followed by 2 sets of spherical harmonic coefficients. The first of these gives the gravity potential and the second gives the zeta-to-N corrections to the geoid height. The format for each set of coefficients is:

• N, the maximum degree of the sum stored as a 4-byte signed integer. This must satisfy N ≥ −1.
• M, the maximum order of the sum stored as a 4-byte signed integer. This must satisfy NM ≥ −1.
• Cnm, the coefficients of the cosine coefficients of the sum in column (i.e., m) major order. There are (M + 1) (2N - M + 2) / 2 elements which are stored as IEEE doubles (8 bytes). For example for N = M = 3, there are 10 coefficients arranged as C00, C10, C20, C30, C11, C21, C31, C22, C32, C33.
• Snm, the coefficients of the sine coefficients of the sum in column (i.e., m) major order starting at m = 1. There are M (2N - M + 1) / 2 elements which are stored as IEEE doubles (8 bytes). For example for N = M = 3, there are 6 coefficients arranged as S11, S21, S31, S22, S32, S33.

Although the coefficient file is in little endian order, GeographicLib can read it on big endian machines. It can only be read on machines which store doubles in IEEE format.

As an illustration, here is egm2008.egm:

EGMF-1
# An Earth Gravity Model (Format 1) file.  For documentation on the
# format of this file see
# http://geographiclib.sourceforge.net/html/gravity.html#gravityformat
Name            egm2008
Publisher       National Geospatial Intelligence Agency
Description     Earth Gravity Model 2008
URL             http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008
ReleaseDate     2008-06-01
ConversionDate  2011-11-19
DataVersion     1
ModelMass       3986004.415e8
AngularVelocity 7292115e-11
ReferenceMass   3986004.418e8
Flattening      1/298.257223563
HeightOffset    -0.41

# Gravitational and correction coefficients taken from
# EGM2008_to2190_TideFree and Zeta-to-N_to2160_egm2008 from
# the egm2008 distribution.
ID              EGM2008A


The code to produce the coefficient files for the wgs84 and grs80 models is

// Write the coefficient files needed for approximating the normal gravity
// field with a GravityModel. WARNING: this creates files, wgs84.egm.cof and
// grs80.egm.cof, in the current directory.
#include <cmath>
#include <fstream>
#include <iostream>
using namespace std;
using namespace GeographicLib;
int main() {
try {
const char* filenames[] = {"wgs84.egm.cof", "grs80.egm.cof"};
const char* ids[] = {"WGS1984A", "GRS1980A"};
for (int grs80 = 0; grs80 < 2; ++grs80) {
ofstream file(filenames[grs80], ios::binary);
Utility::writearray<char, char, false>(file, ids[grs80], 8);
const int N = 20, M = 0,
cnum = (M + 1) * (2 * N - M + 2) / 2; // cnum = N + 1
vector<int> num(2);
num[0] = N; num[1] = M;
Utility::writearray<int, int, false>(file, num);
vector<Math::real> c(cnum, 0);
const NormalGravity& earth(grs80 ? NormalGravity::GRS80() :
for (int n = 2; n <= N; n += 2)
c[n] = - earth.DynamicalFormFactor(n) / sqrt(Math::real(2*n + 1));
Utility::writearray<double, Math::real, false>(file, c);
num[0] = num[1] = -1;
Utility::writearray<int, int, false>(file, num);
}
}
catch (const exception& e) {
cerr << "Caught exception: " << e.what() << "\n";
return 1;
}
catch (...) {
cerr << "Caught unknown exception\n";
return 1;
}
}

# Comments on the NGA harmonic synthesis code

GravityModel attempts to reproduce the results of NGA's harmonic synthesis code for EGM2008, hsynth_WGS84.f. Listed here are issues that I encountered using the NGA code:

1. A compiler which allocates local variables on the stack produces an executable with just returns NaNs. The problem here is a missing SAVE statement in subroutine latf.
2. Because the model and references masses for egm2008 differ (by about 1 part in 109), there should be a 1/r contribution to the disturbing potential T. However, this term is set to zero in hsynth_WGS84 (effectively altering the normal potential). This shifts the surface W = U0 outward by about 5 mm. Note too that the reference ellipsoid is no longer a level surface of the effective normal potential.
3. Subroutine radgrav computes the ellipsoidal coordinate β incorrectly. This leads to small errors in the deflection of the vertical, ξ and η, when the height above the ellipsoid, h, is non-zero (about 10−7 arcsec at h = 400 km).
4. There are several expressions which will return inaccurate results due to cancellation. For example, subroutine grs computes the flattening using f = 1 − sqrt(1 − e2). Much better is to use f = e2/(1 + sqrt(1 − e2)). The expressions for q and q' in grs and radgrav suffer from similar problems. The resulting errors are tiny (about 50 pm in the geoid height); however, given that's there's essentially no cost to using more accurate expressions, it's preferable to do so.
5. hsynth_WGS84 returns an "undefined" value for xi and eta at the poles. Better would be to return the value obtained by taking the limit lat → ± 90°.

Issues 1–4 have been reported to the authors of hsynth_WGS84. Issue 1 is peculiar to Fortran and is not encountered in C++ code and GravityModel corrects issues 3–5. On issue 2, GravityModel neglects the 1/r term in T in GravityModel::GeoidHeight and GravityModel::SphericalAnomaly in order to produce results which match NGA's for these quantities. On the other hand, GravityModel::Disturbance and GravityModel::T do include this term.

# Details of the geoid height and anomaly calculations

Ideally, the geoid represents a surface of constant gravitational potential which approximates mean sea level. In reality some approximations are taken in determining this surface. The steps taking by GravityModel in computing the geoid height are described here (in the context of EGM2008). This mimics NGA's code hsynth_WGS84 closely because most users of EGM2008 use the gridded data generated by this code (e.g., Geoid) and it is desirable to use a consistent definition of the geoid height.

• The model potential is band limited; the minimum wavelength that is represented is 360°/2160 = 10' (i.e., about 10NM or 18.5km). The maximum degree for the spherical harmonic sum is 2190; however the model was derived using ellipsoidal harmonics of degree and order 2160 and the degree was increased to 2190 in order to capture the ellipsoidal harmonics faithfully with spherical harmonics.
• The 1/r term is omitted from the T (this is issue 2 in Comments on the NGA harmonic synthesis code). This moves the equipotential surface by about 5mm.
• The surface W = U0 is determined by Bruns' formula, which is roughly equivalent to a single iteration of Newton's method. The RMS error in this approximation is about 1.5mm with a maximum error of about 10 mm.
• The model potential is only valid above the earth's surface. A correction therefore needs to be included where the geoid lies beneath the terrain. This is NGA's "zeta-to-N" correction, which is represented by a spherical harmonic sum of degree and order 2160 (and so it misses short wavelength terrain variations). In addition, it entails estimating the isostatic equilibrium of the earth's crust. The correction lies in the range [−5.05 m, 0.05 m], however for 99.9% of the earth's surface the correction is less than 10 mm in magnitude.
• The resulting surface lies above the observed mean sea level, so −0.41m is added to the geoid height. (Better would be to change the potential used to define the geoid; but this would only change the result by about 2mm.)

A useful discussion of the problems with defining a geoid is given by Dru A. Smith in There is no such thing as "The" EGM96 geoid: Subtle points on the use of a global geopotential model, IGeS Bulletin No. 8, International Geoid Service, Milan, Italy, pp. 17–28 (1998).

GravityModel::GeoidHeight reproduces the results of the several NGA codes for harmonic synthesis with the following maximum discrepancies:

• egm84 = 1.1mm. This is probably due to inconsistent parameters for the reference ellipsoid in the NGA code. (In particular, the value of mass constant excludes the atmosphere; however, it's not clear whether the other parameters have been correspondingly adjusted.) Note that geoid heights predicted by egm84 differ from those of more recent gravity models by about 1 meter.
• egm96 = 23nm.
• egm2008 = 78pm. After addressing some of the issues alluded to in issue 4 in Comments on the NGA harmonic synthesis code, the maximum discrepancy becomes 23pm.

The formula for the gravity anomaly vector involves computing gravity and normal gravity at two different points (with the displacement between the points unknown ab initio). Since the gravity anomaly is already a small quantity it is sometimes acceptable to employ approximations that change the quantities by O(f). The NGA code uses the spherical approximation described by Heiskanen and Moritz, Sec. 2-14 and GravityModel::SphericalAnomaly uses the same approximation for compatibility. In this approximation, the gravity disturbance delta = grad T is calculated. Here, T once again excludes the 1/r term (this is issue 2 in Comments on the NGA harmonic synthesis code and is consistent with the computation of the geoid height). Note that delta compares the gravity and the normal gravity at the same point; the gravity anomaly vector is then found by estimating the gradient of the normal gravity in the limit that the earth is spherically symmetric. delta is expressed in spherical coordinates as deltax, deltay, deltaz where, for example, deltaz is the radial component of delta (not the component perpendicular to the ellipsoid) and deltay is similarly slightly different from the usual northerly component. The components of the anomaly are then given by

• gravity anomaly, Dg01 = deltaz − 2T/R, where R distance to the center of the earth;
• northerly component of the deflection of the vertical, xi = − deltay/gamma, where gamma is the magnitude of the normal gravity;
• easterly component of the deflection of the vertical, eta = − deltax/gamma.

NormalGravity computes the normal gravity accurately and avoids issue 3 of Comments on the NGA harmonic synthesis code. Thus while GravityModel::SphericalAnomaly reproduces the results for xi and eta at h = 0, there is a slight discrepancy if h is non-zero.

# The effect of the mass of the atmosphere

All of the supported models use WGS84 for the reference ellipsoid. This has (see TR8350.2, table 3.1)

• a = 6378137 m
• f = 1/298.257223563
• ω = 7292115 × 10−11 rad s−1
• GM = 3986004.418 × 108 m3 s−2.

The value of GM includes the mass of the atmosphere and so strictly only applies above the earth's atmosphere. Near the surface of the earth, the value of g will be less (in absolute value) than the value predicted by these models by about δg = (4πG/g) A = 8.552 × 10−11 A m2/kg, where G is the gravitational constant, g is the earth's gravity, and A is the pressure of the atmosphere. At sea level we have A = 101.3 kPa, and δg = 8.7−×−10−6 m s−2, approximately. (In other words the effect is about 1 part in a million; by way of comparison, buoyancy effects are about 100 times larger.)

# Geoid heights on a multi-processor system

The egm2008 model includes many terms (over 2 million spherical harmonics). For that reason computations using this model may be slow; for example it takes about 78 ms to compute the geoid height at a single point. There are two ways to speed up this computation:

• Use a GravityCircle to compute the geoid height at several points on a circle of latitude. This reduces the cost per point to about 92 μs (a reduction by a factor of over 800).
• Compute the values on several circles of latitude in parallel. One of the simplest ways of doing this is with OpenMP; on an 8-processor system, this can speed up the computation by another factor of 8.

Both of these techniques are illustrated by the following code, which computes a table of geoid heights on a regular grid and writes on the result in a .gtx file. On an 8-processor Intel 2.66 GHz machine using OpenMP (-DHAVE_OPENMP=1), it takes about 40 minutes of elapsed time to compute the geoid height for EGM2008 on a 1' grid. (Without these optimizations, the computation would have taken about 200 days!)

// Write out a gtx file of geoid heights above the ellipsoid. For egm2008 at
// 1' resolution this takes about 40 mins on a 8-processor Intel 2.66 GHz
// machine using OpenMP.
//
// For the format of gtx files, see
// http://vdatum.noaa.gov/dev/gtx_info.html#dev_gtx_binary
//
// data is binary big-endian:
// south latitude edge (degrees double)
// west longitude edge (degrees double)
// delta latitude (degrees double)
// delta longitude (degrees double)
// nlat = number of latitude rows (integer)
// nlong = number of longitude columns (integer)
// nlat * nlong geoid heights (meters float)
#include <vector>
#include <iostream>
#include <fstream>
#include <string>
#include <algorithm>
#if defined(_OPENMP)
#define HAVE_OPENMP 1
#else
#define HAVE_OPENMP 0
#endif
#if HAVE_OPENMP
# include <omp.h>
#endif
using namespace std;
using namespace GeographicLib;
int main(int argc, const char* const argv[]) {
// Hardwired for 3 args:
// 1 = the gravity model (e.g., egm2008)
// 2 = intervals per degree
// 3 = output GTX file
if (argc != 4) {
cerr << "Usage: " << argv[0]
<< " gravity-model intervals-per-degree output.gtx\n";
return 1;
}
try {
// Will need to set the precision for each thread, so save return value
int ndigits = Utility::set_digits();
string model(argv[1]);
// Number of intervals per degree
int ndeg = Utility::val<int>(string(argv[2]));
string filename(argv[3]);
GravityModel g(model);
int
nlat = 180 * ndeg + 1,
nlon = 360 * ndeg;
delta = 1 / Math::real(ndeg), // Grid spacing
latorg = -90,
lonorg = -180;
// Write results as floats in binary mode
ofstream file(filename.c_str(), ios::binary);
{
Math::real transform[] = {latorg, lonorg, delta, delta};
unsigned sizes[] = {unsigned(nlat), unsigned(nlon)};
Utility::writearray<double, Math::real, true>(file, transform, 4);
Utility::writearray<unsigned, unsigned, true>(file, sizes, 2);
}
// Compute and store results for nbatch latitudes at a time
const int nbatch = 64;
vector< vector<float> > N(nbatch, vector<float>(nlon));
for (int ilat0 = 0; ilat0 < nlat; ilat0 += nbatch) { // Loop over batches
int nlat0 = min(nlat, ilat0 + nbatch);
#if HAVE_OPENMP
# pragma omp parallel for
#endif
for (int ilat = ilat0; ilat < nlat0; ++ilat) { // Loop over latitudes
Utility::set_digits(ndigits); // Set the precision
lat = latorg + (ilat / ndeg) + delta * (ilat - ndeg * (ilat / ndeg)),
h = 0;
for (int ilon = 0; ilon < nlon; ++ilon) { // Loop over longitudes
Math::real lon = lonorg
+ (ilon / ndeg) + delta * (ilon - ndeg * (ilon / ndeg));
N[ilat - ilat0][ilon] = float(c.GeoidHeight(lon));
} // longitude loop
} // latitude loop -- end of parallel section
for (int ilat = ilat0; ilat < nlat0; ++ilat) // write out data
Utility::writearray<float, float, true>(file, N[ilat - ilat0]);
} // batch loop
}
catch (const exception& e) {
cerr << "Caught exception: " << e.what() << "\n";
return 1;
}
catch (...) {
cerr << "Caught unknown exception\n";
return 1;
}
}

cmake will add in support for OpenMP for examples/GeoidToGTX.cpp, if it is available.

# Normal gravity

The NormalGravity class computes "normal gravity" which refers to the exact (classical) solution of an idealised system consisting of an ellipsoid of revolution with the following properties:

• M the mass interior to the ellipsoid,
• b the polar semi-axis,
• ω the rotation rate about the polar axis.

(N.B. The mass always appears in the combination GM, with units m3/s2, where G is the gravtitational constant.) The distribution of the mass M within the ellipsoid is such that the surface of the ellipsoid is at a constant normal potential where the normal potential is the sum of the gravitational potential (due to the gravitional attraction) and the centrifugal potention (due to the rotation). The resulting field exterior to the ellipsoid is called normal gravity and was found by Somigliana (1929). Because the potential is constant on the ellipsoid it is sometimes referred to as the level ellipsoid.

References:

• C. Somigliana, Teoria generale del campo gravitazionale dell'ellissoide di rotazione, Mem. Soc. Astron. Ital, 4, 541–599 (1929).
• W. A. Heiskanen and H. Moritz, Physical Geodesy (Freeman, San Francisco, 1967), Secs. 1-19, 2-7, 2-8 (2-9, 2-10), 6-2 (6-3).
• B. Hofmann-Wellenhof, H. Moritz, Physical Geodesy (Second edition, Springer, 2006). https://doi.org/10.1007/978-3-211-33545-1
• H. Moritz, Geodetic Reference System 1980, J. Geodesy 54(3), 395–405 (1980). https://doi.org/10.1007/BF02521480
• H. Moritz, The Figure of the Earth (Wichmann, 1990).
• M. Chasles, Solution nouvelle du problème de l’attraction d’un ellipsoïde hétérogène sur un point exterieur, Jour. Liouville 5, 465–488 (1840), Note 2. http://sites.mathdoc.fr/JMPA/PDF/JMPA_1840_1_5_A41_0.pdf

## Ellipsoidal coordinates

Two set of formulas are presented: those of Heiskanen and Moritz (1967) which are applicable to an oblate ellipsoid and a second set where the variables are distinguished with primes which apply to a prolate ellipsoid. The primes are omitted from those variables which are the same in the two cases. In the text, the parenthetical "resp." clauses apply to prolate ellipsoids.

Cylindrical coordinates $$R,Z$$ are expressed in terms of ellipsoidal coordinates

\begin{align} R &= \sqrt{u^2 + E^2} \cos\beta = u' \cos\beta,\\ Z &= u \sin\beta = \sqrt{u'^2 + E'^2} \sin\beta, \end{align}

where

\begin{align} E^2 = a^2 - b^2,&{} \qquad u^2 = u'^2 + E'^2,\\ E'^2 = b^2 - a^2,&{} \qquad u'^2 = u^2 + E^2. \end{align}

Surfaces of constant $$u$$ (or $$u'$$) are confocal ellipsoids. The surface $$u = 0$$ (resp. $$u' = 0$$) corresponds to the focal disc of diameter $$2E$$ (resp. focal rod of length $$2E'$$). The level ellipsoid is given by $$u = b$$ (resp. $$u' = a$$). On the level ellipsoid, $$\beta$$ is the familiar parametric latitude.

In writing the potential and the gravity, it is useful to introduce the functions

\begin{align} Q(z) &= \frac1{2z^3} \biggl[\biggl(1 + \frac3{z^2}\biggr)\tan^{-1}z - \frac3z\biggr],\\ Q'(z') &= \frac{(1+z'^2)^{3/2}}{2z'^3} \biggl[\biggl(2 + \frac3{z'^2}\biggr)\sinh^{-1}z' - \frac{3\sqrt{1+z'^2}}{z'}\biggr],\\ G(z) &= \biggl(3Q(z)+z\frac{dQ(z)}{dz}\biggr)(1+z^2)\\ &= \frac1{z^4}\biggl[3(1+z^2) \biggl(1-\frac{\tan^{-1}z}z\biggr)-z^2\biggr],\\ G'(z') &= \frac{3Q'(z')}{1+z'^2}+z'\frac{dQ'(z')}{dz'}\\ &= \frac{1+z'^2}{z'^4} \biggl[3\biggl(1-\frac{\sqrt{1+z'^2}\sinh^{-1}z'}{z'}\biggr) +z'^2\biggr]. \end{align}

The function arguments are $$z = E/u$$ (resp. $$z' = E'/u'$$) and the unprimed and primed quantities are related by

\begin{align} Q'(z') = Q(z),&{} \qquad G'(z') = G(z),\\ z'^2 = -\frac{z^2}{1 + z^2},&{} \qquad z^2 = -\frac{z'^2}{1 + z'^2}. \end{align}

The functions $$q(u)$$ and $$q'(u)$$ used by Heiskanen and Moritz are given by

$q(u) = \frac{E^3}{u^3}Q\biggl(\frac Eu\biggr),\qquad q'(u) = \frac{E^2}{u^2}G\biggl(\frac Eu\biggr).$

The functions $$Q(z)$$, $$Q'(z')$$, $$G(z)$$, and $$G'(z')$$ are more convenient for use in numerical codes because they are finite in the spherical limit $$E \rightarrow 0$$, i.e., $$Q(0) = Q'(0) = \frac2{15}$$, and $$G(0) = G'(0) = \frac25$$.

## The normal potential

The normal potential is the sum of three components, a mass term, a quadrupole term and a centrifugal term,

$U = U_m + U_q + U_r.$

The mass term is

$U_m = \frac {GM}E \tan^{-1}\frac Eu = \frac {GM}{E'} \sinh^{-1}\frac{E'}{u'};$

\begin{align} U_q &= \frac{\omega^2}2 \frac{a^2 b^3}{u^3} \frac{Q(E/u)}{Q(E/b)} \biggl(\sin^2\beta-\frac13\biggr)\\ &= \frac{\omega^2}2 \frac{a^2 b^3}{(u'^2+E'^2)^{3/2}} \frac{Q'(E'/u')}{Q'(E'/a)} \biggl(\sin^2\beta-\frac13\biggr); \end{align}

finally, the rotational term is

$U_r = \frac{\omega^2}2 R^2 = \frac{\omega^2}2 (u^2 + E^2) \cos^2\beta = \frac{\omega^2}2 u'^2 \cos^2\beta.$

$$U_m$$ and $$U_q$$ are both gravitational potentials (due to mass within the ellipsoid). The total mass contributing to $$U_m$$ is $$M$$; the total mass contributing to $$U_q$$ vanishes (far from the ellipsoid, the $$U_q$$ decays inversely as the cube of the distance).

$$U_m$$ and $$U_q + U_r$$ are separately both constant on the level ellipsoid. $$U_m$$ is the normal potential for a massive non-rotating ellipsoid. $$U_q + U_r$$ is the potential for a massless rotating ellipsoid. Combining all the terms, $$U$$ is the normal potential for a massive rotating ellipsoid.

## The mass distribution

Typically, the normal potential, $$U$$, is only of interest for outside the ellipsoid $$u \ge b$$ (resp. $$u' \ge a$$). In planetary applications a open problem is finding a mass distribution which is in hydrostatic equilibrium (the mass density is non-negative and a non-decreasing function of the potential interior to the ellipsoid).

However it is possible to give singular mass distributions consistent with the normal potential.

For a non-rotating body, the potential $$U = U_m$$ is generated by a sandwiching the mass $$M$$ uniformly between the level ellipsoid with semi-axes $$a$$ and $$b$$ and a close similar ellipsoid with semi-axes $$(1-\epsilon)a$$ and $$(1-\epsilon)b$$. Chasles (1840) extends a theorem of Newton to show that the field interior to such an ellipsoidal shell vanishes. Thus the potential on the ellipsoid is constant, i.e., it is indeed a level ellipsoid. This result also holds for a non-rotating triaxial ellipsoid.

Observing that $$U_m$$ and $$U_q$$ as defined above obey $$\nabla^2 U_m = \nabla^2 U_q = 0$$ everywhere for $$u > 0$$ (resp. $$u' > 0$$), we see that these potentials correspond to masses concentrated at $$u = 0$$ (resp. $$u' = 0$$).

In the oblate case, $$U_m$$ is generated by a massive disc at $$Z = 0$$, $$R < E$$, with mass density (mass per unit area) $$\rho_m$$ and moments of inertia about the equatorial (resp. polar) axis of $$A_m$$ (resp. $$C_m$$) given by

\begin{align} \rho_m &= \frac M{2\pi E\sqrt{E^2 - R^2}},\\ A_m &= \frac {ME^2}3, \\ C_m &= \frac {2ME^2}3, \\ C_m-A_m &= \frac {ME^2}3. \end{align}

This mass distribution is the same as that produced by projecting a uniform spheric shell of mass $$M$$ and radius $$E$$ onto the equatorial plane.

In the prolate case, $$U_m$$ is generated by a massive rod at $$R = 0$$, $$Z < E'$$ and now the mass density $$\rho'_m$$ has units mass per unit length,

\begin{align} \rho'_m &= \frac M{2E'},\\ A_m &= \frac {ME'^2}3, \\ C_m &= 0, \\ C_m-A_m &= -\frac {ME'^2}3. \end{align}

This mass distribution is the same as that produced by projecting a uniform spheric shell of mass $$M$$ and radius $$E'$$ onto the polar axis.

Similarly, $$U_q$$ is generated in the oblate case by

\begin{align} \rho_q &= \frac{a^2 b^3 \omega^2}G \frac{2E^2 - 3R^2}{6\pi E^5 \sqrt{E^2 - R^2} Q(E/b)}, \\ A_q &= -\frac{a^2 b^3 \omega^2}G \frac2{45 Q(E/b)}, \\ C_q &= -\frac{a^2 b^3 \omega^2}G \frac4{45 Q(E/b)}, \\ C_q-A_q &= -\frac{a^2 b^3 \omega^2}G \frac2{45 Q(E/b)}. \end{align}

The corresponding results for a prolate ellipsoid are

\begin{align} \rho_q' &= \frac{a^2 b^3 \omega^2}G \frac{3Z^2 - E'^2}{12 E'^5 Q'(E'/a)}, \\ A_q &= \frac{a^2 b^3 \omega^2}G \frac2{45 Q'(E'/a)}, \\ C_q &= 0, \\ C_q-A_q &= -\frac{a^2 b^3 \omega^2}G \frac2{45 Q'(E'/a)}. \end{align}

Summing up the mass and quadrupole terms, we have

\begin{align} A &= A_m + A_q, \\ C &= C_m + C_q, \\ J_2 & = \frac{C - A}{Ma^2}, \end{align}

where $$J_2$$ is the dynamical form factor.

## The surface gravity

Each term in the potential contributes to the gravity on the surface of the ellipsoid

$\gamma = \gamma_m + \gamma_q + \gamma_r;$

These are the components of gravity normal to the ellipsoid and, by convention, $$\gamma$$ is positive downwards. The tangential components of the total gravity and that due to $$U_m$$ vanish. Those tangential components of the gravity due to $$U_q$$ and $$U_r$$ cancel one another.

The gravity $$\gamma$$ has the following dependence on latitude

\begin{align} \gamma &= \frac{b\gamma_a\cos^2\beta + a\gamma_b\sin^2\beta} {\sqrt{b^2\cos^2\beta + a^2\sin^2\beta}}\\ &= \frac{a\gamma_a\cos^2\phi + b\gamma_b\sin^2\phi} {\sqrt{a^2\cos^2\phi + b^2\sin^2\phi}}, \end{align}

and the individual components, $$\gamma_m$$, $$\gamma_q$$, and $$\gamma_r$$, have the same dependence on latitude. The equatorial and polar gravities are

\begin{align} \gamma_a &= \gamma_{ma} + \gamma_{qa} + \gamma_{ra},\\ \gamma_b &= \gamma_{mb} + \gamma_{qb} + \gamma_{rb}, \end{align}

where

\begin{align} \gamma_{ma} &= \frac{GM}{ab},\qquad \gamma_{mb} = \frac{GM}{a^2},\\ \gamma_{qa} &= -\frac{\omega^2 a}6 \frac{G(E/b)}{Q(E/b)} = -\frac{\omega^2 a}6 \frac{G'(E'/a)}{Q'(E'/a)},\\ \gamma_{qb} &= -\frac{\omega^2 b}3 \frac{G(E/b)}{Q(E/b)} = -\frac{\omega^2 b}3 \frac{G'(E'/a)}{Q'(E'/a)},\\ \gamma_{ra} &= -\omega^2 a,\qquad \gamma_{rb} = 0. \end{align}

## The mean gravity

Performing an average of the surface gravity over the area of the ellipsoid gives

$\langle \gamma \rangle = \frac {4\pi a^2 b}A \biggl(\frac{2\gamma_a}{3a} + \frac{\gamma_b}{3b}\biggr),$

where $$A$$ is the area of the ellipsoid

\begin{align} A &= 2\pi\biggl( a^2 + ab\frac{\sinh^{-1}(E/b)}{E/b} \biggr)\\ &= 2\pi\biggl( a^2 + b^2\frac{\tan^{-1}(E'/a)}{E'/a} \biggr). \end{align}

The contributions to the mean gravity are

\begin{align} \langle \gamma_m \rangle &= \frac{4\pi}A GM, \\ \langle \gamma_q \rangle &= 0 \quad \text{(as expected)}, \\ \langle \gamma_r \rangle &= -\frac{4\pi}A \frac{2\omega^2 a^2b}3,\\ \end{align}

resulting in

$\langle \gamma \rangle = \frac{4\pi}A \biggl(GM - \frac{2\omega^2 a^2b}3\biggr).$

## Possible values of the dynamical form factor

The solution for the normal gravity is well defined for arbitrary $$M$$, $$\omega$$, $$a > 0$$, and $$f < 1$$. (Note that arbitrary oblate and prolate ellipsoids are possible, although hydrostatic equilibrium would not result in a prolate ellipsoid.) However, if is much easier to measure the dynamical form factor $$J_2$$ (from the motion of artificial satellites) than the flattening $$f$$. (Note too that $$GM$$ is also typically measured from from satellite or astronomical observations and so it includes the mass of the atmosphere.)

So a question for the software developer is: given values of $$M > 0$$, $$\omega$$, and $$a > 0$$, what are the allowed values of $$J_2$$? We restrict the question to $$M > 0$$. The (unphysical) case $$M = 0$$ is problematic because $$M$$ appears in the denominator in the definition of $$J_2$$. In the (also unphysical) case $$M < 0$$, a given $$J_2$$ can result from two distinct values of $$f$$.

Holding $$M > 0$$, $$\omega$$, and $$a > 0$$ fixed and varying $$f$$ from $$-\infty$$ to $$1$$, we find that $$J_2$$ monotonically increases from $$-\infty$$ to

$\frac13 - \frac8{45\pi} \frac{\omega^2 a^3}{GM}.$

Thus any value of $$J_2$$ less that this value is permissible (but some of these values may be unphysical). In obtaining this limiting value, we used the result $$Q(z \rightarrow \infty) \rightarrow \pi/(4 z^3)$$. The value

$J_2 = -\frac13 \frac{\omega^2 a^3}{GM}$

results in a sphere ( $$f = 0$$).

Back to Geoid height. Forward to Magnetic models. Up to Contents.