GeographicLib  1.35
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  Geodesic _earth;
73  GeodesicLine _meridian;
74  real _sbet0, _cbet0;
75  static const real eps1_;
76  static const real tiny_;
77  static const unsigned maxit_ = 10;
78 
79  // The following private helper functions are copied from Geodesic.
80  static inline real AngRound(real x) throw() {
81  // The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57
82  // for reals = 0.7 pm on the earth if x is an angle in degrees. (This
83  // is about 1000 times more resolution than we get with angles around 90
84  // degrees.) We use this to avoid having to deal with near singular
85  // cases when x is non-zero but tiny (e.g., 1.0e-200).
86  const real z = 1/real(16);
87  volatile real y = std::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) throw() {
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) throw()
109  : _earth(earth) {}
110 
111  /**
112  * Constructor for CassiniSoldner specifying a center point.
113  *
114  * @param[in] lat0 latitude of center point of projection (degrees).
115  * @param[in] lon0 longitude of center point of projection (degrees).
116  * @param[in] earth the Geodesic object to use for geodesic calculations.
117  * By default this uses the WGS84 ellipsoid.
118  *
119  * \e lat0 should be in the range [&minus;90&deg;, 90&deg;] and \e
120  * lon0 should be in the range [&minus;540&deg;, 540&deg;).
121  **********************************************************************/
122  CassiniSoldner(real lat0, real lon0,
123  const Geodesic& earth = Geodesic::WGS84) throw()
124  : _earth(earth) {
125  Reset(lat0, lon0);
126  }
127 
128  /**
129  * Set the central point of the projection
130  *
131  * @param[in] lat0 latitude of center point of projection (degrees).
132  * @param[in] lon0 longitude of center point of projection (degrees).
133  *
134  * \e lat0 should be in the range [&minus;90&deg;, 90&deg;] and \e
135  * lon0 should be in the range [&minus;540&deg;, 540&deg;).
136  **********************************************************************/
137  void Reset(real lat0, real lon0) throw();
138 
139  /**
140  * Forward projection, from geographic to Cassini-Soldner.
141  *
142  * @param[in] lat latitude of point (degrees).
143  * @param[in] lon longitude of point (degrees).
144  * @param[out] x easting of point (meters).
145  * @param[out] y northing of point (meters).
146  * @param[out] azi azimuth of easting direction at point (degrees).
147  * @param[out] rk reciprocal of azimuthal northing scale at point.
148  *
149  * \e lat should be in the range [&minus;90&deg;, 90&deg;] and \e
150  * lon should be in the range [&minus;540&deg;, 540&deg;). A call
151  * to Forward followed by a call to Reverse will return the original (\e
152  * lat, \e lon) (to within roundoff). The routine does nothing if the
153  * origin has not been set.
154  **********************************************************************/
155  void Forward(real lat, real lon,
156  real& x, real& y, real& azi, real& rk) const throw();
157 
158  /**
159  * Reverse projection, from Cassini-Soldner to geographic.
160  *
161  * @param[in] x easting of point (meters).
162  * @param[in] y northing of point (meters).
163  * @param[out] lat latitude of point (degrees).
164  * @param[out] lon longitude of point (degrees).
165  * @param[out] azi azimuth of easting direction at point (degrees).
166  * @param[out] rk reciprocal of azimuthal northing scale at point.
167  *
168  * A call to Reverse followed by a call to Forward will return the original
169  * (\e x, \e y) (to within roundoff), provided that \e x and \e y are
170  * sufficiently small not to "wrap around" the earth. The routine does
171  * nothing if the origin has not been set.
172  **********************************************************************/
173  void Reverse(real x, real y,
174  real& lat, real& lon, real& azi, real& rk) const throw();
175 
176  /**
177  * CassiniSoldner::Forward without returning the azimuth and scale.
178  **********************************************************************/
179  void Forward(real lat, real lon,
180  real& x, real& y) const throw() {
181  real azi, rk;
182  Forward(lat, lon, x, y, azi, rk);
183  }
184 
185  /**
186  * CassiniSoldner::Reverse without returning the azimuth and scale.
187  **********************************************************************/
188  void Reverse(real x, real y,
189  real& lat, real& lon) const throw() {
190  real azi, rk;
191  Reverse(x, y, lat, lon, azi, rk);
192  }
193 
194  /** \name Inspector functions
195  **********************************************************************/
196  ///@{
197  /**
198  * @return true if the object has been initialized.
199  **********************************************************************/
200  bool Init() const throw() { return _meridian.Init(); }
201 
202  /**
203  * @return \e lat0 the latitude of origin (degrees).
204  **********************************************************************/
205  Math::real LatitudeOrigin() const throw()
206  { return _meridian.Latitude(); }
207 
208  /**
209  * @return \e lon0 the longitude of origin (degrees).
210  **********************************************************************/
211  Math::real LongitudeOrigin() const throw()
212  { return _meridian.Longitude(); }
213 
214  /**
215  * @return \e a the equatorial radius of the ellipsoid (meters). This is
216  * the value inherited from the Geodesic object used in the constructor.
217  **********************************************************************/
218  Math::real MajorRadius() const throw() { return _earth.MajorRadius(); }
219 
220  /**
221  * @return \e f the flattening of the ellipsoid. This is the value
222  * inherited from the Geodesic object used in the constructor.
223  **********************************************************************/
224  Math::real Flattening() const throw() { return _earth.Flattening(); }
225  ///@}
226 
227  /// \cond SKIP
228  /**
229  * <b>DEPRECATED</b>
230  * @return \e r the inverse flattening of the ellipsoid.
231  **********************************************************************/
232  Math::real InverseFlattening() const throw()
233  { return _earth.InverseFlattening(); }
234  /// \endcond
235  };
236 
237 } // namespace GeographicLib
238 
239 #endif // GEOGRAPHICLIB_CASSINISOLDNER_HPP