00001 /** 00002 * \file CassiniSoldner.hpp 00003 * \brief Header for GeographicLib::CassiniSoldner class 00004 * 00005 * Copyright (c) Charles Karney (2009, 2010) <charles@karney.com> 00006 * and licensed under the LGPL. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #if !defined(GEOGRAPHICLIB_CASSINISOLDNER_HPP) 00011 #define GEOGRAPHICLIB_CASSINISOLDNER_HPP "$Id: CassiniSoldner.hpp 6827 2010-05-20 19:56:18Z karney $" 00012 00013 #include "GeographicLib/Geodesic.hpp" 00014 #include "GeographicLib/Constants.hpp" 00015 00016 namespace GeographicLib { 00017 00018 /** 00019 * \brief Cassini-Soldner Projection. 00020 * 00021 * Cassini-Soldner projection centered at an arbitrary position, \e lat0, \e 00022 * lon0, on the ellipsoid. This projection is a transverse cylindrical 00023 * equidistant projection. The projection from (\e lat, \e lon) to easting 00024 * and northing (\e x, \e y) is defined by geodesics as follows. Go north 00025 * along a geodesic a distance \e y from the central point; then turn 00026 * clockwise 90<sup>o</sup> and go a distance \e x along a geodesic. 00027 * (Although the initial heading is north, this changes to south if the pole 00028 * is crossed.) This procedure uniquely defines the reverse projection. The 00029 * forward projection is constructed as follows. Find the point (\e lat1, \e 00030 * lon1) on the meridian closest to (\e lat, \e lon). Here we consider the 00031 * full meridian so that \e lon1 may be either \e lon0 or \e lon0 + 00032 * 180<sup>o</sup>. \e x is the geodesic distance from (\e lat1, \e lon1) to 00033 * (\e lat, \e lon), appropriately signed according to which side of the 00034 * central meridian (\e lat, \e lon) lies. \e y is the shortest distance 00035 * along the meridian from (\e lat0, \e lon0) to (\e lat1, \e lon1), again, 00036 * appropriately signed according to the initial heading. [Note that, in the 00037 * case of prolate ellipsoids, the shortest meridional path from (\e lat0, \e 00038 * lon0) to (\e lat1, \e lon1) may not be the shortest path.] This procedure 00039 * uniquely defines the forward projection except for a small class of points 00040 * for which there may be two equally short routes for either leg of the 00041 * path. 00042 * 00043 * Because of the properties of geodesics, the (\e x, \e y) grid is 00044 * orthogonal. The scale in the easting direction is unity. The scale, \e 00045 * k, in the northing direction is unity on the central meridian and 00046 * increases away from the central meridian. The projection routines return 00047 * \e azi, the true bearing of the easting direction, and \e rk = 1/\e k, the 00048 * reciprocal of the scale in the northing direction. 00049 * 00050 * The conversions all take place using a GeographicLib::Geodesic object (by 00051 * default GeographicLib::Geodesic::WGS84). For more information on 00052 * geodesics see \ref geodesic. The determination of (\e lat1, \e lon1) in 00053 * the forward projection is by solving the inverse geodesic problem for (\e 00054 * lat, \e lon) and its twin obtained by reflection in the meridional plane. 00055 * The scale is found by determining where two neighboring geodesics 00056 * intersecting the central meridan at \e lat1 and \e lat1 + \e dlat1 00057 * intersect and taking the ratio of the reduced lengths for the two 00058 * geodesics between that point and, respectively, (\e lat1, \e lon1) and (\e 00059 * lat, \e lon). 00060 **********************************************************************/ 00061 00062 class CassiniSoldner { 00063 private: 00064 typedef Math::real real; 00065 const Geodesic _earth; 00066 GeodesicLine _meridian; 00067 real _sbet0, _cbet0; 00068 static const real eps1, eps2; 00069 static const unsigned maxit = 10; 00070 00071 static inline real sq(real x) throw() { return x * x; } 00072 // The following private helper functions are copied from Geodesic. 00073 static inline real AngNormalize(real x) throw() { 00074 // Place angle in [-180, 180). Assumes x is in [-540, 540). 00075 return x >= 180 ? x - 360 : x < -180 ? x + 360 : x; 00076 } 00077 static inline real AngRound(real x) throw() { 00078 // The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57 00079 // for reals = 0.7 pm on the earth if x is an angle in degrees. (This 00080 // is about 1000 times more resolution than we get with angles around 90 00081 // degrees.) We use this to avoid having to deal with near singular 00082 // cases when x is non-zero but tiny (e.g., 1.0e-200). 00083 const real z = real(0.0625); // 1/16 00084 volatile real y = std::abs(x); 00085 // The compiler mustn't "simplify" z - (z - y) to y 00086 y = y < z ? z - (z - y) : y; 00087 return x < 0 ? -y : y; 00088 } 00089 static inline void SinCosNorm(real& sinx, real& cosx) throw() { 00090 real r = Math::hypot(sinx, cosx); 00091 sinx /= r; 00092 cosx /= r; 00093 } 00094 public: 00095 00096 /** 00097 * Constructor for CassiniSoldner setting the Geodesic object, \e earth, to 00098 * use for geodesic calculations. By default this uses the WGS84 00099 * ellipsoid. This constructor makes an "uninitialized" object. Call 00100 * Reset to set the central latitude and longuitude, prior to calling 00101 * Forward and Reverse. 00102 **********************************************************************/ 00103 explicit CassiniSoldner(const Geodesic& earth = Geodesic::WGS84) throw() 00104 : _earth(earth) {} 00105 00106 /** 00107 * Constructor for CassiniSoldner setting the center point, \e lat0, \e 00108 * lon0 (degrees) of the projection and the Geodesic object, \e earth, to 00109 * use for geodesic calculations. By default this uses the WGS84 00110 * ellipsoid. 00111 **********************************************************************/ 00112 CassiniSoldner(real lat0, real lon0, 00113 const Geodesic& earth = Geodesic::WGS84) throw() 00114 : _earth(earth) { 00115 Reset(lat0, lon0); 00116 } 00117 00118 /** 00119 * Set the central latititude to \e lat0 and central longitude to \e lon0 00120 * (degrees). \e lat0 should be in the range [-90, 90] and \e lon0 should 00121 * be in the range [-180, 360]. 00122 **********************************************************************/ 00123 void Reset(real lat0, real lon0) throw(); 00124 00125 /** 00126 * Convert from latitude \e lat (degrees) and longitude \e lon (degrees) to 00127 * Cassini-Soldner easting \e x (meters) and northing \e y (meters). Also 00128 * return the azimuth of the easting direction \e azi (degrees) and the 00129 * reciprocal of the northing scale \e rk. \e lat should be in the range 00130 * [-90, 90] and \e lon should be in the range [-180, 360]. A call to 00131 * Forward followed by a call to Reverse will return the original (\e lat, 00132 * \e lon) (to within roundoff). The routine does nothing if the origin 00133 * has not been set. 00134 **********************************************************************/ 00135 void Forward(real lat, real lon, 00136 real& x, real& y, real& azi, real& rk) const throw(); 00137 00138 /** 00139 * Convert from Cassini-Soldner easting \e x (meters) and northing \e y 00140 * (meters) to latitude \e lat (degrees) and longitude \e lon (degrees). 00141 * Also return the azimuth of the easting direction \e azi (degrees) and 00142 * the reciprocal of the northing scale \e rk. A call to Reverse followed 00143 * by a call to Forward will return the original (\e x, \e y) (to within 00144 * roundoff), provided that \e x and \e y are sufficiently small not to 00145 * "wrap around" the earth. The routine does nothing if the origin has not 00146 * been set. 00147 **********************************************************************/ 00148 void Reverse(real x, real y, 00149 real& lat, real& lon, real& azi, real& rk) const throw(); 00150 00151 /** 00152 * Has this object been initialized with an origin? 00153 **********************************************************************/ 00154 bool Init() const throw() { return _meridian.Init(); } 00155 00156 /** 00157 * Return the latitude of the origin (degrees). 00158 **********************************************************************/ 00159 Math::real LatitudeOrigin() const throw() 00160 { return _meridian.Latitude(); } 00161 00162 /** 00163 * Return the longitude of the origin (degrees). 00164 **********************************************************************/ 00165 Math::real LongitudeOrigin() const throw() 00166 { return _meridian.Longitude(); } 00167 00168 /** 00169 * The major radius of the ellipsoid (meters). This is that value of \e a 00170 * inherited from the Geodesic object used in the constructor. 00171 **********************************************************************/ 00172 Math::real MajorRadius() const throw() { return _earth.MajorRadius(); } 00173 00174 /** 00175 * The inverse flattening of the ellipsoid. This is that value of \e r 00176 * inherited from the Geodesic object used in the constructor. A value of 00177 * 0 is returned for a sphere (infinite inverse flattening). 00178 **********************************************************************/ 00179 Math::real InverseFlattening() const throw() 00180 { return _earth.InverseFlattening(); } 00181 }; 00182 00183 } // namespace GeographicLib 00184 00185 #endif