/** * \file Geocentric.hpp * \brief Header for GeographicLib::Geocentric class * * Copyright (c) Charles Karney (2008-2020) and licensed * under the MIT/X11 License. For more information, see * https://geographiclib.sourceforge.io/ **********************************************************************/ #if !defined(GEOGRAPHICLIB_GEOCENTRIC_HPP) #define GEOGRAPHICLIB_GEOCENTRIC_HPP 1 #include #include namespace GeographicLib { /** * \brief %Geocentric coordinates * * Convert between geodetic coordinates latitude = \e lat, longitude = \e * lon, height = \e h (measured vertically from the surface of the ellipsoid) * to geocentric coordinates (\e X, \e Y, \e Z). The origin of geocentric * coordinates is at the center of the earth. The \e Z axis goes thru the * north pole, \e lat = 90°. The \e X axis goes thru \e lat = 0, * \e lon = 0. %Geocentric coordinates are also known as earth centered, * earth fixed (ECEF) coordinates. * * The conversion from geographic to geocentric coordinates is * straightforward. For the reverse transformation we use * - H. Vermeille, * Direct * transformation from geocentric coordinates to geodetic coordinates, * J. Geodesy 76, 451--454 (2002). * . * Several changes have been made to ensure that the method returns accurate * results for all finite inputs (even if \e h is infinite). The changes are * described in Appendix B of * - C. F. F. Karney, * Geodesics * on an ellipsoid of revolution, * Feb. 2011; * preprint * arxiv:1102.1215v1. * . * Vermeille similarly updated his method in * - H. Vermeille, * * An analytical method to transform geocentric into * geodetic coordinates, J. Geodesy 85, 105--117 (2011). * . * See \ref geocentric for more information. * * The errors in these routines are close to round-off. Specifically, for * points within 5000 km of the surface of the ellipsoid (either inside or * outside the ellipsoid), the error is bounded by 7 nm (7 nanometers) for * the WGS84 ellipsoid. See \ref geocentric for further information on the * errors. * * Example of use: * \include example-Geocentric.cpp * * CartConvert is a command-line utility * providing access to the functionality of Geocentric and LocalCartesian. **********************************************************************/ class GEOGRAPHICLIB_EXPORT Geocentric { private: typedef Math::real real; friend class LocalCartesian; friend class MagneticCircle; // MagneticCircle uses Rotation friend class MagneticModel; // MagneticModel uses IntForward friend class GravityCircle; // GravityCircle uses Rotation friend class GravityModel; // GravityModel uses IntForward friend class NormalGravity; // NormalGravity uses IntForward static const size_t dim_ = 3; static const size_t dim2_ = dim_ * dim_; real _a, _f, _e2, _e2m, _e2a, _e4a, _maxrad; static void Rotation(real sphi, real cphi, real slam, real clam, real M[dim2_]); static void Rotate(const real M[dim2_], real x, real y, real z, real& X, real& Y, real& Z) { // Perform [X,Y,Z]^t = M.[x,y,z]^t // (typically local cartesian to geocentric) X = M[0] * x + M[1] * y + M[2] * z; Y = M[3] * x + M[4] * y + M[5] * z; Z = M[6] * x + M[7] * y + M[8] * z; } static void Unrotate(const real M[dim2_], real X, real Y, real Z, real& x, real& y, real& z) { // Perform [x,y,z]^t = M^t.[X,Y,Z]^t // (typically geocentric to local cartesian) x = M[0] * X + M[3] * Y + M[6] * Z; y = M[1] * X + M[4] * Y + M[7] * Z; z = M[2] * X + M[5] * Y + M[8] * Z; } void IntForward(real lat, real lon, real h, real& X, real& Y, real& Z, real M[dim2_]) const; void IntReverse(real X, real Y, real Z, real& lat, real& lon, real& h, real M[dim2_]) const; public: /** * Constructor for a ellipsoid with * * @param[in] a equatorial radius (meters). * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere. * Negative \e f gives a prolate ellipsoid. * @exception GeographicErr if \e a or (1 − \e f) \e a is not * positive. **********************************************************************/ Geocentric(real a, real f); /** * A default constructor (for use by NormalGravity). **********************************************************************/ Geocentric() : _a(-1) {} /** * Convert from geodetic to geocentric coordinates. * * @param[in] lat latitude of point (degrees). * @param[in] lon longitude of point (degrees). * @param[in] h height of point above the ellipsoid (meters). * @param[out] X geocentric coordinate (meters). * @param[out] Y geocentric coordinate (meters). * @param[out] Z geocentric coordinate (meters). * * \e lat should be in the range [−90°, 90°]. **********************************************************************/ void Forward(real lat, real lon, real h, real& X, real& Y, real& Z) const { if (Init()) IntForward(lat, lon, h, X, Y, Z, NULL); } /** * Convert from geodetic to geocentric coordinates and return rotation * matrix. * * @param[in] lat latitude of point (degrees). * @param[in] lon longitude of point (degrees). * @param[in] h height of point above the ellipsoid (meters). * @param[out] X geocentric coordinate (meters). * @param[out] Y geocentric coordinate (meters). * @param[out] Z geocentric coordinate (meters). * @param[out] M if the length of the vector is 9, fill with the rotation * matrix in row-major order. * * Let \e v be a unit vector located at (\e lat, \e lon, \e h). We can * express \e v as \e column vectors in one of two ways * - in east, north, up coordinates (where the components are relative to a * local coordinate system at (\e lat, \e lon, \e h)); call this * representation \e v1. * - in geocentric \e X, \e Y, \e Z coordinates; call this representation * \e v0. * . * Then we have \e v0 = \e M ⋅ \e v1. **********************************************************************/ void Forward(real lat, real lon, real h, real& X, real& Y, real& Z, std::vector& M) const { if (!Init()) return; if (M.end() == M.begin() + dim2_) { real t[dim2_]; IntForward(lat, lon, h, X, Y, Z, t); std::copy(t, t + dim2_, M.begin()); } else IntForward(lat, lon, h, X, Y, Z, NULL); } /** * Convert from geocentric to geodetic to coordinates. * * @param[in] X geocentric coordinate (meters). * @param[in] Y geocentric coordinate (meters). * @param[in] Z geocentric coordinate (meters). * @param[out] lat latitude of point (degrees). * @param[out] lon longitude of point (degrees). * @param[out] h height of point above the ellipsoid (meters). * * In general, there are multiple solutions and the result which minimizes * |h |is returned, i.e., (lat, lon) corresponds to * the closest point on the ellipsoid. If there are still multiple * solutions with different latitudes (applies only if \e Z = 0), then the * solution with \e lat > 0 is returned. If there are still multiple * solutions with different longitudes (applies only if \e X = \e Y = 0) * then \e lon = 0 is returned. The value of \e h returned satisfies \e h * ≥ − \e a (1 − e2) / sqrt(1 − * e2 sin2\e lat). The value of \e lon * returned is in the range [−180°, 180°]. **********************************************************************/ void Reverse(real X, real Y, real Z, real& lat, real& lon, real& h) const { if (Init()) IntReverse(X, Y, Z, lat, lon, h, NULL); } /** * Convert from geocentric to geodetic to coordinates. * * @param[in] X geocentric coordinate (meters). * @param[in] Y geocentric coordinate (meters). * @param[in] Z geocentric coordinate (meters). * @param[out] lat latitude of point (degrees). * @param[out] lon longitude of point (degrees). * @param[out] h height of point above the ellipsoid (meters). * @param[out] M if the length of the vector is 9, fill with the rotation * matrix in row-major order. * * Let \e v be a unit vector located at (\e lat, \e lon, \e h). We can * express \e v as \e column vectors in one of two ways * - in east, north, up coordinates (where the components are relative to a * local coordinate system at (\e lat, \e lon, \e h)); call this * representation \e v1. * - in geocentric \e X, \e Y, \e Z coordinates; call this representation * \e v0. * . * Then we have \e v1 = MT ⋅ \e v0, where * MT is the transpose of \e M. **********************************************************************/ void Reverse(real X, real Y, real Z, real& lat, real& lon, real& h, std::vector& M) const { if (!Init()) return; if (M.end() == M.begin() + dim2_) { real t[dim2_]; IntReverse(X, Y, Z, lat, lon, h, t); std::copy(t, t + dim2_, M.begin()); } else IntReverse(X, Y, Z, lat, lon, h, NULL); } /** \name Inspector functions **********************************************************************/ ///@{ /** * @return true if the object has been initialized. **********************************************************************/ bool Init() const { return _a > 0; } /** * @return \e a the equatorial radius of the ellipsoid (meters). This is * the value used in the constructor. **********************************************************************/ Math::real EquatorialRadius() const { return Init() ? _a : Math::NaN(); } /** * @return \e f the flattening of the ellipsoid. This is the * value used in the constructor. **********************************************************************/ Math::real Flattening() const { return Init() ? _f : Math::NaN(); } /** * \deprecated An old name for EquatorialRadius(). **********************************************************************/ GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()") Math::real MajorRadius() const { return EquatorialRadius(); } ///@} /** * A global instantiation of Geocentric with the parameters for the WGS84 * ellipsoid. **********************************************************************/ static const Geocentric& WGS84(); }; } // namespace GeographicLib #endif // GEOGRAPHICLIB_GEOCENTRIC_HPP