GeographicLib  1.35
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-2011) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
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="http://dx.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="http://arxiv.org/abs/1102.1215v1">Geodesics
41  * on an ellipsoid of revolution</a>,
42  * Feb. 2011;
43  * preprint
44  * <a href="http://arxiv.org/abs/1102.1215v1">arxiv:1102.1215v1</a>.
45  * .
46  * See \ref geocentric for more information.
47  *
48  * The errors in these routines are close to round-off. Specifically, for
49  * points within 5000 km of the surface of the ellipsoid (either inside or
50  * outside the ellipsoid), the error is bounded by 7 nm (7 nanometers) for
51  * the WGS84 ellipsoid. See \ref geocentric for further information on the
52  * errors.
53  *
54  * Example of use:
55  * \include example-Geocentric.cpp
56  *
57  * <a href="CartConvert.1.html">CartConvert</a> is a command-line utility
58  * providing access to the functionality of Geocentric and LocalCartesian.
59  **********************************************************************/
60 
62  private:
63  typedef Math::real real;
64  friend class LocalCartesian;
65  friend class MagneticCircle; // MagneticCircle uses Rotation
66  friend class MagneticModel; // MagneticModel uses IntForward
67  friend class GravityCircle; // GravityCircle uses Rotation
68  friend class GravityModel; // GravityModel uses IntForward
69  friend class NormalGravity; // NormalGravity uses IntForward
70  friend class SphericalHarmonic;
71  friend class SphericalHarmonic1;
72  friend class SphericalHarmonic2;
73  static const size_t dim_ = 3;
74  static const size_t dim2_ = dim_ * dim_;
75  real _a, _f, _e2, _e2m, _e2a, _e4a, _maxrad;
76  static void Rotation(real sphi, real cphi, real slam, real clam,
77  real M[dim2_]) throw();
78  static void Rotate(real M[dim2_], real x, real y, real z,
79  real& X, real& Y, real& Z) throw() {
80  // Perform [X,Y,Z]^t = M.[x,y,z]^t
81  // (typically local cartesian to geocentric)
82  X = M[0] * x + M[1] * y + M[2] * z;
83  Y = M[3] * x + M[4] * y + M[5] * z;
84  Z = M[6] * x + M[7] * y + M[8] * z;
85  }
86  static void Unrotate(real M[dim2_], real X, real Y, real Z,
87  real& x, real& y, real& z) throw() {
88  // Perform [x,y,z]^t = M^t.[X,Y,Z]^t
89  // (typically geocentric to local cartesian)
90  x = M[0] * X + M[3] * Y + M[6] * Z;
91  y = M[1] * X + M[4] * Y + M[7] * Z;
92  z = M[2] * X + M[5] * Y + M[8] * Z;
93  }
94  void IntForward(real lat, real lon, real h, real& X, real& Y, real& Z,
95  real M[dim2_]) const throw();
96  void IntReverse(real X, real Y, real Z, real& lat, real& lon, real& h,
97  real M[dim2_]) const throw();
98 
99  public:
100 
101  /**
102  * Constructor for a ellipsoid with
103  *
104  * @param[in] a equatorial radius (meters).
105  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
106  * Negative \e f gives a prolate ellipsoid. If \e f > 1, set flattening
107  * to 1/\e f.
108  * @exception GeographicErr if \e a or (1 &minus; \e f ) \e a is not
109  * positive.
110  **********************************************************************/
111  Geocentric(real a, real f);
112 
113  /**
114  * A default constructor (for use by NormalGravity).
115  **********************************************************************/
116  Geocentric() : _a(-1) {}
117 
118  /**
119  * Convert from geodetic to geocentric coordinates.
120  *
121  * @param[in] lat latitude of point (degrees).
122  * @param[in] lon longitude of point (degrees).
123  * @param[in] h height of point above the ellipsoid (meters).
124  * @param[out] X geocentric coordinate (meters).
125  * @param[out] Y geocentric coordinate (meters).
126  * @param[out] Z geocentric coordinate (meters).
127  *
128  * \e lat should be in the range [&minus;90&deg;, 90&deg;]; \e lon
129  * should be in the range [&minus;540&deg;, 540&deg;).
130  **********************************************************************/
131  void Forward(real lat, real lon, real h, real& X, real& Y, real& Z)
132  const throw() {
133  if (Init())
134  IntForward(lat, lon, h, X, Y, Z, NULL);
135  }
136 
137  /**
138  * Convert from geodetic to geocentric coordinates and return rotation
139  * matrix.
140  *
141  * @param[in] lat latitude of point (degrees).
142  * @param[in] lon longitude of point (degrees).
143  * @param[in] h height of point above the ellipsoid (meters).
144  * @param[out] X geocentric coordinate (meters).
145  * @param[out] Y geocentric coordinate (meters).
146  * @param[out] Z geocentric coordinate (meters).
147  * @param[out] M if the length of the vector is 9, fill with the rotation
148  * matrix in row-major order.
149  *
150  * Let \e v be a unit vector located at (\e lat, \e lon, \e h). We can
151  * express \e v as \e column vectors in one of two ways
152  * - in east, north, up coordinates (where the components are relative to a
153  * local coordinate system at (\e lat, \e lon, \e h)); call this
154  * representation \e v1.
155  * - in geocentric \e X, \e Y, \e Z coordinates; call this representation
156  * \e v0.
157  * .
158  * Then we have \e v0 = \e M &sdot; \e v1.
159  **********************************************************************/
160  void Forward(real lat, real lon, real h, real& X, real& Y, real& Z,
161  std::vector<real>& M)
162  const throw() {
163  if (!Init())
164  return;
165  if (M.end() == M.begin() + dim2_) {
166  real t[dim2_];
167  IntForward(lat, lon, h, X, Y, Z, t);
168  copy(t, t + dim2_, M.begin());
169  } else
170  IntForward(lat, lon, h, X, Y, Z, NULL);
171  }
172 
173  /**
174  * Convert from geocentric to geodetic to coordinates.
175  *
176  * @param[in] X geocentric coordinate (meters).
177  * @param[in] Y geocentric coordinate (meters).
178  * @param[in] Z geocentric coordinate (meters).
179  * @param[out] lat latitude of point (degrees).
180  * @param[out] lon longitude of point (degrees).
181  * @param[out] h height of point above the ellipsoid (meters).
182  *
183  * In general there are multiple solutions and the result which maximizes
184  * \e h is returned. If there are still multiple solutions with different
185  * latitudes (applies only if \e Z = 0), then the solution with \e lat > 0
186  * is returned. If there are still multiple solutions with different
187  * longitudes (applies only if \e X = \e Y = 0) then \e lon = 0 is
188  * returned. The value of \e h returned satisfies \e h &ge; &minus; \e a
189  * (1 &minus; <i>e</i><sup>2</sup>) / sqrt(1 &minus; <i>e</i><sup>2</sup>
190  * sin<sup>2</sup>\e lat). The value of \e lon returned is in the range
191  * [&minus;180&deg;, 180&deg;).
192  **********************************************************************/
193  void Reverse(real X, real Y, real Z, real& lat, real& lon, real& h)
194  const throw() {
195  if (Init())
196  IntReverse(X, Y, Z, lat, lon, h, NULL);
197  }
198 
199  /**
200  * Convert from geocentric to geodetic to coordinates.
201  *
202  * @param[in] X geocentric coordinate (meters).
203  * @param[in] Y geocentric coordinate (meters).
204  * @param[in] Z geocentric coordinate (meters).
205  * @param[out] lat latitude of point (degrees).
206  * @param[out] lon longitude of point (degrees).
207  * @param[out] h height of point above the ellipsoid (meters).
208  * @param[out] M if the length of the vector is 9, fill with the rotation
209  * matrix in row-major order.
210  *
211  * Let \e v be a unit vector located at (\e lat, \e lon, \e h). We can
212  * express \e v as \e column vectors in one of two ways
213  * - in east, north, up coordinates (where the components are relative to a
214  * local coordinate system at (\e lat, \e lon, \e h)); call this
215  * representation \e v1.
216  * - in geocentric \e X, \e Y, \e Z coordinates; call this representation
217  * \e v0.
218  * .
219  * Then we have \e v1 = \e M<sup>T</sup> &sdot; \e v0, where \e
220  * M<sup>T</sup> is the transpose of \e M.
221  **********************************************************************/
222  void Reverse(real X, real Y, real Z, real& lat, real& lon, real& h,
223  std::vector<real>& M)
224  const throw() {
225  if (!Init())
226  return;
227  if (M.end() == M.begin() + dim2_) {
228  real t[dim2_];
229  IntReverse(X, Y, Z, lat, lon, h, t);
230  copy(t, t + dim2_, M.begin());
231  } else
232  IntReverse(X, Y, Z, lat, lon, h, NULL);
233  }
234 
235  /** \name Inspector functions
236  **********************************************************************/
237  ///@{
238  /**
239  * @return true if the object has been initialized.
240  **********************************************************************/
241  bool Init() const throw() { return _a > 0; }
242  /**
243  * @return \e a the equatorial radius of the ellipsoid (meters). This is
244  * the value used in the constructor.
245  **********************************************************************/
246  Math::real MajorRadius() const throw()
247  { return Init() ? _a : Math::NaN<real>(); }
248 
249  /**
250  * @return \e f the flattening of the ellipsoid. This is the
251  * value used in the constructor.
252  **********************************************************************/
253  Math::real Flattening() const throw()
254  { return Init() ? _f : Math::NaN<real>(); }
255  ///@}
256 
257  /// \cond SKIP
258  /**
259  * <b>DEPRECATED</b>
260  * @return \e r the inverse flattening of the ellipsoid.
261  **********************************************************************/
262  Math::real InverseFlattening() const throw()
263  { return Init() ? 1/_f : Math::NaN<real>(); }
264  /// \endcond
265 
266  /**
267  * A global instantiation of Geocentric with the parameters for the WGS84
268  * ellipsoid.
269  **********************************************************************/
270  static const Geocentric WGS84;
271  };
272 
273 } // namespace GeographicLib
274 
275 #endif // GEOGRAPHICLIB_GEOCENTRIC_HPP
void Forward(real lat, real lon, real h, real &X, real &Y, real &Z, std::vector< real > &M) const
Definition: Geocentric.hpp:160
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:52
GeographicLib::Math::real real
Definition: GeodSolve.cpp:40
The normal gravity of the earth.
void Forward(real lat, real lon, real h, real &X, real &Y, real &Z) const
Definition: Geocentric.hpp:131
Model of the earth's magnetic field.
Geomagnetic field on a circle of latitude.
Geocentric coordinates
Definition: Geocentric.hpp:61
static const Geocentric WGS84
Definition: Geocentric.hpp:270
Math::real Flattening() const
Definition: Geocentric.hpp:253
Spherical harmonic series with two corrections to the coefficients.
Math::real MajorRadius() const
Definition: Geocentric.hpp:246
Local cartesian coordinates.
Model of the earth's gravity field.
void Reverse(real X, real Y, real Z, real &lat, real &lon, real &h) const
Definition: Geocentric.hpp:193
void Reverse(real X, real Y, real Z, real &lat, real &lon, real &h, std::vector< real > &M) const
Definition: Geocentric.hpp:222
Header for GeographicLib::Constants class.
Spherical harmonic series with a correction to the coefficients.
Spherical harmonic series.
Gravity on a circle of latitude.