GeographicLib 2.2
AuxLatitude.hpp
Go to the documentation of this file.
1/**
2 * \file AuxLatitude.hpp
3 * \brief Header for the GeographicLib::AuxLatitude class
4 *
5 * This file is an implementation of the methods described in
6 * - C. F. F. Karney,
7 * On auxiliary latitudes,
8 * Technical Report, SRI International, December 2022.
9 * https://arxiv.org/abs/2212.05818
10 * .
11 * Copyright (c) Charles Karney (2022-2023) <charles@karney.com> and licensed
12 * under the MIT/X11 License. For more information, see
13 * https://geographiclib.sourceforge.io/
14 **********************************************************************/
15
16#if !defined(GEOGRAPHICLIB_AUXLATITUDE_HPP)
17#define GEOGRAPHICLIB_AUXLATITUDE_HPP 1
18
19#include <utility>
22
23#if !defined(GEOGRAPHICLIB_AUXLATITUDE_ORDER)
24/**
25 * The order of the series approximation used in AuxLatitude.
26 * GEOGRAPHICLIB_AUXLATITUDE_ORDER can be set to one of [4, 6, 8]. Use order
27 * appropriate for double precision, 6, also for GEOGRAPHICLIB_PRECISION == 5
28 * to enable truncation errors to be measured easily.
29 **********************************************************************/
30# define GEOGRAPHICLIB_AUXLATITUDE_ORDER \
31 (GEOGRAPHICLIB_PRECISION == 2 || GEOGRAPHICLIB_PRECISION == 5 ? 6 : \
32 (GEOGRAPHICLIB_PRECISION == 1 ? 4 : 8))
33#endif
34
35namespace GeographicLib {
36
37 /**
38 * \brief Conversions between auxiliary latitudes.
39 *
40 * This class is an implementation of the methods described in
41 * - C. F. F. Karney,<br>
42 * <a href="https://arxiv.org/abs/2212.05818">On
43 * auxiliary latitudes</a>,<br>
44 * Technical Report, SRI International, December 2023.<br>
45 * <a href="https://arxiv.org/abs/2212.05818">arxiv:2212.05818</a>
46 *
47 * The provides accurate conversions between geographic (\e phi, &phi;),
48 * parametric (\e beta, &beta;), geocentric (\e theta, &theta;), rectifying
49 * (\e mu, &mu;), conformal (\e chi, &chi;), and authalic (\e xi, &xi;)
50 * latitudes for an ellipsoid of revolution. A latitude is represented by
51 * the class AuxAngle in order to maintain precision close to the poles.
52 *
53 * The class implements two methods for the conversion:
54 * - Direct evaluation of the defining equations, the \e exact method. These
55 * equations are formulated so as to preserve relative accuracy of the
56 * tangent of the latitude, ensuring high accuracy near the equator and the
57 * poles. Newton's method is used for those conversions that can't be
58 * expressed in closed form.
59 * - Expansions in powers of \e n, the third flattening, the \e series
60 * method. This delivers full accuracy for abs(\e f) &le; 1/150. Here, \e
61 * f is the flattening of the ellipsoid.
62 *
63 * The series method is the preferred method of conversion for any conversion
64 * involving &mu;, &chi;, or &xi;, with abs(\e f) &le; 1/150. The equations
65 * for the conversions between &phi;, &beta;, and &theta; are sufficiently
66 * simple that the exact method should be used for such conversions and also
67 * for conversions with with abs(\e f) &gt; 1/150.
68 *
69 * Example of use:
70 * \include example-AuxLatitude.cpp
71 **********************************************************************/
73 typedef Math::real real;
74 AuxLatitude(const std::pair<real, real>& axes);
75 public:
76 /**
77 * The floating-point type for real numbers. This just connects to the
78 * template parameters for the class.
79 **********************************************************************/
80 /**
81 * The different auxiliary latitudes.
82 **********************************************************************/
83 enum aux {
84 /**
85 * Geographic latitude, \e phi, &phi;
86 * @hideinitializer
87 **********************************************************************/
88 GEOGRAPHIC = 0,
89 /**
90 * Parametric latitude, \e beta, &beta;
91 * @hideinitializer
92 **********************************************************************/
93 PARAMETRIC = 1,
94 /**
95 * %Geocentric latitude, \e theta, &theta;
96 * @hideinitializer
97 **********************************************************************/
98 GEOCENTRIC = 2,
99 /**
100 * Rectifying latitude, \e mu, &mu;
101 * @hideinitializer
102 **********************************************************************/
103 RECTIFYING = 3,
104 /**
105 * Conformal latitude, \e chi, &chi;
106 * @hideinitializer
107 **********************************************************************/
108 CONFORMAL = 4,
109 /**
110 * Authalic latitude, \e xi, &xi;
111 * @hideinitializer
112 **********************************************************************/
113 AUTHALIC = 5,
114 /**
115 * The total number of auxiliary latitudes
116 * @hideinitializer
117 **********************************************************************/
118 AUXNUMBER = 6,
119 /**
120 * An alias for GEOGRAPHIC
121 * @hideinitializer
122 **********************************************************************/
123 PHI = GEOGRAPHIC,
124 /**
125 * An alias for PARAMETRIC
126 * @hideinitializer
127 **********************************************************************/
128 BETA = PARAMETRIC,
129 /**
130 * An alias for GEOCENTRIC
131 * @hideinitializer
132 **********************************************************************/
133 THETA = GEOCENTRIC,
134 /**
135 * An alias for RECTIFYING
136 * @hideinitializer
137 **********************************************************************/
138 MU = RECTIFYING,
139 /**
140 * An alias for CONFORMAL
141 * @hideinitializer
142 **********************************************************************/
143 CHI = CONFORMAL,
144 /**
145 * An alias for AUTHALIC
146 * @hideinitializer
147 **********************************************************************/
148 XI = AUTHALIC,
149 /**
150 * An alias for GEOGRAPHIC
151 * @hideinitializer
152 **********************************************************************/
153 COMMON = GEOGRAPHIC,
154 /**
155 * An alias for GEOGRAPHIC
156 * @hideinitializer
157 **********************************************************************/
158 GEODETIC = GEOGRAPHIC,
159 /**
160 * An alias for PARAMETRIC
161 * @hideinitializer
162 **********************************************************************/
163 REDUCED = PARAMETRIC,
164 };
165 /**
166 * Constructor
167 *
168 * @param[in] a equatorial radius.
169 * @param[in] f flattening of ellipsoid. Setting \e f = 0 gives a sphere.
170 * Negative \e f gives a prolate ellipsoid.
171 * @exception GeographicErr if \e a or (1 &minus; \e f) \e a is not
172 * positive.
173 *
174 * \note the constructor does not precompute the coefficients for the
175 * Fourier series for the series conversions. These are computed and saved
176 * when first needed.
177 **********************************************************************/
178 AuxLatitude(real a, real f);
179 /**
180 * Construct and return an AuxLatitude object specified in terms of the
181 * semi-axes
182 *
183 * @param[in] a equatorial radius.
184 * @param[in] b polar semi-axis.
185 * @exception GeographicErr if \e a or \e b is not positive.
186 *
187 * This allows a new AuxAngle to be initialized as an angle in radians with
188 * @code
189 * AuxLatitude aux(AuxLatitude::axes(a, b));
190 * @endcode
191 **********************************************************************/
192 static AuxLatitude axes(real a, real b) {
193 return AuxLatitude(std::pair<real, real>(a, b));
194 }
195 /**
196 * Convert between any two auxiliary latitudes specified as AuxAngle.
197 *
198 * @param[in] auxin an AuxLatitude::aux indicating the type of
199 * auxiliary latitude \e zeta.
200 * @param[in] auxout an AuxLatitude::aux indicating the type of
201 * auxiliary latitude \e eta.
202 * @param[in] zeta the input auxiliary latitude as an AuxAngle
203 * @param[in] exact if true use the exact equations instead of the Taylor
204 * series [default false].
205 * @return the output auxiliary latitude \e eta as an AuxAngle.
206 *
207 * With \e exact = false, the Fourier coefficients for a specific \e auxin
208 * and \e auxout are computed and saved on the first call; the saved
209 * coefficients are used on subsequent calls. The series method is
210 * accurate for abs(\e f) &le; 1/150; for other \e f, the exact method
211 * should be used.
212 **********************************************************************/
213 AuxAngle Convert(int auxin, int auxout, const AuxAngle& zeta,
214 bool exact = false) const;
215 /**
216 * Convert between any two auxiliary latitudes specified in degrees.
217 *
218 * @param[in] auxin an AuxLatitude::aux indicating the type of
219 * auxiliary latitude \e zeta.
220 * @param[in] auxout an AuxLatitude::aux indicating the type of
221 * auxiliary latitude \e eta.
222 * @param[in] zeta the input auxiliary latitude in degrees.
223 * @param[in] exact if true use the exact equations instead of the Taylor
224 * series [default false].
225 * @return the output auxiliary latitude \e eta in degrees.
226 *
227 * With \e exact = false, the Fourier coefficients for a specific \e auxin
228 * and \e auxout are computed and saved on the first call; the saved
229 * coefficients are used on subsequent calls. The series method is
230 * accurate for abs(\e f) &le; 1/150; for other \e f, the exact method
231 * should be used.
232 **********************************************************************/
233 Math::real Convert(int auxin, int auxout, real zeta, bool exact = false)
234 const;
235 /**
236 * Convert geographic latitude to an auxiliary latitude \e eta.
237 *
238 * @param[in] auxout an AuxLatitude::aux indicating the auxiliary
239 * latitude returned.
240 * @param[in] phi the geographic latitude.
241 * @param[out] diff optional pointer to the derivative d tan(\e eta) / d
242 * tan(\e phi).
243 * @return the auxiliary latitude \e eta.
244 *
245 * This uses the exact equations.
246 **********************************************************************/
247 AuxAngle ToAuxiliary(int auxout, const AuxAngle& phi, real* diff = nullptr)
248 const;
249 /**
250 * Convert an auxiliary latitude \e zeta to geographic latitude.
251 *
252 * @param[in] auxin an AuxLatitude::aux indicating the type of
253 * auxiliary latitude \e zeta.
254 * @param[in] zeta the input auxiliary latitude.
255 * @param[out] niter optional pointer to the number of iterations.
256 * @return the geographic latitude \e phi.
257 *
258 * This uses the exact equations.
259 **********************************************************************/
260 AuxAngle FromAuxiliary(int auxin, const AuxAngle& zeta,
261 int* niter = nullptr) const;
262 /**
263 * Return the rectifying radius.
264 *
265 * @param[in] exact if true use the exact expression instead of the Taylor
266 * series [default false].
267 * @return the rectifying radius in the same units as \e a.
268 **********************************************************************/
269 Math::real RectifyingRadius(bool exact = false) const;
270 /**
271 * Return the authalic radius squared.
272 *
273 * @param[in] exact if true use the exact expression instead of the Taylor
274 * series [default false].
275 * @return the authalic radius squared in the same units as \e a.
276 **********************************************************************/
277 Math::real AuthalicRadiusSquared(bool exact = false) const;
278 /**
279 * @return \e a the equatorial radius of the ellipsoid (meters).
280 **********************************************************************/
281 Math::real EquatorialRadius() const { return _a; }
282 /**
283 * @return \e b the polar semi-axis of the ellipsoid (meters).
284 **********************************************************************/
285 Math::real PolarSemiAxis() const { return _b; }
286 /**
287 * @return \e f, the flattening of the ellipsoid.
288 **********************************************************************/
289 Math::real Flattening() const { return _f; }
290 /**
291 * Use Clenshaw to sum a Fouier series.
292 *
293 * @param[in] sinp if true sum the sine series, else sum the cosine series.
294 * @param[in] szeta sin(\e zeta).
295 * @param[in] czeta cos(\e zeta).
296 * @param[in] c the array of coefficients.
297 * @param[in] K the number of coefficients.
298 * @return the Clenshaw sum.
299 *
300 * The result returned is \f$ \sum_0^{K-1} c_k \sin (2k+2)\zeta \f$, if \e
301 * sinp is true; with \e sinp false, replace sin by cos.
302 **********************************************************************/
303 // Clenshaw applied to sum(c[k] * sin( (2*k+2) * zeta), i, 0, K-1);
304 // if !sinp then subst sine->cosine.
305 static Math::real Clenshaw(bool sinp, real szeta, real czeta,
306 const real c[], int K);
307 /**
308 * The order of the series expansions. This is set at compile time to
309 * either 4, 6, or 8, by the preprocessor macro
310 * GEOGRAPHICLIB_AUXLATITUDE_ORDER.
311 * @hideinitializer
312 **********************************************************************/
313 static const int Lmax = GEOGRAPHICLIB_AUXLATITUDE_ORDER;
314 /**
315 * A global instantiation of Ellipsoid with the parameters for the WGS84
316 * ellipsoid.
317 **********************************************************************/
318 static const AuxLatitude& WGS84();
319 private:
320 // Maximum number of iterations for Newton's method
321 static const int numit_ = 1000;
322 real tol_, bmin_, bmax_; // Static consts for Newton's method
323 // the function atanh(e * sphi)/e + sphi / (1 - (e * sphi)^2);
324 protected:
325 /**
326 * Convert geographic latitude to parametric latitude
327 *
328 * @param[in] phi geographic latitude.
329 * @param[out] diff optional pointer to the derivative d tan(\e beta) / d
330 * tan(\e phi).
331 * @return \e beta, the parametric latitude
332 **********************************************************************/
333 AuxAngle Parametric(const AuxAngle& phi, real* diff = nullptr) const;
334 /**
335 * Convert geographic latitude to geocentric latitude
336 *
337 * @param[in] phi geographic latitude.
338 * @param[out] diff optional pointer to the derivative d tan(\e theta) / d
339 * tan(\e phi).
340 * @return \e theta, the geocentric latitude.
341 **********************************************************************/
342 AuxAngle Geocentric(const AuxAngle& phi, real* diff = nullptr) const;
343 /**
344 * Convert geographic latitude to rectifying latitude
345 *
346 * @param[in] phi geographic latitude.
347 * @param[out] diff optional pointer to the derivative d tan(\e mu) / d
348 * tan(\e phi).
349 * @return \e mu, the rectifying latitude.
350 **********************************************************************/
351 AuxAngle Rectifying(const AuxAngle& phi, real* diff = nullptr) const;
352 /**
353 * Convert geographic latitude to conformal latitude
354 *
355 * @param[in] phi geographic latitude.
356 * @param[out] diff optional pointer to the derivative d tan(\e chi) / d
357 * tan(\e phi).
358 * @return \e chi, the conformal latitude.
359 **********************************************************************/
360 AuxAngle Conformal(const AuxAngle& phi, real* diff = nullptr) const;
361 /**
362 * Convert geographic latitude to authalic latitude
363 *
364 * @param[in] phi geographic latitude.
365 * @param[out] diff optional pointer to the derivative d tan(\e xi) / d
366 * tan(\e phi).
367 * @return \e xi, the authalic latitude.
368 **********************************************************************/
369 AuxAngle Authalic(const AuxAngle& phi, real* diff = nullptr) const;
370 /// \cond SKIP
371 // Ellipsoid parameters
372 real _a, _b, _f, _fm1, _e2, _e2m1, _e12, _e12p1, _n, _e, _e1, _n2, _q;
373 // To hold computed Fourier coefficients
374 mutable real _c[Lmax * AUXNUMBER * AUXNUMBER];
375 // 1d index into AUXNUMBER x AUXNUMBER data
376 static int ind(int auxout, int auxin) {
377 return (auxout >= 0 && auxout < AUXNUMBER &&
378 auxin >= 0 && auxin < AUXNUMBER) ?
379 AUXNUMBER * auxout + auxin : -1;
380 }
381 // the function sqrt(1 + tphi^2), convert tan to sec
382 static real sc(real tphi)
383 { using std::hypot; return hypot(real(1), tphi); }
384 // the function tphi / sqrt(1 + tphi^2), convert tan to sin
385 static real sn(real tphi) {
386 using std::isinf; using std::copysign;
387 return isinf(tphi) ? copysign(real(1), tphi) : tphi / sc(tphi);
388 }
389 // Populate [_c[Lmax * k], _c[Lmax * (k + 1)])
390 void fillcoeff(int auxin, int auxout, int k) const;
391 // the function atanh(e * sphi)/e; works for e^2 = 0 and e^2 < 0
392 real atanhee(real tphi) const;
393 /// \endcond
394 private:
395 // the function atanh(e * sphi)/e + sphi / (1 - (e * sphi)^2);
396 real q(real tphi) const;
397 // The divided difference of (q(1) - q(sphi)) / (1 - sphi)
398 real Dq(real tphi) const;
399 };
400
401} // namespace GeographicLib
402
403#endif // GEOGRAPHICLIB_AUXLATITUDE_HPP
Header for the GeographicLib::AuxAngle class.
#define GEOGRAPHICLIB_AUXLATITUDE_ORDER
Definition: AuxLatitude.hpp:30
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:67
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
Header for GeographicLib::Math class.
An accurate representation of angles.
Definition: AuxAngle.hpp:43
Conversions between auxiliary latitudes.
Definition: AuxLatitude.hpp:72
Math::real EquatorialRadius() const
Math::real Flattening() const
Math::real PolarSemiAxis() const
static AuxLatitude axes(real a, real b)
static Math::real Clenshaw(bool sinp, real szeta, real czeta, const real c[], int K)
Geocentric coordinates
Definition: Geocentric.hpp:67
Namespace for GeographicLib.
Definition: Accumulator.cpp:12