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 6867 2010-09-11 13:04:26Z 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 Cassini-Soldner Projection.
00021    *
00022    * Cassini-Soldner projection centered at an arbitrary position, \e lat0, \e
00023    * lon0, on the ellipsoid.  This projection is a transverse cylindrical
00024    * equidistant projection.  The projection from (\e lat, \e lon) to easting
00025    * and northing (\e x, \e y) is defined by geodesics as follows.  Go north
00026    * along a geodesic a distance \e y from the central point; then turn
00027    * clockwise 90<sup>o</sup> and go a distance \e x along a geodesic.
00028    * (Although the initial heading is north, this changes to south if the pole
00029    * is crossed.)  This procedure uniquely defines the reverse projection.  The
00030    * forward projection is constructed as follows.  Find the point (\e lat1, \e
00031    * lon1) on the meridian closest to (\e lat, \e lon).  Here we consider the
00032    * full meridian so that \e lon1 may be either \e lon0 or \e lon0 +
00033    * 180<sup>o</sup>.  \e x is the geodesic distance from (\e lat1, \e lon1) to
00034    * (\e lat, \e lon), appropriately signed according to which side of the
00035    * central meridian (\e lat, \e lon) lies.  \e y is the shortest distance
00036    * along the meridian from (\e lat0, \e lon0) to (\e lat1, \e lon1), again,
00037    * appropriately signed according to the initial heading.  [Note that, in the
00038    * case of prolate ellipsoids, the shortest meridional path from (\e lat0, \e
00039    * lon0) to (\e lat1, \e lon1) may not be the shortest path.]  This procedure
00040    * uniquely defines the forward projection except for a small class of points
00041    * for which there may be two equally short routes for either leg of the
00042    * path.
00043    *
00044    * Because of the properties of geodesics, the (\e x, \e y) grid is
00045    * orthogonal.  The scale in the easting direction is unity.  The scale, \e
00046    * k, in the northing direction is unity on the central meridian and
00047    * increases away from the central meridian.  The projection routines return
00048    * \e azi, the true bearing of the easting direction, and \e rk = 1/\e k, the
00049    * reciprocal of the scale in the northing direction.
00050    *
00051    * The conversions all take place using a Geodesic object (by default
00052    * Geodesic::WGS84).  For more information on geodesics see \ref geodesic.
00053    * The determination of (\e lat1, \e lon1) in the forward projection is by
00054    * solving the inverse geodesic problem for (\e lat, \e lon) and its twin
00055    * obtained by reflection in the meridional plane.  The scale is found by
00056    * determining where two neighboring geodesics intersecting the central
00057    * meridan at \e lat1 and \e lat1 + \e dlat1 intersect and taking the ratio
00058    * of the reduced lengths for the two geodesics between that point and,
00059    * respectively, (\e lat1, \e lon1) and (\e 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.
00098      *
00099      * @param[in] earth the Geodesic object to use for geodesic calculations.
00100      *   By default this uses the WGS84 ellipsoid.
00101      *
00102      * This constructor makes an "uninitialized" object.  Call Reset to set the
00103      * central latitude and longuitude, prior to calling Forward and Reverse.
00104      **********************************************************************/
00105     explicit CassiniSoldner(const Geodesic& earth = Geodesic::WGS84) throw()
00106       : _earth(earth) {}
00107 
00108     /**
00109      * Constructor for CassiniSoldner specifying a center point.
00110      *
00111      * @param[in] lat0 latitude of center point of projection (degrees).
00112      * @param[in] lon0 longitude of center point of projection (degrees).
00113      * @param[in] earth the Geodesic object to use for geodesic calculations.
00114      *   By default this uses the WGS84 ellipsoid.
00115      *
00116      * \e lat0 should be in the range [-90, 90] and \e lon0 should be in the
00117      * range [-180, 360].
00118      **********************************************************************/
00119     CassiniSoldner(real lat0, real lon0,
00120                    const Geodesic& earth = Geodesic::WGS84) throw()
00121       : _earth(earth) {
00122       Reset(lat0, lon0);
00123     }
00124 
00125     /**
00126      * Set the central point of the projection
00127      *
00128      * @param[in] lat0 latitude of center point of projection (degrees).
00129      * @param[in] lon0 longitude of center point of projection (degrees).
00130      *
00131      * \e lat0 should be in the range [-90, 90] and \e lon0 should be in the
00132      * range [-180, 360].
00133      **********************************************************************/
00134     void Reset(real lat0, real lon0) throw();
00135 
00136     /**
00137      * Forward projection, from geographic to Cassini-Soldner.
00138      *
00139      * @param[in] lat latitude of point (degrees).
00140      * @param[in] lon longitude of point (degrees).
00141      * @param[out] x easting of point (meters).
00142      * @param[out] y northing of point (meters).
00143      * @param[out] azi azimuth of easting direction at point (degrees).
00144      * @param[out] rk reciprocal of azimuthal northing scale at point.
00145      *
00146      * \e lat should be in the range [-90, 90] and \e lon should be in the
00147      * range [-180, 360].  A call to Forward followed by a call to Reverse will
00148      * return the original (\e lat, \e lon) (to within roundoff).  The routine
00149      * does nothing if the origin has not been set.
00150      **********************************************************************/
00151     void Forward(real lat, real lon,
00152                  real& x, real& y, real& azi, real& rk) const throw();
00153 
00154     /**
00155      * Reverse projection, from Cassini-Soldner to geographic.
00156      *
00157      * @param[in] x easting of point (meters).
00158      * @param[in] y northing of point (meters).
00159      * @param[out] lat latitude of point (degrees).
00160      * @param[out] lon longitude of point (degrees).
00161      * @param[out] azi azimuth of easting direction at point (degrees).
00162      * @param[out] rk reciprocal of azimuthal northing scale at point.
00163      *
00164      * A call to Reverse followed by a call to Forward will return the original
00165      * (\e x, \e y) (to within roundoff), provided that \e x and \e y are
00166      * sufficiently small not to "wrap around" the earth.  The routine does
00167      * nothing if the origin has not been set.
00168      **********************************************************************/
00169     void Reverse(real x, real y,
00170                  real& lat, real& lon, real& azi, real& rk) const throw();
00171 
00172     /**
00173      * CassiniSoldner::Forward without returning the azimuth and scale.
00174      **********************************************************************/
00175     void Forward(real lat, real lon,
00176                  real& x, real& y) const throw() {
00177       real azi, rk;
00178       Forward(lat, lon, x, y, azi, rk);
00179     }
00180 
00181     /**
00182      * CassiniSoldner::Reverse without returning the azimuth and scale.
00183      **********************************************************************/
00184     void Reverse(real x, real y,
00185                  real& lat, real& lon) const throw() {
00186       real azi, rk;
00187       Reverse(x, y, lat, lon, azi, rk);
00188     }
00189 
00190     /** \name Inspector functions
00191      **********************************************************************/
00192     ///@{
00193     /**
00194      * @return true if the object has been initialized.
00195      **********************************************************************/
00196     bool Init() const throw() { return _meridian.Init(); }
00197 
00198     /**
00199      * @return \e lat0 the latitude of origin (degrees).
00200      **********************************************************************/
00201     Math::real LatitudeOrigin() const throw()
00202     { return _meridian.Latitude(); }
00203 
00204     /**
00205      * @return \e lon0 the longitude of origin (degrees).
00206      **********************************************************************/
00207     Math::real LongitudeOrigin() const throw()
00208     { return _meridian.Longitude(); }
00209 
00210     /**
00211      * @return \e a the equatorial radius of the ellipsoid (meters).  This is
00212      *   the value inherited from the Geodesic object used in the constructor.
00213      **********************************************************************/
00214     Math::real MajorRadius() const throw() { return _earth.MajorRadius(); }
00215 
00216     /**
00217      * @return \e r the inverse flattening of the ellipsoid.  This is the
00218      *   value inherited from the Geodesic object used in the constructor.  A
00219      *   value of 0 is returned for a sphere (infinite inverse flattening).
00220      **********************************************************************/
00221     Math::real InverseFlattening() const throw()
00222     { return _earth.InverseFlattening(); }
00223     ///@}
00224   };
00225 
00226 } // namespace GeographicLib
00227 
00228 #endif