/usr/share/doc/geographiclib/html/geocentric.html is in geographiclib-doc 1.45-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 | <!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.9.1"/>
<title>GeographicLib: Geocentric coordinates</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 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 style="padding-left: 0.5em;">
<div id="projectname">GeographicLib
 <span id="projectnumber">1.45</span>
</div>
</td>
</tr>
</tbody>
</table>
</div>
<!-- end header part -->
<!-- Generated by Doxygen 1.8.9.1 -->
<div id="navrow1" class="tabs">
<ul class="tablist">
<li><a href="index.html"><span>Main Page</span></a></li>
<li class="current"><a href="pages.html"><span>Related Pages</span></a></li>
<li><a href="namespaces.html"><span>Namespaces</span></a></li>
<li><a href="annotated.html"><span>Classes</span></a></li>
<li><a href="files.html"><span>Files</span></a></li>
</ul>
</div>
</div><!-- top -->
<div class="header">
<div class="headertitle">
<div class="title"><a class="el" href="classGeographicLib_1_1Geocentric.html" title="Geocentric coordinates ">Geocentric</a> coordinates </div> </div>
</div><!--header-->
<div class="contents">
<div class="textblock"><center> Back to <a class="el" href="transversemercator.html">Transverse Mercator projection</a>. Forward to <a class="el" href="auxlat.html">Auxiliary latitudes</a>. Up to <a class="el" href="index.html#contents">Contents</a>. </center><p>The implementation of <a class="el" href="classGeographicLib_1_1Geocentric.html#a1907735ce8f5f915a14a5f7a8b3adfea">Geocentric::Reverse</a> is adapted from</p><ul>
<li>H. Vermeille, <a href="https://dx.doi.org/10.1007/s00190-002-0273-6">Direct transformation from geocentric coordinates to geodetic coordinates</a>, J. Geodesy 76, 451–454 (2002).</li>
</ul>
<p>This provides a closed-form solution but can't directly be applied close to the center of the earth. Several changes have been made to remove this restriction and to improve the numerical accuracy. Now the method is accurate for all inputs (even if <em>h</em> is infinite). The changes are described in Appendix B of</p><ul>
<li>C. F. F. Karney, <a href="http://arxiv.org/abs/1102.1215v1">Geodesics on an ellipsoid of revolution</a>, Feb. 2011; preprint <a href="http://arxiv.org/abs/1102.1215v1">arxiv:1102.1215v1</a>.</li>
</ul>
<p>Vermeille similarly updated his method in</p><ul>
<li>H. Vermeille, <a href="https://dx.doi.org/10.1007/s00190-010-0419-x">An analytical method to transform geocentric into geodetic coordinates</a>, J. Geodesy 85, 105–117 (2011).</li>
</ul>
<p>The problems encountered near the center of the ellipsoid are:</p><ul>
<li>There's a potential division by zero in the definition of <em>s</em>. The equations are easily reformulated to avoid this problem.</li>
<li><em>t</em><sup>3</sup> may be negative. This is OK; we just take the real root.</li>
<li>The solution for <em>t</em> may be complex. However this leads to 3 real roots for <em>u</em>/<em>r</em>. It's then just a matter of picking the one that computes the geodetic result which minimizes |<em>h</em>| and which avoids large round-off errors.</li>
<li>Some of the equations result in a large loss of accuracy due to subtracting nearly equal quantities. E.g., <em>k</em> = sqrt(<em>u</em> + <em>v</em> + <em>w</em><sup>2</sup>) − <em>w</em> is inaccurate if <em>u</em> + <em>v</em> is small; we can fix this by writing <em>k</em> = (<em>u</em> + <em>v</em>)/(sqrt(<em>u</em> + <em>v</em> + <em>w</em><sup>2</sup>) + <em>w</em>).</li>
</ul>
<p>The error is computed as follows. Write a version of Geocentric::WGS84.Forward which uses long doubles (including using long doubles for the WGS84 parameters). Generate random (long double) geodetic coordinates (<em>lat0</em>, <em>lon0</em>, <em>h0</em>) and use the "long double" WGS84.Forward to obtain the corresponding (long double) geocentric coordinates (<em>x0</em>, <em>y0</em>, <em>z0</em>). [We restrict <em>h0</em> so that <em>h0</em> ≥ − <em>a</em> (1 − <em>e</em><sup>2</sup>) / sqrt(1 − <em>e</em><sup>2</sup> sin<sup>2</sup><em>lat0</em>), which ensures that (<em>lat0</em>, <em>lon0</em>, <em>h0</em>) is the principal geodetic inverse of (<em>x0</em>, <em>y0</em>, <em>z0</em>).] Because the forward calculation is numerically stable and because long doubles (on Linux systems using g++) provide 11 bits additional accuracy (about 3.3 decimal digits), we regard this set of test data as exact.</p>
<p>Apply the double version of WGS84.Reverse to (<em>x0</em>, <em>y0</em>, <em>z0</em>) to compute the approximate geodetic coordinates (<em>lat1</em>, <em>lon1</em>, <em>h1</em>). Convert (<em>lat1</em> − <em>lat0</em>, <em>lon1</em> − <em>lon0</em>) to a distance, <em>ds</em>, on the surface of the ellipsoid and define <em>err</em> = hypot(<em>ds</em>, <em>h1</em> − <em>h0</em>). For |<em>h0</em>| < 5000 km, we have <em>err</em> < 7 nm (7 nanometers).</p>
<p>This methodology is not very useful very far from the globe, because the absolute errors in the approximate geodetic height become large, or within 50 km of the center of the earth, because of errors in computing the approximate geodetic latitude. To illustrate the second issue, the maximum value of <em>err</em> for <em>h0</em> < 0 is about 80 mm. The error is maximum close to the circle given by geocentric coordinates satisfying hypot(<em>x</em>, <em>y</em>) = <em>a</em> <em>e</em><sup>2</sup> (= 42.7 km), <em>z</em> = 0. (This is the center of meridional curvature for <em>lat</em> = 0.) The geodetic latitude for these points is <em>lat</em> = 0. However, if we move 1 nm towards the center of the earth, the geodetic latitude becomes 0.04", a distance of 1.4 m from the equator. If, instead, we move 1 nm up, the geodetic latitude becomes 7.45", a distance of 229 m from the equator. In light of this, Reverse does quite well in this vicinity.</p>
<p>To obtain a practical measure of the error for the general case we define</p><ul>
<li><em>err</em><sub>h</sub> = |<em>h1</em> − <em>h0</em>| / max(1, <em>h0</em> / <em>a</em>)</li>
<li>for <em>h0</em> > 0, <em>err</em><sub>out</sub> = <em>ds</em> </li>
<li>for <em>h0</em> < 0, apply the long double version of WGS84.Forward to (<em>lat1</em>, <em>lon1</em>, <em>h1</em>) to give (<em>x1</em>, <em>y1</em>, <em>z1</em>) and compute <em>err</em><sub>in</sub> = hypot(<em>x1</em> − <em>x0</em>, <em>y1</em> − <em>y0</em>, <em>z1</em> − <em>z0</em>).</li>
</ul>
<p>We then find <em>err</em><sub>h</sub> < 8 nm, <em>err</em><sub>out</sub> < 4 nm, and <em>err</em><sub>in</sub> < 7 nm. (1 nm = 1 nanometer.)</p>
<p>The testing has been confined to the WGS84 ellipsoid. The method will work for all ellipsoids used in terrestrial geodesy. However, the central region, which leads to multiple real roots for the cubic equation in Reverse, pokes outside the ellipsoid (at the poles) for ellipsoids with <em>e</em> > 1/sqrt(2); Reverse has not been analyzed for this case. Similarly ellipsoids which are very nearly spherical may yield inaccurate results due to underflow; in the other hand, the case of the sphere, <em>f</em> = 0, is treated specially and gives accurate results.</p>
<p>Other comparable methods are K. M. Borkowski, <a href="https://dx.doi.org/10.1007/BF00643807">Transformation of geocentric to geodetic coordinates without approximations</a>, Astrophys. Space Sci. 139, 1–4 (1987) (<a href="https://dx.doi.org/10.1007/BF00656995">erratum</a>) and T. Fukushima, <a href="https://dx.doi.org/10.1007/s001900050271">Fast transform from geocentric to geodetic coordinates</a>, J. Geodesy 73, 603–610 (1999). However the choice of independent variables in these methods leads to a loss of accuracy for points near the equatorial plane.</p>
<center> Back to <a class="el" href="transversemercator.html">Transverse Mercator projection</a>. Forward to <a class="el" href="auxlat.html">Auxiliary latitudes</a>. Up to <a class="el" href="index.html#contents">Contents</a>.</center><center></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.9.1
</small></address>
</body>
</html>
|