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, 2009, 2010, 2011) <charles@karney.com>
00006  * and licensed under the LGPL.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  **********************************************************************/
00009 
00010 #if !defined(GEOGRAPHICLIB_GEOCENTRIC_HPP)
00011 #define GEOGRAPHICLIB_GEOCENTRIC_HPP "$Id: Geocentric.hpp 6967 2011-02-19 15:53:41Z karney $"
00012 
00013 #include "GeographicLib/Constants.hpp"
00014 #include <vector>
00015 #include <algorithm>
00016 
00017 namespace GeographicLib {
00018 
00019   /**
00020    * \brief %Geocentric coordinates
00021    *
00022    * Convert between geodetic coordinates latitude = \e lat, longitude = \e
00023    * lon, height = \e h (measured vertically from the surface of the ellipsoid)
00024    * to geocentric coordinates (\e x, \e y, \e z).  The origin of geocentric
00025    * coordinates is at the center of the earth.  The \e z axis goes thru the
00026    * north pole, \e lat = 90<sup>o</sup>.  The \e x axis goes thru \e lat = 0,
00027    * \e lon = 0.  %Geocentric coordinates are also known as earth centered,
00028    * earth fixed (ECEF) coordinates.
00029    *
00030    * The conversion from geographic to geocentric coordinates is
00031    * straightforward.  For the reverse transformation we use
00032    * - H. Vermeille,
00033    *   <a href="http://dx.doi.org/10.1007/s00190-002-0273-6"> Direct
00034    *   transformation from geocentric coordinates to geodetic coordinates</a>,
00035    *   J. Geodesy 76, 451&ndash;454 (2002).
00036    * .
00037    * Several changes have been made to ensure that the method returns accurate
00038    * results for all finite inputs (even if \e h is infinite).  The changes are
00039    * described in Appendix B of
00040    * - C. F. F. Karney,
00041    *   <a href="http://arxiv.org/abs/1102.1215">Geodesics
00042    *   on an ellipsoid of revolution</a>,
00043    *   Feb. 2011;
00044    *   preprint
00045    *   <a href="http://arxiv.org/abs/1102.1215">arxiv:1102.1215</a>.
00046    * .
00047    * See \ref geocentric for more information.
00048    *
00049    * The errors in these routines are close to round-off.  Specifically, for
00050    * points within 5000 km of the surface of the ellipsoid (either inside or
00051    * outside the ellipsoid), the error is bounded by 7 nm for the WGS84
00052    * ellipsoid.  See \ref geocentric for further information on the errors.
00053    **********************************************************************/
00054 
00055   class Geocentric {
00056   private:
00057     typedef Math::real real;
00058     friend class LocalCartesian;
00059     static const size_t dim = 3, dim2 = dim * dim;
00060     const real _a, _r, _f, _e2, _e2m, _e2a, _e4a, _maxrad;
00061     static inline real sq(real x) throw() { return x * x; }
00062     // Actually this can be static because it doesn't depend on the ellipsoid.
00063     // But let's be more general than that.
00064     void Rotation(real sphi, real cphi, real slam, real clam,
00065                   real M[dim2]) const throw();
00066     void IntForward(real lat, real lon, real h, real& x, real& y, real& z,
00067                     real M[dim2]) const throw();
00068     void IntReverse(real x, real y, real z, real& lat, real& lon, real& h,
00069                     real M[dim2]) const throw();
00070   public:
00071 
00072     /**
00073      * Constructor for a ellipsoid with
00074      *
00075      * @param[in] a equatorial radius (meters)
00076      * @param[in] r reciprocal flattening.  Setting \e r = 0 implies \e r = inf
00077      *   or flattening = 0 (i.e., a sphere).  Negative \e r indicates a prolate
00078      *   ellipsoid.
00079      *
00080      * An exception is thrown if either of the axes of the ellipsoid is
00081      * non-positive.
00082      **********************************************************************/
00083     Geocentric(real a, real r);
00084 
00085     /**
00086      * Convert from geodetic to geocentric coordinates.
00087      *
00088      * @param[in] lat latitude of point (degrees).
00089      * @param[in] lon longitude of point (degrees).
00090      * @param[in] h height of point above the ellipsoid (meters).
00091      * @param[out] x geocentric coordinate (meters).
00092      * @param[out] y geocentric coordinate (meters).
00093      * @param[out] z geocentric coordinate (meters).
00094      *
00095      * \e lat should be in the range [-90, 90]; \e lon and \e lon0 should be in
00096      * the range [-180, 360].
00097      **********************************************************************/
00098     void Forward(real lat, real lon, real h, real& x, real& y, real& z)
00099       const throw() {
00100       IntForward(lat, lon, h, x, y, z, NULL);
00101     }
00102 
00103     /**
00104      * Convert from geodetic to geocentric coordinates and return rotation
00105      * matrix.
00106      *
00107      * @param[in] lat latitude of point (degrees).
00108      * @param[in] lon longitude of point (degrees).
00109      * @param[in] h height of point above the ellipsoid (meters).
00110      * @param[out] x geocentric coordinate (meters).
00111      * @param[out] y geocentric coordinate (meters).
00112      * @param[out] z geocentric coordinate (meters).
00113      * @param[out] M if the length of the vector is 9, fill with the rotation
00114      *   matrix in row-major order.
00115      *
00116      * Pre-multiplying a unit vector in local cartesian coordinates (east,
00117      * north, up) by \e M transforms the vector to geocentric coordinates.
00118      **********************************************************************/
00119     void Forward(real lat, real lon, real h, real& x, real& y, real& z,
00120                  std::vector<real>& M)
00121       const throw() {
00122       real t[dim2];
00123       IntForward(lat, lon, h, x, y, z, t);
00124       if (M.end() == M.begin() + dim2)
00125         copy(t, t + dim2, M.begin());
00126     }
00127 
00128     /**
00129      * Convert from geocentric to geodetic to coordinates.
00130      *
00131      * @param[in] x geocentric coordinate (meters).
00132      * @param[in] y geocentric coordinate (meters).
00133      * @param[in] z geocentric coordinate (meters).
00134      * @param[out] lat latitude of point (degrees).
00135      * @param[out] lon longitude of point (degrees).
00136      * @param[out] h height of point above the ellipsoid (meters).
00137      *
00138      * In general there are multiple solutions and the result which maximizes
00139      * \e h is returned.  If there are still multiple solutions with different
00140      * latitutes (applies only if \e z = 0), then the solution with \e lat > 0
00141      * is returned.  If there are still multiple solutions with different
00142      * longitudes (applies only if \e x = \e y = 0) then \e lon = 0 is
00143      * returned.  The value of \e h returned satisfies \e h >= - \e a (1 - \e
00144      * e<sup>2</sup>) / sqrt(1 - \e e<sup>2</sup> sin<sup>2</sup>\e lat).  The
00145      * value of \e lon returned is in the range [-180, 180).
00146      **********************************************************************/
00147     void Reverse(real x, real y, real z, real& lat, real& lon, real& h)
00148       const throw() {
00149       IntReverse(x, y, z, lat, lon, h, NULL);
00150     }
00151 
00152     /**
00153      * Convert from geocentric to geodetic to coordinates.
00154      *
00155      * @param[in] x geocentric coordinate (meters).
00156      * @param[in] y geocentric coordinate (meters).
00157      * @param[in] z geocentric coordinate (meters).
00158      * @param[out] lat latitude of point (degrees).
00159      * @param[out] lon longitude of point (degrees).
00160      * @param[out] h height of point above the ellipsoid (meters).
00161      * @param[out] M if the length of the vector is 9, fill with the rotation
00162      *   matrix in row-major order.
00163      *
00164      * Pre-multiplying a unit vector in geocentric coordinates by the transpose
00165      * of \e M transforms the vector to local cartesian coordinates (east,
00166      * north, up).
00167      **********************************************************************/
00168     void Reverse(real x, real y, real z, real& lat, real& lon, real& h,
00169                  std::vector<real>& M)
00170       const throw() {
00171       real t[dim2];
00172       IntReverse(x, y, z, lat, lon, h, t);
00173       if (M.end() == M.begin() + dim2)
00174         copy(t, t + dim2, M.begin());
00175     }
00176 
00177     /** \name Inspector functions
00178      **********************************************************************/
00179     ///@{
00180     /**
00181      * @return \e a the equatorial radius of the ellipsoid (meters).  This is
00182      *   the value used in the constructor.
00183      **********************************************************************/
00184     Math::real MajorRadius() const throw() { return _a; }
00185 
00186     /**
00187      * @return \e r the inverse flattening of the ellipsoid.  This is the
00188      *   value used in the constructor.  A value of 0 is returned for a sphere
00189      *   (infinite inverse flattening).
00190      **********************************************************************/
00191     Math::real InverseFlattening() const throw() { return _r; }
00192     ///@}
00193 
00194     /**
00195      * A global instantiation of Geocentric with the parameters for the WGS84
00196      * ellipsoid.
00197      **********************************************************************/
00198     static const Geocentric WGS84;
00199   };
00200 
00201 } // namespace GeographicLib
00202 #endif