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