/usr/share/doc/geographiclib/html/gravity.html is in geographiclib-doc 1.49-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 | <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="Content-Type" content="text/xhtml;charset=UTF-8"/>
<meta http-equiv="X-UA-Compatible" content="IE=9"/>
<meta name="generator" content="Doxygen 1.8.13"/>
<meta name="viewport" content="width=device-width, initial-scale=1"/>
<title>GeographicLib: Gravity models</title>
<link href="tabs.css" rel="stylesheet" type="text/css"/>
<script type="text/javascript" src="jquery.js"></script>
<script type="text/javascript" src="dynsections.js"></script>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
extensions: ["tex2jax.js"],
jax: ["input/TeX","output/HTML-CSS"],
});
</script><script type="text/javascript" src="/usr/share/javascript/mathjax/MathJax.js/MathJax.js"></script>
<link href="doxygen.css" rel="stylesheet" type="text/css" />
</head>
<body>
<div id="top"><!-- do not remove this div, it is closed by doxygen! -->
<div id="titlearea">
<table cellspacing="0" cellpadding="0">
<tbody>
<tr style="height: 56px;">
<td id="projectalign" style="padding-left: 0.5em;">
<div id="projectname">GeographicLib
 <span id="projectnumber">1.49</span>
</div>
</td>
</tr>
</tbody>
</table>
</div>
<!-- end header part -->
<!-- Generated by Doxygen 1.8.13 -->
<script type="text/javascript" src="menudata.js"></script>
<script type="text/javascript" src="menu.js"></script>
<script type="text/javascript">
$(function() {
initMenu('',false,false,'search.php','Search');
});
</script>
<div id="main-nav"></div>
</div><!-- top -->
<div class="header">
<div class="headertitle">
<div class="title">Gravity models </div> </div>
</div><!--header-->
<div class="contents">
<div class="textblock"><center> Back to <a class="el" href="geoid.html">Geoid height</a>. Forward to <a class="el" href="normalgravity.html">Normal gravity</a>. Up to <a class="el" href="index.html#contents">Contents</a>. </center><p><a class="el" href="namespaceGeographicLib.html" title="Namespace for GeographicLib. ">GeographicLib</a> can compute the earth's gravitational field with an earth gravity model using the <a class="el" href="classGeographicLib_1_1GravityModel.html" title="Model of the earth's gravity field. ">GravityModel</a> and <a class="el" href="classGeographicLib_1_1GravityCircle.html" title="Gravity on a circle of latitude. ">GravityCircle</a> classes and with the <a href="Gravity.1.html">Gravity</a> 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 <em>normal gravity</em> which is the exact solution for an idealized rotating ellipsoidal body. This is implemented with the <a class="el" href="classGeographicLib_1_1NormalGravity.html" title="The normal gravity of the earth. ">NormalGravity</a> class and some notes on are provided in section <a class="el" href="normalgravity.html">Normal gravity</a></p>
<p>Go to</p><ul>
<li><a class="el" href="gravity.html#gravityinst">Installing the gravity models</a></li>
<li><a class="el" href="gravity.html#gravityformat">The format of the gravity model files</a></li>
<li><a class="el" href="gravity.html#gravitynga">Comments on the NGA harmonic synthesis code</a></li>
<li><a class="el" href="gravity.html#gravitygeoid">Details of the geoid height and anomaly calculations</a></li>
<li><a class="el" href="gravity.html#gravityatmos">The effect of the mass of the atmosphere</a></li>
<li><a class="el" href="gravity.html#gravityparallel">Geoid heights on a multi-processor system</a></li>
<li><a class="el" href="normalgravity.html">Normal gravity</a> (now on its own page)</li>
</ul>
<p>The supported models are</p><ul>
<li><b>egm84</b>, the <a href="http://earth-info.nga.mil/GandG/wgs84/gravitymod/wgs84_180/wgs84_180.html">Earth Gravity Model 1984</a>, which includes terms up to degree 180.</li>
<li><b>egm96</b>, the <a href="http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm96/egm96.html">Earth Gravity Model 1996</a>, which includes terms up to degree 360.</li>
<li><b>egm2008</b>, the <a href="http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008">Earth Gravity Model 2008</a>, which includes terms up to degree 2190.</li>
<li><b>wgs84</b>, the <a href="http://earth-info.nga.mil/GandG/publications/tr8350.2/tr8350_2.html">WGS84 Reference Ellipsoid</a>. This is just reproduces the normal gravitational field for the reference ellipsoid. This includes the zonal coefficients up to order 20. Usually <a class="el" href="classGeographicLib_1_1NormalGravity.html#aa07cdd6fe8ce8a92f00b74ea93908a78">NormalGravity::WGS84()</a> should be used instead.</li>
<li><b>grs80</b>, the GRS80 Reference <a class="el" href="classGeographicLib_1_1Ellipsoid.html" title="Properties of an ellipsoid. ">Ellipsoid</a>. This is just reproduces the normal gravitational field for the GRS80 ellipsoid. This includes the zonal coefficients up to order 20. Usually <a class="el" href="classGeographicLib_1_1NormalGravity.html#adbe37628ec483ff0437ee13eab089b68">NormalGravity::GRS80()</a> should be used instead.</li>
</ul>
<p>See</p><ul>
<li>W. A. Heiskanen and H. Moritz, Physical Geodesy (Freeman, San Francisco, 1967).</li>
</ul>
<p>for more information.</p>
<p><b>Acknowledgment:</b> 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.</p>
<h1><a class="anchor" id="gravityinst"></a>
Installing the gravity models</h1>
<p>These gravity models are available for download: </p><center> <a class="anchor" id=""></a>
<table class="doxtable">
<caption>Available gravity models</caption>
<tr>
<th rowspan="2">name </th><th rowspan="2">max<br />
degree </th><th rowspan="2">size<br />
(kB) </th><th colspan="3"><center>Download Links (size, kB)</center> </th></tr>
<tr>
<th>tar file </th><th>Windows<br />
installer </th><th>zip file </th></tr>
<tr>
<td>egm84 </td><td><center>180</center> </td><td><center>27</center> </td><td><center> <a href="https://sourceforge.net/projects/geographiclib/files/gravity-distrib/egm84.tar.bz2">link</a> (26)</center> </td><td><center> <a href="https://sourceforge.net/projects/geographiclib/files/gravity-distrib/egm84.exe">link</a> (55)</center> </td><td><center> <a href="https://sourceforge.net/projects/geographiclib/files/gravity-distrib/egm84.zip">link</a> (26)</center> </td></tr>
<tr>
<td>egm96 </td><td><center>360</center> </td><td><center>2100</center> </td><td><center> <a href="https://sourceforge.net/projects/geographiclib/files/gravity-distrib/egm96.tar.bz2">link</a> (2100)</center> </td><td><center> <a href="https://sourceforge.net/projects/geographiclib/files/gravity-distrib/egm96.exe">link</a> (2300)</center> </td><td><center> <a href="https://sourceforge.net/projects/geographiclib/files/gravity-distrib/egm96.zip">link</a> (2100)</center> </td></tr>
<tr>
<td>egm2008 </td><td><center>2190</center> </td><td><center>76000</center> </td><td><center> <a href="https://sourceforge.net/projects/geographiclib/files/gravity-distrib/egm2008.tar.bz2">link</a> (75000)</center> </td><td><center> <a href="https://sourceforge.net/projects/geographiclib/files/gravity-distrib/egm2008.exe">link</a> (72000)</center> </td><td><center> <a href="https://sourceforge.net/projects/geographiclib/files/gravity-distrib/egm2008.zip">link</a> (73000)</center> </td></tr>
<tr>
<td>wgs84 </td><td><center>20</center> </td><td><center>1</center> </td><td><center> <a href="https://sourceforge.net/projects/geographiclib/files/gravity-distrib/wgs84.tar.bz2">link</a> (1)</center> </td><td><center> <a href="https://sourceforge.net/projects/geographiclib/files/gravity-distrib/wgs84.exe">link</a> (30)</center> </td><td><center> <a href="https://sourceforge.net/projects/geographiclib/files/gravity-distrib/wgs84.zip">link</a> (1)</center> </td></tr>
</table>
</center><p> The "size" column is the size of the uncompressed data.</p>
<p>For Linux and Unix systems, <a class="el" href="namespaceGeographicLib.html" title="Namespace for GeographicLib. ">GeographicLib</a> 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 </p><pre class="fragment"> geographiclib-get-gravity all # to install egm84, egm96, egm2008, wgs84
geographiclib-get-gravity -h # for help
</pre><p> 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.</p>
<p>Windows users should download and run the Windows installers. These will prompt for an installation directory with the default being </p><pre class="fragment"> C:/ProgramData/GeographicLib
</pre><p> (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.)</p>
<p>Otherwise download <em>either</em> the tar.bz2 file <em>or</em> the zip file (they have the same contents). To unpack these, run, for example </p><pre class="fragment"> 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.
</pre><p> and, again, the data will be placed in the "gravity" subdirectory.</p>
<p>However you install the gravity models, all the datasets should be installed in the same directory. <a class="el" href="classGeographicLib_1_1GravityModel.html" title="Model of the earth's gravity field. ">GravityModel</a> and <a href="Gravity.1.html">Gravity</a> uses a compile time default to locate the datasets. This is</p><ul>
<li>/usr/local/share/GeographicLib/gravity, for non-Windows systems</li>
<li>C:/ProgramData/GeographicLib/gravity, for Windows systems</li>
</ul>
<p>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 <a class="el" href="classGeographicLib_1_1GravityModel.html#a0fdf62e41828ae7ae183d9e876f37954">GravityModel::DefaultGravityPath()</a> for details. Finally, the path may be set using the optional second argument to the <a class="el" href="classGeographicLib_1_1GravityModel.html" title="Model of the earth's gravity field. ">GravityModel</a> constructor or with the "-d" flag to <a href="Gravity.1.html">Gravity</a>. Supplying the "-h" flag to <a href="Gravity.1.html">Gravity</a> 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.</p>
<h1><a class="anchor" id="gravityformat"></a>
The format of the gravity model files</h1>
<p>The constructor for <a class="el" href="classGeographicLib_1_1GravityModel.html" title="Model of the earth's gravity field. ">GravityModel</a> 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.</p>
<p>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").</p>
<p>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.</p>
<p><a class="el" href="classGeographicLib_1_1GravityModel.html" title="Model of the earth's gravity field. ">GravityModel</a> only pays attention to the following keywords</p><ul>
<li>keywords that affect the field calculation, namely:<ul>
<li><b>ModelRadius</b> (required), the normalizing radius for the model in meters.</li>
<li><b>ReferenceRadius</b> (required), the equatorial radius <em>a</em> for the reference ellipsoid meters.</li>
<li><b>ModelMass</b> (required), the mass constant <em>GM</em> for the model in meters<sup>3</sup>/seconds<sup>2</sup>.</li>
<li><b>ReferenceMass</b> (required), the mass constant <em>GM</em> for the reference ellipsoid in meters<sup>3</sup>/seconds<sup>2</sup>.</li>
<li><b>AngularVelocity</b> (required), the angular velocity ω for the model <em>and</em> the reference ellipsoid in rad seconds<sup>−1</sup>.</li>
<li><b>Flattening</b>, the flattening of the reference ellipsoid; this can be given a fraction, e.g., 1/298.257223563. One of <b>Flattening</b> and <b>DynamicalFormFactor</b> is required.</li>
<li><b>DynamicalFormFactor</b>, the dynamical form factor <em>J</em><sub>2</sub> for the reference ellipsoid. One of <b>Flattening</b> and <b>DynamicalFormFactor</b> is required.</li>
<li><b>HeightOffset</b> (default 0), the constant height offset (meters) added to obtain the geoid height.</li>
<li><b>CorrectionMultiplier</b> (default 1), the multiplier for the "zeta-to-N" correction terms for the geoid height to convert them to meters.</li>
<li><b>Normalization</b> (default full), the normalization used for the associated Legendre functions (full or schmidt).</li>
<li><b>ID</b> (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).</li>
</ul>
The parameters <b>ModelRadius</b>, <b>ModelMass</b>, and <b>AngularVelocity</b> apply to the gravity model, while <b>ReferenceRadius</b>, <b>ReferenceMass</b>, <b>AngularVelocity</b>, and either <b>Flattening</b> or <b>DynamicalFormFactor</b> characterize the reference ellipsoid. <b>AngularVelocity</b> (because it can be measured precisely) is the same for the gravity model and the reference ellipsoid. <b>ModelRadius</b> 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). <b>ModelMass</b> and <b>ReferenceMass</b> need not be the same and, indeed, they are slightly difference for egm2008. As a result the disturbing potential <em>T</em> has a 1/<em>r</em> dependence at large distances.</li>
<li>keywords that store data that the user can query:<ul>
<li><b>Name</b>, the name of the model.</li>
<li><b>Description</b>, a more descriptive name of the model.</li>
<li><b>ReleaseDate</b>, when the model was created.</li>
</ul>
</li>
<li>keywords that are examined to verify that their values are valid:<ul>
<li><b>ByteOrder</b> (default little), the order of bytes in the .egm.cof file. Only little endian is supported at present.</li>
</ul>
</li>
</ul>
<p>Other keywords are ignored.</p>
<p>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:</p><ul>
<li><em>N</em>, the maximum degree of the sum stored as a 4-byte signed integer. This must satisfy <em>N</em> ≥ −1.</li>
<li><em>M</em>, the maximum order of the sum stored as a 4-byte signed integer. This must satisfy <em>N</em> ≥ <em>M</em> ≥ −1.</li>
<li><em>C</em><sub><em>nm</em></sub>, the coefficients of the cosine coefficients of the sum in column (i.e., <em>m</em>) major order. There are (<em>M</em> + 1) (2<em>N</em> - <em>M</em> + 2) / 2 elements which are stored as IEEE doubles (8 bytes). For example for <em>N</em> = <em>M</em> = 3, there are 10 coefficients arranged as <em>C</em><sub>00</sub>, <em>C</em><sub>10</sub>, <em>C</em><sub>20</sub>, <em>C</em><sub>30</sub>, <em>C</em><sub>11</sub>, <em>C</em><sub>21</sub>, <em>C</em><sub>31</sub>, <em>C</em><sub>22</sub>, <em>C</em><sub>32</sub>, <em>C</em><sub>33</sub>.</li>
<li><em>S</em><sub><em>nm</em></sub>, the coefficients of the sine coefficients of the sum in column (i.e., <em>m</em>) major order starting at <em>m</em> = 1. There are <em>M</em> (2<em>N</em> - <em>M</em> + 1) / 2 elements which are stored as IEEE doubles (8 bytes). For example for <em>N</em> = <em>M</em> = 3, there are 6 coefficients arranged as <em>S</em><sub>11</sub>, <em>S</em><sub>21</sub>, <em>S</em><sub>31</sub>, <em>S</em><sub>22</sub>, <em>S</em><sub>32</sub>, <em>S</em><sub>33</sub>.</li>
</ul>
<p>Although the coefficient file is in little endian order, <a class="el" href="namespaceGeographicLib.html" title="Namespace for GeographicLib. ">GeographicLib</a> can read it on big endian machines. It can only be read on machines which store doubles in IEEE format.</p>
<p>As an illustration, here is egm2008.egm: </p><pre class="fragment">EGMF-1
# An Earth Gravity Model (Format 1) file. For documentation on the
# format of this file see
# https://geographiclib.sourceforge.io/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
ModelRadius 6378136.3
ModelMass 3986004.415e8
AngularVelocity 7292115e-11
ReferenceRadius 6378137
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
</pre><p>The code to produce the coefficient files for the wgs84 and grs80 models is </p><div class="fragment"><div class="line"><span class="comment">// Write the coefficient files needed for approximating the normal gravity</span></div><div class="line"><span class="comment">// field with a GravityModel. WARNING: this creates files, wgs84.egm.cof and</span></div><div class="line"><span class="comment">// grs80.egm.cof, in the current directory.</span></div><div class="line"></div><div class="line"><span class="preprocessor">#include <cmath></span></div><div class="line"><span class="preprocessor">#include <fstream></span></div><div class="line"><span class="preprocessor">#include <iostream></span></div><div class="line"><span class="preprocessor">#include <<a class="code" href="NormalGravity_8hpp.html">GeographicLib/NormalGravity.hpp</a>></span></div><div class="line"><span class="preprocessor">#include <<a class="code" href="Utility_8hpp.html">GeographicLib/Utility.hpp</a>></span></div><div class="line"></div><div class="line"><span class="keyword">using namespace </span><a class="code" href="namespacestd.html">std</a>;</div><div class="line"><span class="keyword">using namespace </span><a class="code" href="namespaceGeographicLib.html">GeographicLib</a>;</div><div class="line"></div><div class="line"><span class="keywordtype">int</span> <a class="code" href="CartConvert_8cpp.html#a814244b02f6701fd8ead246678bb3dae">main</a>() {</div><div class="line"> <span class="keywordflow">try</span> {</div><div class="line"> <a class="code" href="classGeographicLib_1_1Utility.html#ab95e26081d2c69019703adb8a4e4d3d0">Utility::set_digits</a>();</div><div class="line"> <span class="keyword">const</span> <span class="keywordtype">char</span>* filenames[] = {<span class="stringliteral">"wgs84.egm.cof"</span>, <span class="stringliteral">"grs80.egm.cof"</span>};</div><div class="line"> <span class="keyword">const</span> <span class="keywordtype">char</span>* ids[] = {<span class="stringliteral">"WGS1984A"</span>, <span class="stringliteral">"GRS1980A"</span>};</div><div class="line"> <span class="keywordflow">for</span> (<span class="keywordtype">int</span> grs80 = 0; grs80 < 2; ++grs80) {</div><div class="line"> ofstream file(filenames[grs80], ios::binary);</div><div class="line"> Utility::writearray<char, char, false>(file, ids[grs80], 8);</div><div class="line"> <span class="keyword">const</span> <span class="keywordtype">int</span> N = 20, M = 0,</div><div class="line"> cnum = (M + 1) * (2 * N - M + 2) / 2; <span class="comment">// cnum = N + 1</span></div><div class="line"> vector<int> num(2);</div><div class="line"> num[0] = N; num[1] = M;</div><div class="line"> Utility::writearray<int, int, false>(file, num);</div><div class="line"> vector<Math::real> c(cnum, 0);</div><div class="line"> <span class="keyword">const</span> <a class="code" href="classGeographicLib_1_1NormalGravity.html">NormalGravity</a>& earth(grs80 ? <a class="code" href="classGeographicLib_1_1NormalGravity.html#adbe37628ec483ff0437ee13eab089b68">NormalGravity::GRS80</a>() :</div><div class="line"> <a class="code" href="classGeographicLib_1_1NormalGravity.html#aa07cdd6fe8ce8a92f00b74ea93908a78">NormalGravity::WGS84</a>());</div><div class="line"> <span class="keywordflow">for</span> (<span class="keywordtype">int</span> n = 2; n <= N; n += 2)</div><div class="line"> c[n] = - earth.DynamicalFormFactor(n) / sqrt(<a class="code" href="classGeographicLib_1_1Math.html#aeee4778d7cf2f9fb9648efe4911da59d">Math::real</a>(2*n + 1));</div><div class="line"> Utility::writearray<double, Math::real, false>(file, c);</div><div class="line"> num[0] = num[1] = -1;</div><div class="line"> Utility::writearray<int, int, false>(file, num);</div><div class="line"> }</div><div class="line"> }</div><div class="line"> <span class="keywordflow">catch</span> (<span class="keyword">const</span> exception& e) {</div><div class="line"> cerr << <span class="stringliteral">"Caught exception: "</span> << e.what() << <span class="stringliteral">"\n"</span>;</div><div class="line"> <span class="keywordflow">return</span> 1;</div><div class="line"> }</div><div class="line"> <span class="keywordflow">catch</span> (...) {</div><div class="line"> cerr << <span class="stringliteral">"Caught unknown exception\n"</span>;</div><div class="line"> <span class="keywordflow">return</span> 1;</div><div class="line"> }</div><div class="line">}</div></div><!-- fragment --><h1><a class="anchor" id="gravitynga"></a>
Comments on the NGA harmonic synthesis code</h1>
<p><a class="el" href="classGeographicLib_1_1GravityModel.html" title="Model of the earth's gravity field. ">GravityModel</a> 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:</p><ol type="1">
<li>A compiler which allocates local variables on the stack produces an executable with just returns NaNs. The problem here is a missing <code>SAVE</code> statement in subroutine <code>latf</code>.</li>
<li>Because the model and references masses for egm2008 differ (by about 1 part in 10<sup>9</sup>), there should be a 1/<em>r</em> contribution to the disturbing potential <em>T</em>. However, this term is set to zero in hsynth_WGS84 (effectively altering the normal potential). This shifts the surface <em>W</em> = <em>U</em><sub>0</sub> outward by about 5 mm. Note too that the reference ellipsoid is no longer a level surface of the effective normal potential.</li>
<li>Subroutine <code>radgrav</code> computes the ellipsoidal coordinate β incorrectly. This leads to small errors in the deflection of the vertical, ξ and η, when the height above the ellipsoid, <em>h</em>, is non-zero (about 10<sup>−7</sup> arcsec at <em>h</em> = 400 km).</li>
<li>There are several expressions which will return inaccurate results due to cancellation. For example, subroutine <code>grs</code> computes the flattening using <em>f</em> = 1 − sqrt(1 − <em>e</em><sup>2</sup>). Much better is to use <em>f</em> = <em>e</em><sup>2</sup>/(1 + sqrt(1 − <em>e</em><sup>2</sup>)). The expressions for <em>q</em> and <em>q'</em> in <code>grs</code> and <code>radgrav</code> 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.</li>
<li>hsynth_WGS84 returns an "undefined" value for <em>xi</em> and <em>eta</em> at the poles. Better would be to return the value obtained by taking the limit <em>lat</em> → ± 90°.</li>
</ol>
<p>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 <a class="el" href="classGeographicLib_1_1GravityModel.html" title="Model of the earth's gravity field. ">GravityModel</a> corrects issues 3–5. On issue 2, <a class="el" href="classGeographicLib_1_1GravityModel.html" title="Model of the earth's gravity field. ">GravityModel</a> neglects the 1/<em>r</em> term in <em>T</em> in <a class="el" href="classGeographicLib_1_1GravityModel.html#a454bf53c090f2b85d8a3cfeb8dcc3961">GravityModel::GeoidHeight</a> and <a class="el" href="classGeographicLib_1_1GravityModel.html#a303c916a78467d6a210d8ba8be5c4686">GravityModel::SphericalAnomaly</a> in order to produce results which match NGA's for these quantities. On the other hand, <a class="el" href="classGeographicLib_1_1GravityModel.html#a2842adbe0574e1d34e8c9f4f74826320">GravityModel::Disturbance</a> and <a class="el" href="classGeographicLib_1_1GravityModel.html#a46e1863fe86c7a7dcae0bfb433eae287">GravityModel::T</a> <em>do</em> include this term.</p>
<h1><a class="anchor" id="gravitygeoid"></a>
Details of the geoid height and anomaly calculations</h1>
<p>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 <a class="el" href="classGeographicLib_1_1GravityModel.html" title="Model of the earth's gravity field. ">GravityModel</a> 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., <a class="el" href="classGeographicLib_1_1Geoid.html" title="Looking up the height of the geoid above the ellipsoid. ">Geoid</a>) and it is desirable to use a consistent definition of the geoid height.</p><ul>
<li>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.</li>
<li>The 1/<em>r</em> term is omitted from the <em>T</em> (this is issue 2 in <a class="el" href="gravity.html#gravitynga">Comments on the NGA harmonic synthesis code</a>). This moves the equipotential surface by about 5mm.</li>
<li>The surface <em>W</em> = <em>U</em><sub>0</sub> 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.</li>
<li>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.</li>
<li>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.)</li>
</ul>
<p>A useful discussion of the problems with defining a geoid is given by Dru A. Smith in <a href="https://www.ngs.noaa.gov/PUBS_LIB/EGM96_GEOID_PAPER/egm96_geoid_paper.html">There is no such thing as "The" EGM96 geoid: Subtle points on the use of a global geopotential model</a>, IGeS Bulletin No. 8, International <a class="el" href="classGeographicLib_1_1Geoid.html" title="Looking up the height of the geoid above the ellipsoid. ">Geoid</a> Service, Milan, Italy, pp. 17–28 (1998).</p>
<p><a class="el" href="classGeographicLib_1_1GravityModel.html#a454bf53c090f2b85d8a3cfeb8dcc3961">GravityModel::GeoidHeight</a> reproduces the results of the several NGA codes for harmonic synthesis with the following maximum discrepancies:</p><ul>
<li>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.</li>
<li>egm96 = 23nm.</li>
<li>egm2008 = 78pm. After addressing some of the issues alluded to in issue 4 in <a class="el" href="gravity.html#gravitynga">Comments on the NGA harmonic synthesis code</a>, the maximum discrepancy becomes 23pm.</li>
</ul>
<p>The formula for the gravity anomaly vector involves computing gravity and normal gravity at two different points (with the displacement between the points unknown <em>ab initio</em>). Since the gravity anomaly is already a small quantity it is sometimes acceptable to employ approximations that change the quantities by <em>O</em>(<em>f</em>). The NGA code uses the spherical approximation described by Heiskanen and Moritz, Sec. 2-14 and <a class="el" href="classGeographicLib_1_1GravityModel.html#a303c916a78467d6a210d8ba8be5c4686">GravityModel::SphericalAnomaly</a> uses the same approximation for compatibility. In this approximation, the gravity disturbance <b>delta</b> = <b>grad</b> <em>T</em> is calculated. Here, <em>T</em> once again excludes the 1/<em>r</em> term (this is issue 2 in <a class="el" href="gravity.html#gravitynga">Comments on the NGA harmonic synthesis code</a> and is consistent with the computation of the geoid height). Note that <b>delta</b> compares the gravity and the normal gravity at the <em>same</em> 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. <b>delta</b> is expressed in <em>spherical</em> coordinates as <em>deltax</em>, <em>deltay</em>, <em>deltaz</em> where, for example, <em>deltaz</em> is the <em>radial</em> component of <b>delta</b> (not the component perpendicular to the ellipsoid) and <em>deltay</em> is similarly slightly different from the usual northerly component. The components of the anomaly are then given by</p><ul>
<li>gravity anomaly, <em>Dg01</em> = <em>deltaz</em> − 2<em>T</em>/<em>R</em>, where <em>R</em> distance to the center of the earth;</li>
<li>northerly component of the deflection of the vertical, <em>xi</em> = − <em>deltay</em>/<em>gamma</em>, where <em>gamma</em> is the magnitude of the normal gravity;</li>
<li>easterly component of the deflection of the vertical, <em>eta</em> = − <em>deltax</em>/<em>gamma</em>.</li>
</ul>
<p><a class="el" href="classGeographicLib_1_1NormalGravity.html" title="The normal gravity of the earth. ">NormalGravity</a> computes the normal gravity accurately and avoids issue 3 of <a class="el" href="gravity.html#gravitynga">Comments on the NGA harmonic synthesis code</a>. Thus while <a class="el" href="classGeographicLib_1_1GravityModel.html#a303c916a78467d6a210d8ba8be5c4686">GravityModel::SphericalAnomaly</a> reproduces the results for <em>xi</em> and <em>eta</em> at <em>h</em> = 0, there is a slight discrepancy if <em>h</em> is non-zero.</p>
<h1><a class="anchor" id="gravityatmos"></a>
The effect of the mass of the atmosphere</h1>
<p>All of the supported models use WGS84 for the reference ellipsoid. This has (see <a href="http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf">TR8350.2</a>, table 3.1)</p><ul>
<li><em>a</em> = 6378137 m</li>
<li><em>f</em> = 1/298.257223563</li>
<li>ω = 7292115 × 10<sup>−11</sup> rad s<sup>−1</sup></li>
<li><em>GM</em> = 3986004.418 × 10<sup>8</sup> m<sup>3</sup> s<sup>−2</sup>.</li>
</ul>
<p>The value of <em>GM</em> 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 <em>g</em> will be less (in absolute value) than the value predicted by these models by about δ<em>g</em> = (4π<em>G</em>/<em>g</em>) <em>A</em> = 8.552 × 10<sup>−11</sup> <em>A m<sup>2</sup>/kg</em>, where <em>G</em> is the gravitational constant, <em>g</em> is the earth's gravity, and <em>A</em> is the pressure of the atmosphere. At sea level we have <em>A</em> = 101.3 kPa, and δ<em>g</em> = 8.7−×−10<sup>−6</sup> m s<sup>−2</sup>, approximately. (In other words the effect is about 1 part in a million; by way of comparison, buoyancy effects are about 100 times larger.)</p>
<h1><a class="anchor" id="gravityparallel"></a>
Geoid heights on a multi-processor system</h1>
<p>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:</p><ul>
<li>Use a <a class="el" href="classGeographicLib_1_1GravityCircle.html" title="Gravity on a circle of latitude. ">GravityCircle</a> 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).</li>
<li>Compute the values on several circles of latitude in parallel. One of the simplest ways of doing this is with <a href="http://openmp.org">OpenMP</a>; on an 8-processor system, this can speed up the computation by another factor of 8.</li>
</ul>
<p>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 <a href="https://vdatum.noaa.gov/docs/gtx_info.html#dev_gtx_binary">.gtx</a> 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!) </p><div class="fragment"><div class="line"><span class="comment">// Write out a gtx file of geoid heights above the ellipsoid. For egm2008 at</span></div><div class="line"><span class="comment">// 1' resolution this takes about 40 mins on a 8-processor Intel 2.66 GHz</span></div><div class="line"><span class="comment">// machine using OpenMP.</span></div><div class="line"><span class="comment">//</span></div><div class="line"><span class="comment">// For the format of gtx files, see</span></div><div class="line"><span class="comment">// http://vdatum.noaa.gov/dev/gtx_info.html#dev_gtx_binary</span></div><div class="line"><span class="comment">//</span></div><div class="line"><span class="comment">// data is binary big-endian:</span></div><div class="line"><span class="comment">// south latitude edge (degrees double)</span></div><div class="line"><span class="comment">// west longitude edge (degrees double)</span></div><div class="line"><span class="comment">// delta latitude (degrees double)</span></div><div class="line"><span class="comment">// delta longitude (degrees double)</span></div><div class="line"><span class="comment">// nlat = number of latitude rows (integer)</span></div><div class="line"><span class="comment">// nlong = number of longitude columns (integer)</span></div><div class="line"><span class="comment">// nlat * nlong geoid heights (meters float)</span></div><div class="line"></div><div class="line"><span class="preprocessor">#include <vector></span></div><div class="line"><span class="preprocessor">#include <iostream></span></div><div class="line"><span class="preprocessor">#include <fstream></span></div><div class="line"><span class="preprocessor">#include <string></span></div><div class="line"><span class="preprocessor">#include <algorithm></span></div><div class="line"></div><div class="line"><span class="preprocessor">#if defined(_OPENMP)</span></div><div class="line"><span class="preprocessor">#define HAVE_OPENMP 1</span></div><div class="line"><span class="preprocessor">#else</span></div><div class="line"><span class="preprocessor">#define HAVE_OPENMP 0</span></div><div class="line"><span class="preprocessor">#endif</span></div><div class="line"></div><div class="line"><span class="preprocessor">#if HAVE_OPENMP</span></div><div class="line"><span class="preprocessor"># include <omp.h></span></div><div class="line"><span class="preprocessor">#endif</span></div><div class="line"></div><div class="line"><span class="preprocessor">#include <<a class="code" href="GravityModel_8hpp.html">GeographicLib/GravityModel.hpp</a>></span></div><div class="line"><span class="preprocessor">#include <<a class="code" href="GravityCircle_8hpp.html">GeographicLib/GravityCircle.hpp</a>></span></div><div class="line"><span class="preprocessor">#include <<a class="code" href="Utility_8hpp.html">GeographicLib/Utility.hpp</a>></span></div><div class="line"></div><div class="line"><span class="keyword">using namespace </span><a class="code" href="namespacestd.html">std</a>;</div><div class="line"><span class="keyword">using namespace </span><a class="code" href="namespaceGeographicLib.html">GeographicLib</a>;</div><div class="line"></div><div class="line"><span class="keywordtype">int</span> <a class="code" href="CartConvert_8cpp.html#a814244b02f6701fd8ead246678bb3dae">main</a>(<span class="keywordtype">int</span> argc, <span class="keyword">const</span> <span class="keywordtype">char</span>* <span class="keyword">const</span> argv[]) {</div><div class="line"> <span class="comment">// Hardwired for 3 args:</span></div><div class="line"> <span class="comment">// 1 = the gravity model (e.g., egm2008)</span></div><div class="line"> <span class="comment">// 2 = intervals per degree</span></div><div class="line"> <span class="comment">// 3 = output GTX file</span></div><div class="line"> <span class="keywordflow">if</span> (argc != 4) {</div><div class="line"> cerr << <span class="stringliteral">"Usage: "</span> << argv[0]</div><div class="line"> << <span class="stringliteral">" gravity-model intervals-per-degree output.gtx\n"</span>;</div><div class="line"> <span class="keywordflow">return</span> 1;</div><div class="line"> }</div><div class="line"> <span class="keywordflow">try</span> {</div><div class="line"> <span class="comment">// Will need to set the precision for each thread, so save return value</span></div><div class="line"> <span class="keywordtype">int</span> ndigits = <a class="code" href="classGeographicLib_1_1Utility.html#ab95e26081d2c69019703adb8a4e4d3d0">Utility::set_digits</a>();</div><div class="line"> <span class="keywordtype">string</span> model(argv[1]);</div><div class="line"> <span class="comment">// Number of intervals per degree</span></div><div class="line"> <span class="keywordtype">int</span> ndeg = Utility::val<int>(string(argv[2]));</div><div class="line"> <span class="keywordtype">string</span> filename(argv[3]);</div><div class="line"> <a class="code" href="classGeographicLib_1_1GravityModel.html">GravityModel</a> g(model);</div><div class="line"> <span class="keywordtype">int</span></div><div class="line"> nlat = 180 * ndeg + 1,</div><div class="line"> nlon = 360 * ndeg;</div><div class="line"> <a class="code" href="classGeographicLib_1_1Math.html#aeee4778d7cf2f9fb9648efe4911da59d">Math::real</a></div><div class="line"> delta = 1 / <a class="code" href="classGeographicLib_1_1Math.html#aeee4778d7cf2f9fb9648efe4911da59d">Math::real</a>(ndeg), <span class="comment">// Grid spacing</span></div><div class="line"> latorg = -90,</div><div class="line"> lonorg = -180;</div><div class="line"> <span class="comment">// Write results as floats in binary mode</span></div><div class="line"> ofstream file(filename.c_str(), ios::binary);</div><div class="line"></div><div class="line"> <span class="comment">// Write header</span></div><div class="line"> {</div><div class="line"> <a class="code" href="classGeographicLib_1_1Math.html#aeee4778d7cf2f9fb9648efe4911da59d">Math::real</a> transform[] = {latorg, lonorg, delta, delta};</div><div class="line"> <span class="keywordtype">unsigned</span> sizes[] = {unsigned(nlat), unsigned(nlon)};</div><div class="line"> Utility::writearray<double, Math::real, true>(file, transform, 4);</div><div class="line"> Utility::writearray<unsigned, unsigned, true>(file, sizes, 2);</div><div class="line"> }</div><div class="line"></div><div class="line"> <span class="comment">// Compute and store results for nbatch latitudes at a time</span></div><div class="line"> <span class="keyword">const</span> <span class="keywordtype">int</span> nbatch = 64;</div><div class="line"> vector< vector<float> > N(nbatch, vector<float>(nlon));</div><div class="line"></div><div class="line"> <span class="keywordflow">for</span> (<span class="keywordtype">int</span> ilat0 = 0; ilat0 < nlat; ilat0 += nbatch) { <span class="comment">// Loop over batches</span></div><div class="line"> <span class="keywordtype">int</span> nlat0 = min(nlat, ilat0 + nbatch);</div><div class="line"></div><div class="line"><span class="preprocessor">#if HAVE_OPENMP</span></div><div class="line"><span class="preprocessor"># pragma omp parallel for</span></div><div class="line"><span class="preprocessor">#endif</span></div><div class="line"> <span class="keywordflow">for</span> (<span class="keywordtype">int</span> ilat = ilat0; ilat < nlat0; ++ilat) { <span class="comment">// Loop over latitudes</span></div><div class="line"> <a class="code" href="classGeographicLib_1_1Utility.html#ab95e26081d2c69019703adb8a4e4d3d0">Utility::set_digits</a>(ndigits); <span class="comment">// Set the precision</span></div><div class="line"> <a class="code" href="classGeographicLib_1_1Math.html#aeee4778d7cf2f9fb9648efe4911da59d">Math::real</a></div><div class="line"> lat = latorg + (ilat / ndeg) + delta * (ilat - ndeg * (ilat / ndeg)),</div><div class="line"> h = 0;</div><div class="line"> <a class="code" href="classGeographicLib_1_1GravityCircle.html">GravityCircle</a> c(g.Circle(lat, h, <a class="code" href="classGeographicLib_1_1GravityModel.html#af8691d0f13d6d42278cd1e615903d365a0ac768bffe3f104069c3fd5af1ddaa69">GravityModel::GEOID_HEIGHT</a>));</div><div class="line"> <span class="keywordflow">for</span> (<span class="keywordtype">int</span> ilon = 0; ilon < nlon; ++ilon) { <span class="comment">// Loop over longitudes</span></div><div class="line"> <a class="code" href="classGeographicLib_1_1Math.html#aeee4778d7cf2f9fb9648efe4911da59d">Math::real</a> lon = lonorg</div><div class="line"> + (ilon / ndeg) + delta * (ilon - ndeg * (ilon / ndeg));</div><div class="line"> N[ilat - ilat0][ilon] = float(c.GeoidHeight(lon));</div><div class="line"> } <span class="comment">// longitude loop</span></div><div class="line"> } <span class="comment">// latitude loop -- end of parallel section</span></div><div class="line"></div><div class="line"> <span class="keywordflow">for</span> (<span class="keywordtype">int</span> ilat = ilat0; ilat < nlat0; ++ilat) <span class="comment">// write out data</span></div><div class="line"> Utility::writearray<float, float, true>(file, N[ilat - ilat0]);</div><div class="line"> } <span class="comment">// batch loop</span></div><div class="line"> }</div><div class="line"> <span class="keywordflow">catch</span> (<span class="keyword">const</span> exception& e) {</div><div class="line"> cerr << <span class="stringliteral">"Caught exception: "</span> << e.what() << <span class="stringliteral">"\n"</span>;</div><div class="line"> <span class="keywordflow">return</span> 1;</div><div class="line"> }</div><div class="line"> <span class="keywordflow">catch</span> (...) {</div><div class="line"> cerr << <span class="stringliteral">"Caught unknown exception\n"</span>;</div><div class="line"> <span class="keywordflow">return</span> 1;</div><div class="line"> }</div><div class="line">}</div></div><!-- fragment --><p>cmake will add in support for OpenMP for <code>examples/GeoidToGTX.cpp</code>, if it is available.</p>
<center> Back to <a class="el" href="geoid.html">Geoid height</a>. Forward to <a class="el" href="normalgravity.html">Normal gravity</a>. Up to <a class="el" href="index.html#contents">Contents</a>. </center> </div></div><!-- contents -->
<!-- start footer part -->
<hr class="footer"/><address class="footer"><small>
Generated by  <a href="http://www.doxygen.org/index.html">
<img class="footer" src="doxygen.png" alt="doxygen"/>
</a> 1.8.13
</small></address>
</body>
</html>
|