GeographicLib
1.21
|
00001 /** 00002 * \file CassiniSoldner.hpp 00003 * \brief Header for GeographicLib::CassiniSoldner class 00004 * 00005 * Copyright (c) Charles Karney (2009-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_CASSINISOLDNER_HPP) 00011 #define GEOGRAPHICLIB_CASSINISOLDNER_HPP \ 00012 "$Id: d794ea8a1e64fd9cbb8dcee34755b6dc4fee623a $" 00013 00014 #include <GeographicLib/Geodesic.hpp> 00015 #include <GeographicLib/GeodesicLine.hpp> 00016 #include <GeographicLib/Constants.hpp> 00017 00018 namespace GeographicLib { 00019 00020 /** 00021 * \brief Cassini-Soldner Projection. 00022 * 00023 * Cassini-Soldner projection centered at an arbitrary position, \e lat0, \e 00024 * lon0, on the ellipsoid. This projection is a transverse cylindrical 00025 * equidistant projection. The projection from (\e lat, \e lon) to easting 00026 * and northing (\e x, \e y) is defined by geodesics as follows. Go north 00027 * along a geodesic a distance \e y from the central point; then turn 00028 * clockwise 90<sup>o</sup> and go a distance \e x along a geodesic. 00029 * (Although the initial heading is north, this changes to south if the pole 00030 * is crossed.) This procedure uniquely defines the reverse projection. The 00031 * forward projection is constructed as follows. Find the point (\e lat1, \e 00032 * lon1) on the meridian closest to (\e lat, \e lon). Here we consider the 00033 * full meridian so that \e lon1 may be either \e lon0 or \e lon0 + 00034 * 180<sup>o</sup>. \e x is the geodesic distance from (\e lat1, \e lon1) to 00035 * (\e lat, \e lon), appropriately signed according to which side of the 00036 * central meridian (\e lat, \e lon) lies. \e y is the shortest distance 00037 * along the meridian from (\e lat0, \e lon0) to (\e lat1, \e lon1), again, 00038 * appropriately signed according to the initial heading. [Note that, in the 00039 * case of prolate ellipsoids, the shortest meridional path from (\e lat0, \e 00040 * lon0) to (\e lat1, \e lon1) may not be the shortest path.] This procedure 00041 * uniquely defines the forward projection except for a small class of points 00042 * for which there may be two equally short routes for either leg of the 00043 * path. 00044 * 00045 * Because of the properties of geodesics, the (\e x, \e y) grid is 00046 * orthogonal. The scale in the easting direction is unity. The scale, \e 00047 * k, in the northing direction is unity on the central meridian and 00048 * increases away from the central meridian. The projection routines return 00049 * \e azi, the true bearing of the easting direction, and \e rk = 1/\e k, the 00050 * reciprocal of the scale in the northing direction. 00051 * 00052 * The conversions all take place using a Geodesic object (by default 00053 * Geodesic::WGS84). For more information on geodesics see \ref geodesic. 00054 * The determination of (\e lat1, \e lon1) in the forward projection is by 00055 * solving the inverse geodesic problem for (\e lat, \e lon) and its twin 00056 * obtained by reflection in the meridional plane. The scale is found by 00057 * determining where two neighboring geodesics intersecting the central 00058 * meridian at \e lat1 and \e lat1 + \e dlat1 intersect and taking the ratio 00059 * of the reduced lengths for the two geodesics between that point and, 00060 * respectively, (\e lat1, \e lon1) and (\e lat, \e lon). 00061 * 00062 * Example of use: 00063 * \include example-CassiniSoldner.cpp 00064 * 00065 * <a href="GeodesicProj.1.html">GeodesicProj</a> is a command-line utility 00066 * providing access to the functionality of AzimuthalEquidistant, Gnomonic, 00067 * and CassiniSoldner. 00068 **********************************************************************/ 00069 00070 class GEOGRAPHIC_EXPORT CassiniSoldner { 00071 private: 00072 typedef Math::real real; 00073 Geodesic _earth; 00074 GeodesicLine _meridian; 00075 real _sbet0, _cbet0; 00076 static const real eps1_; 00077 static const real tiny_; 00078 static const unsigned maxit_ = 10; 00079 00080 // The following private helper functions are copied from Geodesic. 00081 static inline real AngNormalize(real x) throw() { 00082 // Place angle in [-180, 180). Assumes x is in [-540, 540). 00083 return x >= 180 ? x - 360 : (x < -180 ? x + 360 : x); 00084 } 00085 static inline real AngRound(real x) throw() { 00086 // The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57 00087 // for reals = 0.7 pm on the earth if x is an angle in degrees. (This 00088 // is about 1000 times more resolution than we get with angles around 90 00089 // degrees.) We use this to avoid having to deal with near singular 00090 // cases when x is non-zero but tiny (e.g., 1.0e-200). 00091 const real z = real(0.0625); // 1/16 00092 volatile real y = std::abs(x); 00093 // The compiler mustn't "simplify" z - (z - y) to y 00094 y = y < z ? z - (z - y) : y; 00095 return x < 0 ? -y : y; 00096 } 00097 static inline void SinCosNorm(real& sinx, real& cosx) throw() { 00098 real r = Math::hypot(sinx, cosx); 00099 sinx /= r; 00100 cosx /= r; 00101 } 00102 public: 00103 00104 /** 00105 * Constructor for CassiniSoldner. 00106 * 00107 * @param[in] earth the Geodesic object to use for geodesic calculations. 00108 * By default this uses the WGS84 ellipsoid. 00109 * 00110 * This constructor makes an "uninitialized" object. Call Reset to set the 00111 * central latitude and longitude, prior to calling Forward and Reverse. 00112 **********************************************************************/ 00113 explicit CassiniSoldner(const Geodesic& earth = Geodesic::WGS84) throw() 00114 : _earth(earth) {} 00115 00116 /** 00117 * Constructor for CassiniSoldner specifying a center point. 00118 * 00119 * @param[in] lat0 latitude of center point of projection (degrees). 00120 * @param[in] lon0 longitude of center point of projection (degrees). 00121 * @param[in] earth the Geodesic object to use for geodesic calculations. 00122 * By default this uses the WGS84 ellipsoid. 00123 * 00124 * \e lat0 should be in the range [-90, 90] and \e lon0 should be in the 00125 * range [-180, 360]. 00126 **********************************************************************/ 00127 CassiniSoldner(real lat0, real lon0, 00128 const Geodesic& earth = Geodesic::WGS84) throw() 00129 : _earth(earth) { 00130 Reset(lat0, lon0); 00131 } 00132 00133 /** 00134 * Set the central point of the projection 00135 * 00136 * @param[in] lat0 latitude of center point of projection (degrees). 00137 * @param[in] lon0 longitude of center point of projection (degrees). 00138 * 00139 * \e lat0 should be in the range [-90, 90] and \e lon0 should be in the 00140 * range [-180, 360]. 00141 **********************************************************************/ 00142 void Reset(real lat0, real lon0) throw(); 00143 00144 /** 00145 * Forward projection, from geographic to Cassini-Soldner. 00146 * 00147 * @param[in] lat latitude of point (degrees). 00148 * @param[in] lon longitude of point (degrees). 00149 * @param[out] x easting of point (meters). 00150 * @param[out] y northing of point (meters). 00151 * @param[out] azi azimuth of easting direction at point (degrees). 00152 * @param[out] rk reciprocal of azimuthal northing scale at point. 00153 * 00154 * \e lat should be in the range [-90, 90] and \e lon should be in the 00155 * range [-180, 360]. A call to Forward followed by a call to Reverse will 00156 * return the original (\e lat, \e lon) (to within roundoff). The routine 00157 * does nothing if the origin has not been set. 00158 **********************************************************************/ 00159 void Forward(real lat, real lon, 00160 real& x, real& y, real& azi, real& rk) const throw(); 00161 00162 /** 00163 * Reverse projection, from Cassini-Soldner to geographic. 00164 * 00165 * @param[in] x easting of point (meters). 00166 * @param[in] y northing of point (meters). 00167 * @param[out] lat latitude of point (degrees). 00168 * @param[out] lon longitude of point (degrees). 00169 * @param[out] azi azimuth of easting direction at point (degrees). 00170 * @param[out] rk reciprocal of azimuthal northing scale at point. 00171 * 00172 * A call to Reverse followed by a call to Forward will return the original 00173 * (\e x, \e y) (to within roundoff), provided that \e x and \e y are 00174 * sufficiently small not to "wrap around" the earth. The routine does 00175 * nothing if the origin has not been set. 00176 **********************************************************************/ 00177 void Reverse(real x, real y, 00178 real& lat, real& lon, real& azi, real& rk) const throw(); 00179 00180 /** 00181 * CassiniSoldner::Forward without returning the azimuth and scale. 00182 **********************************************************************/ 00183 void Forward(real lat, real lon, 00184 real& x, real& y) const throw() { 00185 real azi, rk; 00186 Forward(lat, lon, x, y, azi, rk); 00187 } 00188 00189 /** 00190 * CassiniSoldner::Reverse without returning the azimuth and scale. 00191 **********************************************************************/ 00192 void Reverse(real x, real y, 00193 real& lat, real& lon) const throw() { 00194 real azi, rk; 00195 Reverse(x, y, lat, lon, azi, rk); 00196 } 00197 00198 /** \name Inspector functions 00199 **********************************************************************/ 00200 ///@{ 00201 /** 00202 * @return true if the object has been initialized. 00203 **********************************************************************/ 00204 bool Init() const throw() { return _meridian.Init(); } 00205 00206 /** 00207 * @return \e lat0 the latitude of origin (degrees). 00208 **********************************************************************/ 00209 Math::real LatitudeOrigin() const throw() 00210 { return _meridian.Latitude(); } 00211 00212 /** 00213 * @return \e lon0 the longitude of origin (degrees). 00214 **********************************************************************/ 00215 Math::real LongitudeOrigin() const throw() 00216 { return _meridian.Longitude(); } 00217 00218 /** 00219 * @return \e a the equatorial radius of the ellipsoid (meters). This is 00220 * the value inherited from the Geodesic object used in the constructor. 00221 **********************************************************************/ 00222 Math::real MajorRadius() const throw() { return _earth.MajorRadius(); } 00223 00224 /** 00225 * @return \e f the flattening of the ellipsoid. This is the value 00226 * inherited from the Geodesic object used in the constructor. 00227 **********************************************************************/ 00228 Math::real Flattening() const throw() { return _earth.Flattening(); } 00229 ///@} 00230 00231 /// \cond SKIP 00232 /** 00233 * <b>DEPRECATED</b> 00234 * @return \e r the inverse flattening of the ellipsoid. 00235 **********************************************************************/ 00236 Math::real InverseFlattening() const throw() 00237 { return _earth.InverseFlattening(); } 00238 /// \endcond 00239 }; 00240 00241 } // namespace GeographicLib 00242 00243 #endif // GEOGRAPHICLIB_CASSINISOLDNER_HPP