00001 /** 00002 * \file Gnomonic.hpp 00003 * \brief Header for GeographicLib::Gnomonic class 00004 * 00005 * Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licensed 00006 * under the LGPL. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_GNOMONIC_HPP) 00011 #define GEOGRAPHICLIB_GNOMONIC_HPP "$Id: Gnomonic.hpp 6968 2011-02-19 15:58:56Z karney $" 00012 00013 #include "GeographicLib/Geodesic.hpp" 00014 #include "GeographicLib/GeodesicLine.hpp" 00015 #include "GeographicLib/Constants.hpp" 00016 00017 namespace GeographicLib { 00018 00019 /** 00020 * \brief %Gnomonic Projection. 00021 * 00022 * %Gnomonic projection centered at an arbitrary position \e C on the 00023 * ellipsoid. This projection is derived in Section 13 of 00024 * - C. F. F. Karney, 00025 * <a href="http://arxiv.org/abs/1102.1215">Geodesics 00026 * on an ellipsoid of revolution</a>, 00027 * Feb. 2011; 00028 * preprint 00029 * <a href="http://arxiv.org/abs/1102.1215">arxiv:1102.1215</a>. 00030 * . 00031 * The projection of \e P is defined as follows: compute the 00032 * geodesic line from \e C to \e P; compute the reduced length \e m12, 00033 * geodesic scale \e M12, and \e rho = \e m12/\e M12; finally \e x = \e rho 00034 * sin \e azi1; \e y = \e rho cos \e azi1, where \e azi1 is the azimuth of 00035 * the geodesic at \e C. The Gnomonic::Forward and Gnomonic::Reverse methods 00036 * also return the azimuth \e azi of the geodesic at \e P and reciprocal 00037 * scale \e rk in the azimuthal direction. The scale in the radial direction 00038 * if 1/\e rk<sup>2</sup>. 00039 * 00040 * For a sphere, \e rho is reduces to \e a tan(\e s12/\e a), where \e s12 is 00041 * the length of the geodesic from \e C to \e P, and the gnomonic projection 00042 * has the property that all geodesics appear as straight lines. For an 00043 * ellipsoid, this property holds only for geodesics interesting the centers. 00044 * However geodesic segments close to the center are approximately straight. 00045 * 00046 * Consider a geodesic segment of length \e l. Let \e T be the point on the 00047 * geodesic (extended if necessary) closest to \e C the center of the 00048 * projection and \e t be the distance \e CT. To lowest order, the maximum 00049 * deviation (as a true distance) of the corresponding gnomonic line segment 00050 * (i.e., with the same end points) from the geodesic is<br> 00051 * <br> 00052 * (\e K(T) - \e K(C)) \e l<sup>2</sup> \e t / 32.<br> 00053 * <br> 00054 * where \e K is the Gaussian curvature. 00055 * 00056 * This result applies for any surface. For an allipsoid of revolution, 00057 * consider all geodesics whose end points are within a distance \e r of \e 00058 * C. For a given \e r, the deviation is maximum when the latitude of \e C 00059 * is 45<sup>o</sup>, when endpoints are a distance \e r away, and when their 00060 * azimuths from the center are +/- 45<sup>o</sup> or +/- 135<sup>o</sup>. 00061 * To lowest order in \e r and the flattening \e f, the deviation is \e f 00062 * (\e r/2\e a)<sup>3</sup> \e r. 00063 * 00064 * The conversions all take place using a Geodesic object (by default 00065 * Geodesic::WGS84). For more information on geodesics see \ref geodesic. 00066 * 00067 * <b>CAUTION:</b> The definition of this projection for a sphere is 00068 * standard. However, there is no standard for how it should be extended to 00069 * an ellipsoid. The choices are: 00070 * - Declare that the projection is undefined for an ellipsoid. 00071 * - Project to a tangent plane from the center of the ellipsoid. This 00072 * causes great ellipses to appear as straight lines in the projection; 00073 * i.e., it generalizes the spherical great circle to a great ellipse. 00074 * This was proposed by independently by Bowring and Williams in 1997. 00075 * - Project to the conformal sphere with the constant of integration chosen 00076 * so that the values of the latitude match for the center point and 00077 * perform a central projection onto the plane tangent to the conformal 00078 * sphere at the center point. This causes normal sections through the 00079 * center point to appear as straight lines in the projection; i.e., it 00080 * generalizes the spherical great circle to a normal section. This was 00081 * proposed by I. G. Letoval'tsev, Generalization of the %Gnomonic 00082 * Projection for a Spheroid and the Principal Geodetic Problems Involved 00083 * in the Alignment of Surface Routes, Geodesy and Aerophotography (5), 00084 * 271-274 (1963). 00085 * - The projection given here. This causes geodesics close to the center 00086 * point to appear as straight lines in the projection; i.e., it 00087 * generalizes the spherical great circle to a geodesic. 00088 **********************************************************************/ 00089 00090 class Gnomonic { 00091 private: 00092 typedef Math::real real; 00093 const Geodesic _earth; 00094 real _a, _f; 00095 static const real eps0, eps; 00096 static const int numit = 5; 00097 public: 00098 00099 /** 00100 * Constructor for Gnomonic. 00101 * 00102 * @param[in] earth the Geodesic object to use for geodesic calculations. 00103 * By default this uses the WGS84 ellipsoid. 00104 **********************************************************************/ 00105 explicit Gnomonic(const Geodesic& earth = Geodesic::WGS84) 00106 throw() 00107 : _earth(earth) 00108 , _a(_earth.MajorRadius()) 00109 , _f(_earth.InverseFlattening() ? 00110 1/std::abs(_earth.InverseFlattening()) : 0) 00111 {} 00112 00113 /** 00114 * Forward projection, from geographic to gnomonic. 00115 * 00116 * @param[in] lat0 latitude of center point of projection (degrees). 00117 * @param[in] lon0 longitude of center point of projection (degrees). 00118 * @param[in] lat latitude of point (degrees). 00119 * @param[in] lon longitude of point (degrees). 00120 * @param[out] x easting of point (meters). 00121 * @param[out] y northing of point (meters). 00122 * @param[out] azi azimuth of geodesic at point (degrees). 00123 * @param[out] rk reciprocal of azimuthal scale at point. 00124 * 00125 * \e lat0 and \e lat should be in the range [-90, 90] and \e lon0 and \e 00126 * lon should be in the range [-180, 360]. The scale of the projection is 00127 * 1/\e rk<sup>2</sup> in the "radial" direction, \e azi clockwise from 00128 * true north, and is 1/\e rk in the direction perpendicular to this. If 00129 * the point lies "over the horizon", i.e., if \e rk <= 0, then NaNs are 00130 * returned for \e x and \e y (the correct values are returned for \e azi 00131 * and \e rk). A call to Forward followed by a call to Reverse will return 00132 * the original (\e lat, \e lon) (to within roundoff) provided the point in 00133 * not over the horizon. 00134 **********************************************************************/ 00135 void Forward(real lat0, real lon0, real lat, real lon, 00136 real& x, real& y, real& azi, real& rk) const throw(); 00137 00138 /** 00139 * Reverse projection, from gnomonic to geographic. 00140 * 00141 * @param[in] lat0 latitude of center point of projection (degrees). 00142 * @param[in] lon0 longitude of center point of projection (degrees). 00143 * @param[in] x easting of point (meters). 00144 * @param[in] y northing of point (meters). 00145 * @param[out] lat latitude of point (degrees). 00146 * @param[out] lon longitude of point (degrees). 00147 * @param[out] azi azimuth of geodesic at point (degrees). 00148 * @param[out] rk reciprocal of azimuthal scale at point. 00149 * 00150 * \e lat0 should be in the range [-90, 90] and \e lon0 should be in the 00151 * range [-180, 360]. \e lat will be in the range [-90, 90] and \e lon 00152 * will be in the range [-180, 180). The scale of the projection is 1/\e 00153 * rk<sup>2</sup> in the "radial" direction, \e azi clockwise from true 00154 * north, and is 1/\e rk in the direction perpendicular to this. Even 00155 * though all inputs should return a valid \e lat and \e lon, it's possible 00156 * that the procedure fails to converge for very large \e x or \e y; in 00157 * this case NaNs are returned for all the output arguments. A call to 00158 * Reverse followed by a call to Forward will return the original (\e x, \e 00159 * y) (to roundoff). 00160 **********************************************************************/ 00161 void Reverse(real lat0, real lon0, real x, real y, 00162 real& lat, real& lon, real& azi, real& rk) const throw(); 00163 00164 /** 00165 * Gnomonic::Forward without returning the azimuth and scale. 00166 **********************************************************************/ 00167 void Forward(real lat0, real lon0, real lat, real lon, 00168 real& x, real& y) const throw() { 00169 real azi, rk; 00170 Forward(lat0, lon0, lat, lon, x, y, azi, rk); 00171 } 00172 00173 /** 00174 * Gnomonic::Reverse without returning the azimuth and scale. 00175 **********************************************************************/ 00176 void Reverse(real lat0, real lon0, real x, real y, 00177 real& lat, real& lon) const throw() { 00178 real azi, rk; 00179 Reverse(lat0, lon0, x, y, lat, lon, azi, rk); 00180 } 00181 00182 /** \name Inspector functions 00183 **********************************************************************/ 00184 ///@{ 00185 /** 00186 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00187 * the value inherited from the Geodesic object used in the constructor. 00188 **********************************************************************/ 00189 Math::real MajorRadius() const throw() { return _earth.MajorRadius(); } 00190 00191 /** 00192 * @return \e r the inverse flattening of the ellipsoid. This is the 00193 * value inherited from the Geodesic object used in the constructor. A 00194 * value of 0 is returned for a sphere (infinite inverse flattening). 00195 **********************************************************************/ 00196 Math::real InverseFlattening() const throw() 00197 { return _earth.InverseFlattening(); } 00198 ///@} 00199 }; 00200 00201 } // namespace GeographicLib 00202 00203 #endif