GeographicLib  1.50.1
Geocentric.hpp
Go to the documentation of this file.
1 /**
2  * \file Geocentric.hpp
3  * \brief Header for GeographicLib::Geocentric class
4  *
5  * Copyright (c) Charles Karney (2008-2019) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
8  **********************************************************************/
9 
10 #if !defined(GEOGRAPHICLIB_GEOCENTRIC_HPP)
11 #define GEOGRAPHICLIB_GEOCENTRIC_HPP 1
12 
13 #include <vector>
15 
16 namespace GeographicLib {
17 
18  /**
19  * \brief %Geocentric coordinates
20  *
21  * Convert between geodetic coordinates latitude = \e lat, longitude = \e
22  * lon, height = \e h (measured vertically from the surface of the ellipsoid)
23  * to geocentric coordinates (\e X, \e Y, \e Z). The origin of geocentric
24  * coordinates is at the center of the earth. The \e Z axis goes thru the
25  * north pole, \e lat = 90&deg;. The \e X axis goes thru \e lat = 0,
26  * \e lon = 0. %Geocentric coordinates are also known as earth centered,
27  * earth fixed (ECEF) coordinates.
28  *
29  * The conversion from geographic to geocentric coordinates is
30  * straightforward. For the reverse transformation we use
31  * - H. Vermeille,
32  * <a href="https://doi.org/10.1007/s00190-002-0273-6"> Direct
33  * transformation from geocentric coordinates to geodetic coordinates</a>,
34  * J. Geodesy 76, 451--454 (2002).
35  * .
36  * Several changes have been made to ensure that the method returns accurate
37  * results for all finite inputs (even if \e h is infinite). The changes are
38  * described in Appendix B of
39  * - C. F. F. Karney,
40  * <a href="https://arxiv.org/abs/1102.1215v1">Geodesics
41  * on an ellipsoid of revolution</a>,
42  * Feb. 2011;
43  * preprint
44  * <a href="https://arxiv.org/abs/1102.1215v1">arxiv:1102.1215v1</a>.
45  * .
46  * Vermeille similarly updated his method in
47  * - H. Vermeille,
48  * <a href="https://doi.org/10.1007/s00190-010-0419-x">
49  * An analytical method to transform geocentric into
50  * geodetic coordinates</a>, J. Geodesy 85, 105--117 (2011).
51  * .
52  * See \ref geocentric for more information.
53  *
54  * The errors in these routines are close to round-off. Specifically, for
55  * points within 5000 km of the surface of the ellipsoid (either inside or
56  * outside the ellipsoid), the error is bounded by 7 nm (7 nanometers) for
57  * the WGS84 ellipsoid. See \ref geocentric for further information on the
58  * errors.
59  *
60  * Example of use:
61  * \include example-Geocentric.cpp
62  *
63  * <a href="CartConvert.1.html">CartConvert</a> is a command-line utility
64  * providing access to the functionality of Geocentric and LocalCartesian.
65  **********************************************************************/
66 
68  private:
69  typedef Math::real real;
70  friend class LocalCartesian;
71  friend class MagneticCircle; // MagneticCircle uses Rotation
72  friend class MagneticModel; // MagneticModel uses IntForward
73  friend class GravityCircle; // GravityCircle uses Rotation
74  friend class GravityModel; // GravityModel uses IntForward
75  friend class NormalGravity; // NormalGravity uses IntForward
76  static const size_t dim_ = 3;
77  static const size_t dim2_ = dim_ * dim_;
78  real _a, _f, _e2, _e2m, _e2a, _e4a, _maxrad;
79  static void Rotation(real sphi, real cphi, real slam, real clam,
80  real M[dim2_]);
81  static void Rotate(real M[dim2_], real x, real y, real z,
82  real& X, real& Y, real& Z) {
83  // Perform [X,Y,Z]^t = M.[x,y,z]^t
84  // (typically local cartesian to geocentric)
85  X = M[0] * x + M[1] * y + M[2] * z;
86  Y = M[3] * x + M[4] * y + M[5] * z;
87  Z = M[6] * x + M[7] * y + M[8] * z;
88  }
89  static void Unrotate(real M[dim2_], real X, real Y, real Z,
90  real& x, real& y, real& z) {
91  // Perform [x,y,z]^t = M^t.[X,Y,Z]^t
92  // (typically geocentric to local cartesian)
93  x = M[0] * X + M[3] * Y + M[6] * Z;
94  y = M[1] * X + M[4] * Y + M[7] * Z;
95  z = M[2] * X + M[5] * Y + M[8] * Z;
96  }
97  void IntForward(real lat, real lon, real h, real& X, real& Y, real& Z,
98  real M[dim2_]) const;
99  void IntReverse(real X, real Y, real Z, real& lat, real& lon, real& h,
100  real M[dim2_]) const;
101 
102  public:
103 
104  /**
105  * Constructor for a ellipsoid with
106  *
107  * @param[in] a equatorial radius (meters).
108  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
109  * Negative \e f gives a prolate ellipsoid.
110  * @exception GeographicErr if \e a or (1 &minus; \e f) \e a is not
111  * positive.
112  **********************************************************************/
113  Geocentric(real a, real f);
114 
115  /**
116  * A default constructor (for use by NormalGravity).
117  **********************************************************************/
118  Geocentric() : _a(-1) {}
119 
120  /**
121  * Convert from geodetic to geocentric coordinates.
122  *
123  * @param[in] lat latitude of point (degrees).
124  * @param[in] lon longitude of point (degrees).
125  * @param[in] h height of point above the ellipsoid (meters).
126  * @param[out] X geocentric coordinate (meters).
127  * @param[out] Y geocentric coordinate (meters).
128  * @param[out] Z geocentric coordinate (meters).
129  *
130  * \e lat should be in the range [&minus;90&deg;, 90&deg;].
131  **********************************************************************/
132  void Forward(real lat, real lon, real h, real& X, real& Y, real& Z)
133  const {
134  if (Init())
135  IntForward(lat, lon, h, X, Y, Z, NULL);
136  }
137 
138  /**
139  * Convert from geodetic to geocentric coordinates and return rotation
140  * matrix.
141  *
142  * @param[in] lat latitude of point (degrees).
143  * @param[in] lon longitude of point (degrees).
144  * @param[in] h height of point above the ellipsoid (meters).
145  * @param[out] X geocentric coordinate (meters).
146  * @param[out] Y geocentric coordinate (meters).
147  * @param[out] Z geocentric coordinate (meters).
148  * @param[out] M if the length of the vector is 9, fill with the rotation
149  * matrix in row-major order.
150  *
151  * Let \e v be a unit vector located at (\e lat, \e lon, \e h). We can
152  * express \e v as \e column vectors in one of two ways
153  * - in east, north, up coordinates (where the components are relative to a
154  * local coordinate system at (\e lat, \e lon, \e h)); call this
155  * representation \e v1.
156  * - in geocentric \e X, \e Y, \e Z coordinates; call this representation
157  * \e v0.
158  * .
159  * Then we have \e v0 = \e M &sdot; \e v1.
160  **********************************************************************/
161  void Forward(real lat, real lon, real h, real& X, real& Y, real& Z,
162  std::vector<real>& M)
163  const {
164  if (!Init())
165  return;
166  if (M.end() == M.begin() + dim2_) {
167  real t[dim2_];
168  IntForward(lat, lon, h, X, Y, Z, t);
169  std::copy(t, t + dim2_, M.begin());
170  } else
171  IntForward(lat, lon, h, X, Y, Z, NULL);
172  }
173 
174  /**
175  * Convert from geocentric to geodetic to coordinates.
176  *
177  * @param[in] X geocentric coordinate (meters).
178  * @param[in] Y geocentric coordinate (meters).
179  * @param[in] Z geocentric coordinate (meters).
180  * @param[out] lat latitude of point (degrees).
181  * @param[out] lon longitude of point (degrees).
182  * @param[out] h height of point above the ellipsoid (meters).
183  *
184  * In general, there are multiple solutions and the result which minimizes
185  * |<i>h</i> |is returned, i.e., (<i>lat</i>, <i>lon</i>) corresponds to
186  * the closest point on the ellipsoid. If there are still multiple
187  * solutions with different latitudes (applies only if \e Z = 0), then the
188  * solution with \e lat > 0 is returned. If there are still multiple
189  * solutions with different longitudes (applies only if \e X = \e Y = 0)
190  * then \e lon = 0 is returned. The value of \e h returned satisfies \e h
191  * &ge; &minus; \e a (1 &minus; <i>e</i><sup>2</sup>) / sqrt(1 &minus;
192  * <i>e</i><sup>2</sup> sin<sup>2</sup>\e lat). The value of \e lon
193  * returned is in the range [&minus;180&deg;, 180&deg;].
194  **********************************************************************/
195  void Reverse(real X, real Y, real Z, real& lat, real& lon, real& h)
196  const {
197  if (Init())
198  IntReverse(X, Y, Z, lat, lon, h, NULL);
199  }
200 
201  /**
202  * Convert from geocentric to geodetic to coordinates.
203  *
204  * @param[in] X geocentric coordinate (meters).
205  * @param[in] Y geocentric coordinate (meters).
206  * @param[in] Z geocentric coordinate (meters).
207  * @param[out] lat latitude of point (degrees).
208  * @param[out] lon longitude of point (degrees).
209  * @param[out] h height of point above the ellipsoid (meters).
210  * @param[out] M if the length of the vector is 9, fill with the rotation
211  * matrix in row-major order.
212  *
213  * Let \e v be a unit vector located at (\e lat, \e lon, \e h). We can
214  * express \e v as \e column vectors in one of two ways
215  * - in east, north, up coordinates (where the components are relative to a
216  * local coordinate system at (\e lat, \e lon, \e h)); call this
217  * representation \e v1.
218  * - in geocentric \e X, \e Y, \e Z coordinates; call this representation
219  * \e v0.
220  * .
221  * Then we have \e v1 = <i>M</i><sup>T</sup> &sdot; \e v0, where
222  * <i>M</i><sup>T</sup> is the transpose of \e M.
223  **********************************************************************/
224  void Reverse(real X, real Y, real Z, real& lat, real& lon, real& h,
225  std::vector<real>& M)
226  const {
227  if (!Init())
228  return;
229  if (M.end() == M.begin() + dim2_) {
230  real t[dim2_];
231  IntReverse(X, Y, Z, lat, lon, h, t);
232  std::copy(t, t + dim2_, M.begin());
233  } else
234  IntReverse(X, Y, Z, lat, lon, h, NULL);
235  }
236 
237  /** \name Inspector functions
238  **********************************************************************/
239  ///@{
240  /**
241  * @return true if the object has been initialized.
242  **********************************************************************/
243  bool Init() const { return _a > 0; }
244  /**
245  * @return \e a the equatorial radius of the ellipsoid (meters). This is
246  * the value used in the constructor.
247  **********************************************************************/
249  { return Init() ? _a : Math::NaN(); }
250 
251  /**
252  * @return \e f the flattening of the ellipsoid. This is the
253  * value used in the constructor.
254  **********************************************************************/
256  { return Init() ? _f : Math::NaN(); }
257 
258  /**
259  * \deprecated An old name for EquatorialRadius().
260  **********************************************************************/
261  // GEOGRAPHICLIB_DEPRECATED("Use EquatorialRadius()")
262  Math::real MajorRadius() const { return EquatorialRadius(); }
263  ///@}
264 
265  /**
266  * A global instantiation of Geocentric with the parameters for the WGS84
267  * ellipsoid.
268  **********************************************************************/
269  static const Geocentric& WGS84();
270  };
271 
272 } // namespace GeographicLib
273 
274 #endif // GEOGRAPHICLIB_GEOCENTRIC_HPP
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:92
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
The normal gravity of the earth.
Math::real EquatorialRadius() const
Definition: Geocentric.hpp:248
Model of the earth&#39;s magnetic field.
Geomagnetic field on a circle of latitude.
Geocentric coordinates
Definition: Geocentric.hpp:67
void Forward(real lat, real lon, real h, real &X, real &Y, real &Z) const
Definition: Geocentric.hpp:132
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
static T NaN()
Definition: Math.cpp:389
Math::real Flattening() const
Definition: Geocentric.hpp:255
void Reverse(real X, real Y, real Z, real &lat, real &lon, real &h, std::vector< real > &M) const
Definition: Geocentric.hpp:224
Local cartesian coordinates.
Model of the earth&#39;s gravity field.
Header for GeographicLib::Constants class.
Math::real MajorRadius() const
Definition: Geocentric.hpp:262
void Forward(real lat, real lon, real h, real &X, real &Y, real &Z, std::vector< real > &M) const
Definition: Geocentric.hpp:161
void Reverse(real X, real Y, real Z, real &lat, real &lon, real &h) const
Definition: Geocentric.hpp:195
Gravity on a circle of latitude.