GeographicLib
1.21
|
00001 /** 00002 * \file Geocentric.hpp 00003 * \brief Header for GeographicLib::Geocentric class 00004 * 00005 * Copyright (c) Charles Karney (2008-2011) <charles@karney.com> and licensed 00006 * under the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_GEOCENTRIC_HPP) 00011 #define GEOGRAPHICLIB_GEOCENTRIC_HPP \ 00012 "$Id: e9f709c85e61f60509c492429061cba04350eea8 $" 00013 00014 #include <vector> 00015 #include <algorithm> 00016 #include <GeographicLib/Constants.hpp> 00017 00018 namespace GeographicLib { 00019 00020 /** 00021 * \brief %Geocentric coordinates 00022 * 00023 * Convert between geodetic coordinates latitude = \e lat, longitude = \e 00024 * lon, height = \e h (measured vertically from the surface of the ellipsoid) 00025 * to geocentric coordinates (\e X, \e Y, \e Z). The origin of geocentric 00026 * coordinates is at the center of the earth. The \e Z axis goes thru the 00027 * north pole, \e lat = 90<sup>o</sup>. The \e X axis goes thru \e lat = 0, 00028 * \e lon = 0. %Geocentric coordinates are also known as earth centered, 00029 * earth fixed (ECEF) coordinates. 00030 * 00031 * The conversion from geographic to geocentric coordinates is 00032 * straightforward. For the reverse transformation we use 00033 * - H. Vermeille, 00034 * <a href="http://dx.doi.org/10.1007/s00190-002-0273-6"> Direct 00035 * transformation from geocentric coordinates to geodetic coordinates</a>, 00036 * J. Geodesy 76, 451–454 (2002). 00037 * . 00038 * Several changes have been made to ensure that the method returns accurate 00039 * results for all finite inputs (even if \e h is infinite). The changes are 00040 * described in Appendix B of 00041 * - C. F. F. Karney, 00042 * <a href="http://arxiv.org/abs/1102.1215v1">Geodesics 00043 * on an ellipsoid of revolution</a>, 00044 * Feb. 2011; 00045 * preprint 00046 * <a href="http://arxiv.org/abs/1102.1215v1">arxiv:1102.1215v1</a>. 00047 * . 00048 * See \ref geocentric for more information. 00049 * 00050 * The errors in these routines are close to round-off. Specifically, for 00051 * points within 5000 km of the surface of the ellipsoid (either inside or 00052 * outside the ellipsoid), the error is bounded by 7 nm (7 nanometers) for 00053 * the WGS84 ellipsoid. See \ref geocentric for further information on the 00054 * errors. 00055 * 00056 * Example of use: 00057 * \include example-Geocentric.cpp 00058 * 00059 * <a href="CartConvert.1.html">CartConvert</a> is a command-line utility 00060 * providing access to the functionality of Geocentric and LocalCartesian. 00061 **********************************************************************/ 00062 00063 class GEOGRAPHIC_EXPORT Geocentric { 00064 private: 00065 typedef Math::real real; 00066 friend class LocalCartesian; 00067 friend class MagneticCircle; // MagneticCircle uses Rotation 00068 friend class MagneticModel; // MagneticModel uses IntForward 00069 friend class GravityCircle; // GravityCircle uses Rotation 00070 friend class GravityModel; // GravityModel uses IntForward 00071 friend class NormalGravity; // NormalGravity uses IntForward 00072 friend class SphericalHarmonic; 00073 friend class SphericalHarmonic1; 00074 friend class SphericalHarmonic2; 00075 static const size_t dim_ = 3; 00076 static const size_t dim2_ = dim_ * dim_; 00077 real _a, _f, _e2, _e2m, _e2a, _e4a, _maxrad; 00078 static void Rotation(real sphi, real cphi, real slam, real clam, 00079 real M[dim2_]) throw(); 00080 static void Rotate(real M[dim2_], real x, real y, real z, 00081 real& X, real& Y, real& Z) throw() { 00082 // Perform [X,Y,Z]^t = M.[x,y,z]^t 00083 // (typically local cartesian to geocentric) 00084 X = M[0] * x + M[1] * y + M[2] * z; 00085 Y = M[3] * x + M[4] * y + M[5] * z; 00086 Z = M[6] * x + M[7] * y + M[8] * z; 00087 } 00088 static void Unrotate(real M[dim2_], real X, real Y, real Z, 00089 real& x, real& y, real& z) throw() { 00090 // Perform [x,y,z]^t = M^t.[X,Y,Z]^t 00091 // (typically geocentric to local cartesian) 00092 x = M[0] * X + M[3] * Y + M[6] * Z; 00093 y = M[1] * X + M[4] * Y + M[7] * Z; 00094 z = M[2] * X + M[5] * Y + M[8] * Z; 00095 } 00096 void IntForward(real lat, real lon, real h, real& X, real& Y, real& Z, 00097 real M[dim2_]) const throw(); 00098 void IntReverse(real X, real Y, real Z, real& lat, real& lon, real& h, 00099 real M[dim2_]) const throw(); 00100 00101 public: 00102 00103 /** 00104 * Constructor for a ellipsoid with 00105 * 00106 * @param[in] a equatorial radius (meters). 00107 * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. 00108 * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flattening 00109 * to 1/\e f. 00110 * 00111 * An exception is thrown if either of the axes of the ellipsoid is 00112 * non-positive. 00113 **********************************************************************/ 00114 Geocentric(real a, real f); 00115 00116 /** 00117 * A default constructor (for use by NormalGravity). 00118 **********************************************************************/ 00119 Geocentric() : _a(-1) {} 00120 00121 /** 00122 * Convert from geodetic to geocentric coordinates. 00123 * 00124 * @param[in] lat latitude of point (degrees). 00125 * @param[in] lon longitude of point (degrees). 00126 * @param[in] h height of point above the ellipsoid (meters). 00127 * @param[out] X geocentric coordinate (meters). 00128 * @param[out] Y geocentric coordinate (meters). 00129 * @param[out] Z geocentric coordinate (meters). 00130 * 00131 * \e lat should be in the range [-90, 90]; \e lon and \e lon0 should be in 00132 * the range [-180, 360]. 00133 **********************************************************************/ 00134 void Forward(real lat, real lon, real h, real& X, real& Y, real& Z) 00135 const throw() { 00136 if (Init()) 00137 IntForward(lat, lon, h, X, Y, Z, NULL); 00138 } 00139 00140 /** 00141 * Convert from geodetic to geocentric coordinates and return rotation 00142 * matrix. 00143 * 00144 * @param[in] lat latitude of point (degrees). 00145 * @param[in] lon longitude of point (degrees). 00146 * @param[in] h height of point above the ellipsoid (meters). 00147 * @param[out] X geocentric coordinate (meters). 00148 * @param[out] Y geocentric coordinate (meters). 00149 * @param[out] Z geocentric coordinate (meters). 00150 * @param[out] M if the length of the vector is 9, fill with the rotation 00151 * matrix in row-major order. 00152 * 00153 * Let \e v be a unit vector located at (\e lat, \e lon, \e h). We can 00154 * express \e v as \e column vectors in one of two ways 00155 * - in east, north, up coordinates (where the components are relative to a 00156 * local coordinate system at (\e lat, \e lon, \e h)); call this 00157 * representation \e v1. 00158 * - in geocentric \e X, \e Y, \e Z coordinates; call this representation 00159 * \e v0. 00160 * . 00161 * Then we have \e v0 = \e M . \e v1. 00162 **********************************************************************/ 00163 void Forward(real lat, real lon, real h, real& X, real& Y, real& Z, 00164 std::vector<real>& M) 00165 const throw() { 00166 if (!Init()) 00167 return; 00168 if (M.end() == M.begin() + dim2_) { 00169 real t[dim2_]; 00170 IntForward(lat, lon, h, X, Y, Z, t); 00171 copy(t, t + dim2_, M.begin()); 00172 } else 00173 IntForward(lat, lon, h, X, Y, Z, NULL); 00174 } 00175 00176 /** 00177 * Convert from geocentric to geodetic to coordinates. 00178 * 00179 * @param[in] X geocentric coordinate (meters). 00180 * @param[in] Y geocentric coordinate (meters). 00181 * @param[in] Z geocentric coordinate (meters). 00182 * @param[out] lat latitude of point (degrees). 00183 * @param[out] lon longitude of point (degrees). 00184 * @param[out] h height of point above the ellipsoid (meters). 00185 * 00186 * In general there are multiple solutions and the result which maximizes 00187 * \e h is returned. If there are still multiple solutions with different 00188 * latitudes (applies only if \e Z = 0), then the solution with \e lat > 0 00189 * is returned. If there are still multiple solutions with different 00190 * longitudes (applies only if \e X = \e Y = 0) then \e lon = 0 is 00191 * returned. The value of \e h returned satisfies \e h >= - \e a (1 - 00192 * <i>e</i><sup>2</sup>) / sqrt(1 - <i>e</i><sup>2</sup> sin<sup>2</sup>\e 00193 * lat). The value of \e lon returned is in the range [-180, 180). 00194 **********************************************************************/ 00195 void Reverse(real X, real Y, real Z, real& lat, real& lon, real& h) 00196 const throw() { 00197 if (Init()) 00198 IntReverse(X, Y, Z, lat, lon, h, NULL); 00199 } 00200 00201 /** 00202 * Convert from geocentric to geodetic to coordinates. 00203 * 00204 * @param[in] X geocentric coordinate (meters). 00205 * @param[in] Y geocentric coordinate (meters). 00206 * @param[in] Z geocentric coordinate (meters). 00207 * @param[out] lat latitude of point (degrees). 00208 * @param[out] lon longitude of point (degrees). 00209 * @param[out] h height of point above the ellipsoid (meters). 00210 * @param[out] M if the length of the vector is 9, fill with the rotation 00211 * matrix in row-major order. 00212 * 00213 * Let \e v be a unit vector located at (\e lat, \e lon, \e h). We can 00214 * express \e v as \e column vectors in one of two ways 00215 * - in east, north, up coordinates (where the components are relative to a 00216 * local coordinate system at (\e lat, \e lon, \e h)); call this 00217 * representation \e v1. 00218 * - in geocentric \e X, \e Y, \e Z coordinates; call this representation 00219 * \e v0. 00220 * . 00221 * Then we have \e v1 = \e M^T . \e v0, where \e M^T is the transpose of \e 00222 * M. 00223 **********************************************************************/ 00224 void Reverse(real X, real Y, real Z, real& lat, real& lon, real& h, 00225 std::vector<real>& M) 00226 const throw() { 00227 if (!Init()) 00228 return; 00229 if (M.end() == M.begin() + dim2_) { 00230 real t[dim2_]; 00231 IntReverse(X, Y, Z, lat, lon, h, t); 00232 copy(t, t + dim2_, M.begin()); 00233 } else 00234 IntReverse(X, Y, Z, lat, lon, h, NULL); 00235 } 00236 00237 /** \name Inspector functions 00238 **********************************************************************/ 00239 ///@{ 00240 /** 00241 * @return true if the object has been initialized. 00242 **********************************************************************/ 00243 bool Init() const throw() { return _a > 0; } 00244 /** 00245 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00246 * the value used in the constructor. 00247 **********************************************************************/ 00248 Math::real MajorRadius() const throw() 00249 { return Init() ? _a : Math::NaN<real>(); } 00250 00251 /** 00252 * @return \e f the flattening of the ellipsoid. This is the 00253 * value used in the constructor. 00254 **********************************************************************/ 00255 Math::real Flattening() const throw() 00256 { return Init() ? _f : Math::NaN<real>(); } 00257 ///@} 00258 00259 /// \cond SKIP 00260 /** 00261 * <b>DEPRECATED</b> 00262 * @return \e r the inverse flattening of the ellipsoid. 00263 **********************************************************************/ 00264 Math::real InverseFlattening() const throw() 00265 { return Init() ? 1/_f : Math::NaN<real>(); } 00266 /// \endcond 00267 00268 /** 00269 * A global instantiation of Geocentric with the parameters for the WGS84 00270 * ellipsoid. 00271 **********************************************************************/ 00272 static const Geocentric WGS84; 00273 }; 00274 00275 } // namespace GeographicLib 00276 00277 #endif // GEOGRAPHICLIB_GEOCENTRIC_HPP