CassiniSoldner.hpp

Go to the documentation of this file.
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

Generated on 21 May 2010 for GeographicLib by  doxygen 1.6.1