GeographicLib  1.21
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-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