GeographicLib  2.1.2
AuxLatitude.hpp
Go to the documentation of this file.
1 /**
2  * \file AuxLatitude.hpp
3  * \brief Header for the GeographicLib::AuxLatitude and GeographicLib::AuxAngle
4  * classes.
5  *
6  * \note This is just sample code. It is not part of GeographicLib itself.
7  *
8  * This file is an implementation of the methods described in
9  * - C. F. F. Karney,
10  * On auxiliary latitudes,
11  * Technical Report, SRI International, December 2022.
12  * https://arxiv.org/abs/2212.05818
13  * .
14  * Copyright (c) Charles Karney (2022) <charles@karney.com> and licensed under
15  * the MIT/X11 License. For more information, see
16  * https://geographiclib.sourceforge.io/
17  **********************************************************************/
18 
19 #if !defined(AUXLATITUDE_HPP)
20 #define AUXLATITUDE_HPP 1
21 
22 #include <GeographicLib/Math.hpp>
23 
24 #if !defined(GEOGRAPHICLIB_AUXLATITUDE_ORDER)
25 /**
26  * The order of the series approximation used in AuxLatitude.
27  * GEOGRAPHICLIB_AUXLATITUDE_ORDER can be set to one of [4, 6, 8].
28  **********************************************************************/
29 # define GEOGRAPHICLIB_AUXLATITUDE_ORDER 6
30 #endif
31 
32 namespace GeographicLib {
33 
34  /**
35  * \brief An accurate representation of angles.
36  *
37  * \note This is just sample code. It is not part of GeographicLib itself.
38  *
39  * This class is an implementation of the methods described in
40  * - C. F. F. Karney,
41  * On auxiliary latitudes,
42  * Technical Report, SRI International, December 2022.
43  * https://arxiv.org/abs/2212.05818
44  *
45  * An angle is represented be the \e y and \e x coordinates of a point in the
46  * 2d plane. The two coordinates are proportional to the sine and cosine of
47  * the angle. This allows angles close to the cardinal points to be
48  * represented accurately. Only angles in [&minus;180&deg;, 180&deg;] can be
49  * represented. (A possible extension would be to keep count of the number
50  * of turns.)
51  *
52  * @tparam T the floating-point type to use for real numbers.
53  **********************************************************************/
54  template<typename T = double>
55  class AuxAngle {
56  public:
57  /**
58  * The floating-point type for real numbers. This just connects to the
59  * template parameters for the class.
60  **********************************************************************/
61  typedef T real;
62  /**
63  * The constructor.
64  *
65  * @param[in] y the \e y coordinate.
66  * @param[in] x the \e x coordinate.
67  *
68  * \note the \e y coordinate is specified \e first.
69  * \warning either \e x or \e y can be infinite, but not both.
70  *
71  * The defaults (\e x = 1 and \e y = 0) are such that
72  * + no arguments gives an angle of 0;
73  * + 1 argument specifies the tangent of the angle.
74  **********************************************************************/
75  AuxAngle(real y = 0, real x = 1) : _y(y), _x(x) {}
76  /**
77  * @return the \e y component. This is the sine of the angle if the
78  * AuxAngle has been normalized.
79  **********************************************************************/
80  real y() const { return _y; }
81  /**
82  * @return the \e x component. This is the cosine of the angle if the
83  * AuxAngle has been normalized.
84  **********************************************************************/
85  real x() const { return _x; }
86  /**
87  * @return a reference to the \e y component. This allows this component
88  * to be altered.
89  **********************************************************************/
90  real& y() { return _y; }
91  /**
92  * @return a reference to the \e x component. This allows this component
93  * to be altered.
94  **********************************************************************/
95  real& x() { return _x; }
96  /**
97  * @return the AuxAngle converted to the conventional angle measured in
98  * degrees.
99  **********************************************************************/
100  real degrees() const;
101  /**
102  * @return the AuxAngle converted to the conventional angle measured in
103  * radians.
104  **********************************************************************/
105  real radians() const;
106  /**
107  * @return the tangent of the angle.
108  **********************************************************************/
109  real tan() const { return _y / _x; }
110  /**
111  * @return a new normalized AuxAngle with the point lying on the unit
112  * circle and the \e y and \e x components are equal to the sine and
113  * cosine of the angle.
114  **********************************************************************/
115  AuxAngle normalized() const;
116  /**
117  * Normalize the AuxAngle in place so that the \e y and \e x components are
118  * equal to the sine and cosine of the angle.
119  **********************************************************************/
120  void normalize() { *this = normalized(); }
121  /**
122  * Set the quadrant for the AuxAngle.
123  *
124  * @param[in] p the AuxAngle from which the quadrant information is taken.
125  * @return the new AuxAngle in the same quadrant as \e p.
126  **********************************************************************/
127  AuxAngle copyquadrant(const AuxAngle& p) const;
128  /**
129  * Add an AuxAngle.
130  *
131  * @param[in] p the AuxAngle to be added.
132  * @return a reference to the new AuxAngle.
133  *
134  * The addition is done in place, altering the current AuxAngle.
135  *
136  * \warning Neither *this nor \e p should have an infinite component. If
137  * necessary, invoke AuxAngle::normalize on these angles first.
138  **********************************************************************/
139  AuxAngle& operator+=(const AuxAngle& p);
140  /**
141  * Convert degrees to an AuxAngle.
142  *
143  * @param[in] d the angle measured in degrees.
144  * @return the corresponding AuxAngle.
145  *
146  * This allows a new AuxAngle to be initialized as an angle in degrees with
147  * @code
148  * AuxAngle<real> phi = AuxAngle<real>::degrees(d);
149  * @endcode
150  * This is the so-called "named constructor" idiom.
151  **********************************************************************/
152  static AuxAngle degrees(real d);
153  /**
154  * Convert radians to an AuxAngle.
155  *
156  * @param[in] r the angle measured in radians.
157  * @return the corresponding AuxAngle.
158  *
159  * This allows a new AuxAngle to be initialized as an angle in radians with
160  * @code
161  * AuxAngle<real> phi = AuxAngle<real>::radians(r);
162  * @endcode
163  * This is the so-called "named constructor" idiom.
164  **********************************************************************/
165  static AuxAngle radians(real r);
166  /**
167  * @return a "NaN" AuxAngle.
168  **********************************************************************/
169  static AuxAngle NaN();
170  /**
171  * Compute the absolute error in another angle.
172  *
173  * @tparam T1 the floating-point type of the other angle.
174  * @param[in] p the other angle
175  * @return the absolute error between p and *this considered as angles in
176  * radians.
177  **********************************************************************/
178  template<typename T1>
179  real AbsError(const AuxAngle<T1>& p) const;
180  /**
181  * Compute the relative error in another angle.
182  *
183  * @tparam T1 the floating-point type of the other angle.
184  * @param[in] p the other angle
185  * @return the relative error between p.tan() and this->tan().
186  **********************************************************************/
187  template<typename T1>
188  real RelError(const AuxAngle<T1>& p) const;
189  private:
190  real _y, _x;
191  };
192 
193  /// \cond SKIP
194  template<typename T>
196  real y, x;
197  Math::sincosd(d, y, x);
198  return AuxAngle(y, x);
199  }
200 
201  template<typename T>
202  inline AuxAngle<T> AuxAngle<T>::radians(real r) {
203  using std::sin; using std::cos;
204  return AuxAngle(sin(r), cos(r));
205  }
206 
207  template<typename T>
208  inline T AuxAngle<T>::degrees() const {
209  return Math::atan2d(_y, _x);
210  }
211 
212  template<typename T>
213  inline T AuxAngle<T>::radians() const {
214  using std::atan2; return atan2(_y, _x);
215  }
216 
217  template<typename T> template<typename T1>
218  inline T AuxAngle<T>::AbsError(const AuxAngle<T1>& p) const {
219  using std::fabs;
220  return fabs((AuxAngle(-T(p.y()), T(p.x())) += *this).radians());
221  }
222 
223  template<typename T> template<typename T1>
224  inline T AuxAngle<T>::RelError(const AuxAngle<T1>& p) const {
225  using std::fabs;
226  return fabs((T(p.y()) / T(p.x()) - tan()) / tan());
227  }
228  /// \endcond
229 
230  /**
231  * \brief Conversions between auxiliary latitudes.
232  *
233  * \note This is just sample code. It is not part of GeographicLib itself.
234  *
235  * This class is an implementation of the methods described in
236  * - C. F. F. Karney,
237  * On auxiliary latitudes,
238  * Technical Report, SRI International, December 2022.
239  * https://arxiv.org/abs/2212.05818
240  *
241  * The provides accurate conversions between geographic (\e phi, &phi;),
242  * parametric (\e beta, &beta;), geocentric (\e theta, &theta;), rectifying
243  * (\e mu, &mu;), conformal (\e chi, &chi;), and authalic (\e xi, &xi;)
244  * latitudes for an ellipsoid of revolution. A latitude is represented by an
245  * AuxAngle in order to maintain precision close to the poles.
246  *
247  * The class implements two methods for the conversion:
248  * - Direct evaluation of the defining equations, the \e exact method. These
249  * equations are formulated so as to preserve relative accuracy of the
250  * tangent of the latitude, ensuring high accuracy near the equator and the
251  * poles. Newton's method is used for those conversions that can't be
252  * expressed in closed form.
253  * - Expansions in powers of &e n, the third flattening, the \e series
254  * method. This delivers full accuracy for abs(\e f) &le; 1/150. Here, \e
255  * f is the flattening of the ellipsoid.
256  *
257  * The series method is the preferred method of conversion for any conversion
258  * involving &mu;, &chi;, or &xi;, with abs(\e f) &le; 1/150. The equations
259  * for the conversions between &phi;, &beta;, and &theta; are sufficiently
260  * simple that the exact method should be used for such conversions and also
261  * for conversions with with abs(\e f) &gt; 1/150.
262  *
263  * @tparam T the floating-point type to use.
264  *
265  * Example of use:
266  * \include example-AuxLatitude.cpp
267  *
268  * For more information on this projection, see \ref auxlat.
269  **********************************************************************/
270  template<typename T = double>
271  class AuxLatitude {
272  public:
273  /**
274  * The floating-point type for real numbers. This just connects to the
275  * template parameters for the class.
276  **********************************************************************/
277  typedef T real;
278  /**
279  * The type used to represent angles.
280  **********************************************************************/
282  /**
283  * The different auxiliary latitudes.
284  **********************************************************************/
285  enum aux {
286  /**
287  * Geographic latitude, \e phi, &phi;
288  * @hideinitializer
289  **********************************************************************/
291  /**
292  * Parametric latitude, \e beta, &beta;
293  * @hideinitializer
294  **********************************************************************/
296  /**
297  * %Geocentric latitude, \e theta, &theta;
298  * @hideinitializer
299  **********************************************************************/
301  /**
302  * Rectifying latitude, \e mu, &mu;
303  * @hideinitializer
304  **********************************************************************/
306  /**
307  * Conformal latitude, \e chi, &chi;
308  * @hideinitializer
309  **********************************************************************/
311  /**
312  * Authalic latitude, \e xi, &xi;
313  * @hideinitializer
314  **********************************************************************/
315  AUTHALIC = 5,
316  /**
317  * The total number of auxiliary latitudes
318  * @hideinitializer
319  **********************************************************************/
321  /**
322  * An alias for GEOGRAPHIC
323  * @hideinitializer
324  **********************************************************************/
326  /**
327  * An alias for GEOGRAPHIC
328  * @hideinitializer
329  **********************************************************************/
331  /**
332  * An alias for PARAMETRIC
333  * @hideinitializer
334  **********************************************************************/
336  };
337  /**
338  * Constructor
339  *
340  * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
341  * Negative \e f gives a prolate ellipsoid.
342  *
343  * \note the constructor does not precompute the coefficients for the
344  * Fourier series for the series conversions. These are computed and saved
345  * when first needed.
346  **********************************************************************/
347  AuxLatitude(real f);
348  /**
349  * Constructor
350  *
351  * @param[in] a equatorial radius.
352  * @param[in] b polar semi-axis.
353  *
354  * \note the constructor does not precompute the coefficients for the
355  * Fourier series for the series conversions. These are computed and saved
356  * when first needed.
357  **********************************************************************/
358  AuxLatitude(real a, real b);
359  /**
360  * Convert between any two auxiliary latitudes.
361  *
362  * @param[in] auxin an AuxLatitude::aux indicating the type of
363  * auxiliary latitude \e zeta.
364  * @param[in] auxout an AuxLatitude::aux indicating the type of
365  * auxiliary latitude \e eta.
366  * @param[in] zeta the input auxiliary latitude.
367  * @param[in] series if true use the Taylor series instead of the exact
368  * equations [default false].
369  * @return the output auxiliary latitude \e eta.
370  *
371  * With \e series = true, the Fourier coefficients for a specific \e auxin
372  * and \e auxout are computed and saved on the first call; the saved
373  * coefficients are used on subsequent calls. The series method is
374  * accurate for abs(\e f) &le; 1/150; for other \e f, the exact method
375  * should be used
376  **********************************************************************/
377  angle Convert(int auxin, int auxout, const angle& zeta,
378  bool series = false) const;
379  /**
380  * Convert geographic latitude to an auxiliary latitude \e eta.
381  *
382  * @param[in] auxout an AuxLatitude::aux indicating the auxiliary
383  * latitude returned.
384  * @param[in] phi the geographic latitude.
385  * @param[out] diff optional pointer to the derivative d tan(\e eta) / d
386  * tan(\e phi).
387  * @return the auxiliary latitude \e eta.
388  *
389  * This uses the exact equations.
390  **********************************************************************/
391  angle ToAuxiliary(int auxout, const angle& phi,
392  real* diff = nullptr) const;
393  /**
394  * Convert an auxiliary latitude \e zeta to geographic latitude.
395  *
396  * @param[in] auxin an AuxLatitude::aux indicating the type of
397  * auxiliary latitude \e zeta.
398  * @param[in] zeta the input auxiliary latitude.
399  * @param[out] niter optional pointer to the number of iterations.
400  * @return the geographic latitude \e phi.
401  *
402  * This uses the exact equations.
403  **********************************************************************/
404  angle FromAuxiliary(int auxin, const angle& zeta,
405  int* niter = nullptr) const;
406  /**
407  * @return \e f, the flattening of the ellipsoid.
408  **********************************************************************/
409  real Flattening() const { return _f; }
410  /**
411  * The order of the series expansions. This is set at compile time to
412  * either 4, 6, or 8, by the preprocessor macro
413  * GEOGRAPHICLIB_AUXLATITUDE_ORDER.
414  * @hideinitializer
415  **********************************************************************/
417  private:
418  /**
419  * Convert geographic latitude to parametric latitude
420  *
421  * @param[in] phi geographic latitude.
422  * @param[out] diff optional pointer to the derivative d tan(\e beta) / d
423  * tan(\e phi).
424  * @return \e beta, the parametric latitude
425  **********************************************************************/
426  angle Parametric(const angle& phi, real* diff = nullptr) const;
427  /**
428  * Convert geographic latitude to geocentric latitude
429  *
430  * @param[in] phi geographic latitude.
431  * @param[out] diff optional pointer to the derivative d tan(\e theta) / d
432  * tan(\e phi).
433  * @return \e theta, the geocentric latitude.
434  **********************************************************************/
435  angle Geocentric(const angle& phi, real* diff = nullptr) const;
436  /**
437  * Convert geographic latitude to rectifying latitude
438  *
439  * @param[in] phi geographic latitude.
440  * @param[out] diff optional pointer to the derivative d tan(\e mu) / d
441  * tan(\e phi).
442  * @return \e mu, the rectifying latitude.
443  **********************************************************************/
444  angle Rectifying(const angle& phi, real* diff = nullptr) const;
445  /**
446  * Convert geographic latitude to conformal latitude
447  *
448  * @param[in] phi geographic latitude.
449  * @param[out] diff optional pointer to the derivative d tan(\e chi) / d
450  * tan(\e phi).
451  * @return \e chi, the conformal latitude.
452  **********************************************************************/
453  angle Conformal(const angle& phi, real* diff = nullptr) const;
454  /**
455  * Convert geographic latitude to authalic latitude
456  *
457  * @param[in] phi geographic latitude.
458  * @param[out] diff optional pointer to the derivative d tan(\e xi) / d
459  * tan(\e phi).
460  * @return \e xi, the authalic latitude.
461  **********************************************************************/
462  angle Authalic(const angle& phi, real* diff = nullptr) const;
463  // Maximum number of iterations for Newton's method
464  static const int numit_ = 1000;
465  real tol_, bmin_, bmax_; // Static consts for Newton's method
466  // Ellipsoid parameters
467  real _f, _fm1, _e2, _e2m1, _e12, _e12p1, _n, _e, _e1, _n2, _q;
468  // To hold computed Fourier coefficients
469  mutable real _c[Lmax * AUXNUMBER * AUXNUMBER];
470  // 1d index into AUXNUMBER x AUXNUMBER data
471  static int ind(int auxout, int auxin) {
472  return (auxout >= 0 && auxout < AUXNUMBER &&
473  auxin >= 0 && auxin < AUXNUMBER) ?
474  AUXNUMBER * auxout + auxin : -1;
475  }
476  // the function sqrt(1 + tphi^2), convert tan to sec
477  static real sc(real tphi)
478  { using std::hypot; return hypot(real(1), tphi); }
479  // the function tphi / sqrt(1 + tphi^2), convert tan to sin
480  static real sn(real tphi) {
481  using std::isfinite; using std::isnan; using std::copysign;
482  return isfinite(tphi) || isnan(tphi) ? tphi / sc(tphi) :
483  copysign(real(1), tphi);
484  }
485  // The symmetric elliptic integral RD
486  static real RD(real x, real y, real z);
487  // The symmetric elliptic integral RF
488  static real RF(real x, real y, real z);
489  // the function atanh(e * sphi)/e; works for e^2 = 0 and e^2 < 0
490  real atanhee(real tphi) const;
491  // the function atanh(e * sphi)/e + sphi / (1 - (e * sphi)^2);
492  real q(real tphi) const;
493  // The divided difference of (q(1) - q(sphi)) / (1 - sphi)
494  real Dq(real tphi) const;
495  // Populate [_c[Lmax * k], _c[Lmax * (k + 1)])
496  void fillcoeff(int auxin, int auxout, int k) const;
497  // Clenshaw applied to sum(c[k] * sin( (2*k+2) * zeta), i, 0, K-1)
498  // if alt, use the Reinsch optimizations
499  static real Clenshaw(real szeta, real czeta, const real c[], int K,
500  bool alt = true);
501  };
502 
503 } // namespace GeographicLib
504 
505 #endif // AUXLATITUDE_HPP
#define GEOGRAPHICLIB_AUXLATITUDE_ORDER
Definition: AuxLatitude.hpp:29
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
Header for GeographicLib::Math class.
An accurate representation of angles.
Definition: AuxLatitude.hpp:55
AuxAngle & operator+=(const AuxAngle &p)
Definition: AuxLatitude.cpp:79
static AuxAngle radians(real r)
AuxAngle(real y=0, real x=1)
Definition: AuxLatitude.hpp:75
real RelError(const AuxAngle< T1 > &p) const
AuxAngle copyquadrant(const AuxAngle &p) const
Definition: AuxLatitude.cpp:74
real AbsError(const AuxAngle< T1 > &p) const
static AuxAngle degrees(real d)
static AuxAngle NaN()
Definition: AuxLatitude.cpp:52
AuxAngle normalized() const
Definition: AuxLatitude.cpp:58
Conversions between auxiliary latitudes.
angle FromAuxiliary(int auxin, const angle &zeta, int *niter=nullptr) const
angle ToAuxiliary(int auxout, const angle &phi, real *diff=nullptr) const
angle Convert(int auxin, int auxout, const angle &zeta, bool series=false) const
Geocentric coordinates
Definition: Geocentric.hpp:67
static void sincosd(T x, T &sinx, T &cosx)
Definition: Math.cpp:106
static T atan2d(T y, T x)
Definition: Math.cpp:183
Namespace for GeographicLib.
Definition: Accumulator.cpp:12