GeographicLib  1.40
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
CassiniSoldner.hpp
Go to the documentation of this file.
1 /**
2  * \file CassiniSoldner.hpp
3  * \brief Header for GeographicLib::CassiniSoldner class
4  *
5  * Copyright (c) Charles Karney (2009-2011) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
8  **********************************************************************/
9 
10 #if !defined(GEOGRAPHICLIB_CASSINISOLDNER_HPP)
11 #define GEOGRAPHICLIB_CASSINISOLDNER_HPP 1
12 
16 
17 namespace GeographicLib {
18 
19  /**
20  * \brief Cassini-Soldner projection
21  *
22  * Cassini-Soldner projection centered at an arbitrary position, \e lat0, \e
23  * lon0, on the ellipsoid. This projection is a transverse cylindrical
24  * equidistant projection. The projection from (\e lat, \e lon) to easting
25  * and northing (\e x, \e y) is defined by geodesics as follows. Go north
26  * along a geodesic a distance \e y from the central point; then turn
27  * clockwise 90&deg; and go a distance \e x along a geodesic.
28  * (Although the initial heading is north, this changes to south if the pole
29  * is crossed.) This procedure uniquely defines the reverse projection. The
30  * forward projection is constructed as follows. Find the point (\e lat1, \e
31  * lon1) on the meridian closest to (\e lat, \e lon). Here we consider the
32  * full meridian so that \e lon1 may be either \e lon0 or \e lon0 +
33  * 180&deg;. \e x is the geodesic distance from (\e lat1, \e lon1) to
34  * (\e lat, \e lon), appropriately signed according to which side of the
35  * central meridian (\e lat, \e lon) lies. \e y is the shortest distance
36  * along the meridian from (\e lat0, \e lon0) to (\e lat1, \e lon1), again,
37  * appropriately signed according to the initial heading. [Note that, in the
38  * case of prolate ellipsoids, the shortest meridional path from (\e lat0, \e
39  * lon0) to (\e lat1, \e lon1) may not be the shortest path.] This procedure
40  * uniquely defines the forward projection except for a small class of points
41  * for which there may be two equally short routes for either leg of the
42  * path.
43  *
44  * Because of the properties of geodesics, the (\e x, \e y) grid is
45  * orthogonal. The scale in the easting direction is unity. The scale, \e
46  * k, in the northing direction is unity on the central meridian and
47  * increases away from the central meridian. The projection routines return
48  * \e azi, the true bearing of the easting direction, and \e rk = 1/\e k, the
49  * reciprocal of the scale in the northing direction.
50  *
51  * The conversions all take place using a Geodesic object (by default
52  * Geodesic::WGS84()). For more information on geodesics see \ref geodesic.
53  * The determination of (\e lat1, \e lon1) in the forward projection is by
54  * solving the inverse geodesic problem for (\e lat, \e lon) and its twin
55  * obtained by reflection in the meridional plane. The scale is found by
56  * determining where two neighboring geodesics intersecting the central
57  * meridian at \e lat1 and \e lat1 + \e dlat1 intersect and taking the ratio
58  * of the reduced lengths for the two geodesics between that point and,
59  * respectively, (\e lat1, \e lon1) and (\e lat, \e lon).
60  *
61  * Example of use:
62  * \include example-CassiniSoldner.cpp
63  *
64  * <a href="GeodesicProj.1.html">GeodesicProj</a> is a command-line utility
65  * providing access to the functionality of AzimuthalEquidistant, Gnomonic,
66  * and CassiniSoldner.
67  **********************************************************************/
68 
70  private:
71  typedef Math::real real;
72  real eps1_, tiny_;
73  Geodesic _earth;
74  GeodesicLine _meridian;
75  real _sbet0, _cbet0;
76  static const unsigned maxit_ = 10;
77 
78  // The following private helper functions are copied from Geodesic.
79  static inline real AngRound(real x) {
80  // The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57
81  // for reals = 0.7 pm on the earth if x is an angle in degrees. (This
82  // is about 1000 times more resolution than we get with angles around 90
83  // degrees.) We use this to avoid having to deal with near singular
84  // cases when x is non-zero but tiny (e.g., 1.0e-200).
85  using std::abs;
86  const real z = 1/real(16);
87  GEOGRAPHICLIB_VOLATILE real y = abs(x);
88  // The compiler mustn't "simplify" z - (z - y) to y
89  y = y < z ? z - (z - y) : y;
90  return x < 0 ? -y : y;
91  }
92  static inline void SinCosNorm(real& sinx, real& cosx) {
93  real r = Math::hypot(sinx, cosx);
94  sinx /= r;
95  cosx /= r;
96  }
97  public:
98 
99  /**
100  * Constructor for CassiniSoldner.
101  *
102  * @param[in] earth the Geodesic object to use for geodesic calculations.
103  * By default this uses the WGS84 ellipsoid.
104  *
105  * This constructor makes an "uninitialized" object. Call Reset to set the
106  * central latitude and longitude, prior to calling Forward and Reverse.
107  **********************************************************************/
108  explicit CassiniSoldner(const Geodesic& earth = Geodesic::WGS84());
109 
110  /**
111  * Constructor for CassiniSoldner specifying a center point.
112  *
113  * @param[in] lat0 latitude of center point of projection (degrees).
114  * @param[in] lon0 longitude of center point of projection (degrees).
115  * @param[in] earth the Geodesic object to use for geodesic calculations.
116  * By default this uses the WGS84 ellipsoid.
117  *
118  * \e lat0 should be in the range [&minus;90&deg;, 90&deg;] and \e
119  * lon0 should be in the range [&minus;540&deg;, 540&deg;).
120  **********************************************************************/
121  CassiniSoldner(real lat0, real lon0,
122  const Geodesic& earth = Geodesic::WGS84());
123 
124  /**
125  * Set the central point of the projection
126  *
127  * @param[in] lat0 latitude of center point of projection (degrees).
128  * @param[in] lon0 longitude of center point of projection (degrees).
129  *
130  * \e lat0 should be in the range [&minus;90&deg;, 90&deg;] and \e
131  * lon0 should be in the range [&minus;540&deg;, 540&deg;).
132  **********************************************************************/
133  void Reset(real lat0, real lon0);
134 
135  /**
136  * Forward projection, from geographic to Cassini-Soldner.
137  *
138  * @param[in] lat latitude of point (degrees).
139  * @param[in] lon longitude of point (degrees).
140  * @param[out] x easting of point (meters).
141  * @param[out] y northing of point (meters).
142  * @param[out] azi azimuth of easting direction at point (degrees).
143  * @param[out] rk reciprocal of azimuthal northing scale at point.
144  *
145  * \e lat should be in the range [&minus;90&deg;, 90&deg;] and \e
146  * lon should be in the range [&minus;540&deg;, 540&deg;). A call
147  * to Forward followed by a call to Reverse will return the original (\e
148  * lat, \e lon) (to within roundoff). The routine does nothing if the
149  * origin has not been set.
150  **********************************************************************/
151  void Forward(real lat, real lon,
152  real& x, real& y, real& azi, real& rk) const;
153 
154  /**
155  * Reverse projection, from Cassini-Soldner to geographic.
156  *
157  * @param[in] x easting of point (meters).
158  * @param[in] y northing of point (meters).
159  * @param[out] lat latitude of point (degrees).
160  * @param[out] lon longitude of point (degrees).
161  * @param[out] azi azimuth of easting direction at point (degrees).
162  * @param[out] rk reciprocal of azimuthal northing scale at point.
163  *
164  * A call to Reverse followed by a call to Forward will return the original
165  * (\e x, \e y) (to within roundoff), provided that \e x and \e y are
166  * sufficiently small not to "wrap around" the earth. The routine does
167  * nothing if the origin has not been set.
168  **********************************************************************/
169  void Reverse(real x, real y,
170  real& lat, real& lon, real& azi, real& rk) const;
171 
172  /**
173  * CassiniSoldner::Forward without returning the azimuth and scale.
174  **********************************************************************/
175  void Forward(real lat, real lon,
176  real& x, real& y) const {
177  real azi, rk;
178  Forward(lat, lon, x, y, azi, rk);
179  }
180 
181  /**
182  * CassiniSoldner::Reverse without returning the azimuth and scale.
183  **********************************************************************/
184  void Reverse(real x, real y,
185  real& lat, real& lon) const {
186  real azi, rk;
187  Reverse(x, y, lat, lon, azi, rk);
188  }
189 
190  /** \name Inspector functions
191  **********************************************************************/
192  ///@{
193  /**
194  * @return true if the object has been initialized.
195  **********************************************************************/
196  bool Init() const { return _meridian.Init(); }
197 
198  /**
199  * @return \e lat0 the latitude of origin (degrees).
200  **********************************************************************/
202  { return _meridian.Latitude(); }
203 
204  /**
205  * @return \e lon0 the longitude of origin (degrees).
206  **********************************************************************/
208  { return _meridian.Longitude(); }
209 
210  /**
211  * @return \e a the equatorial radius of the ellipsoid (meters). This is
212  * the value inherited from the Geodesic object used in the constructor.
213  **********************************************************************/
214  Math::real MajorRadius() const { return _earth.MajorRadius(); }
215 
216  /**
217  * @return \e f the flattening of the ellipsoid. This is the value
218  * inherited from the Geodesic object used in the constructor.
219  **********************************************************************/
220  Math::real Flattening() const { return _earth.Flattening(); }
221  ///@}
222 
223  /// \cond SKIP
224  /**
225  * <b>DEPRECATED</b>
226  * @return \e r the inverse flattening of the ellipsoid.
227  **********************************************************************/
228  Math::real InverseFlattening() const
229  { return _earth.InverseFlattening(); }
230  /// \endcond
231  };
232 
233 } // namespace GeographicLib
234 
235 #endif // GEOGRAPHICLIB_CASSINISOLDNER_HPP
Header for GeographicLib::GeodesicLine class.
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:69
Math::real LatitudeOrigin() const
void Reverse(real x, real y, real &lat, real &lon) const
GeographicLib::Math::real real
Definition: GeodSolve.cpp:32
static const Geodesic & WGS84()
Definition: Geodesic.cpp:90
Cassini-Soldner projection.
#define GEOGRAPHICLIB_VOLATILE
Definition: Math.hpp:84
Header for GeographicLib::Geodesic class.
static T hypot(T x, T y)
Definition: Math.hpp:255
Math::real Flattening() const
Math::real MajorRadius() const
Header for GeographicLib::Constants class.
void Forward(real lat, real lon, real &x, real &y) const
Geodesic calculations
Definition: Geodesic.hpp:171
Math::real LongitudeOrigin() const