GeographicLib  1.21
Geocentric.hpp
Go to the documentation of this file.
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&ndash;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