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>
20
#include <
GeographicLib/Math.hpp
>
21
#include <
GeographicLib/AuxAngle.hpp
>
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
35
namespace
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, φ),
48
* parametric (\e beta, β), geocentric (\e theta, θ), rectifying
49
* (\e mu, μ), conformal (\e chi, χ), and authalic (\e 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) ≤ 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 μ, χ, or ξ, with abs(\e f) ≤ 1/150. The equations
65
* for the conversions between φ, β, and θ are sufficiently
66
* simple that the exact method should be used for such conversions and also
67
* for conversions with with abs(\e f) > 1/150.
68
*
69
* Example of use:
70
* \include example-AuxLatitude.cpp
71
**********************************************************************/
72
class
GEOGRAPHICLIB_EXPORT
AuxLatitude
{
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, φ
86
* @hideinitializer
87
**********************************************************************/
88
GEOGRAPHIC = 0,
89
/**
90
* Parametric latitude, \e beta, β
91
* @hideinitializer
92
**********************************************************************/
93
PARAMETRIC = 1,
94
/**
95
* %Geocentric latitude, \e theta, θ
96
* @hideinitializer
97
**********************************************************************/
98
GEOCENTRIC = 2,
99
/**
100
* Rectifying latitude, \e mu, μ
101
* @hideinitializer
102
**********************************************************************/
103
RECTIFYING = 3,
104
/**
105
* Conformal latitude, \e chi, χ
106
* @hideinitializer
107
**********************************************************************/
108
CONFORMAL = 4,
109
/**
110
* Authalic latitude, \e 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 − \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) ≤ 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) ≤ 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
AuxAngle.hpp
Header for the GeographicLib::AuxAngle class.
GEOGRAPHICLIB_AUXLATITUDE_ORDER
#define GEOGRAPHICLIB_AUXLATITUDE_ORDER
Definition
AuxLatitude.hpp:30
GEOGRAPHICLIB_EXPORT
#define GEOGRAPHICLIB_EXPORT
Definition
Constants.hpp:67
real
GeographicLib::Math::real real
Definition
GeodSolve.cpp:31
Math.hpp
Header for GeographicLib::Math class.
GeographicLib::AuxAngle
An accurate representation of angles.
Definition
AuxAngle.hpp:43
GeographicLib::AuxLatitude
Conversions between auxiliary latitudes.
Definition
AuxLatitude.hpp:72
GeographicLib::AuxLatitude::EquatorialRadius
Math::real EquatorialRadius() const
Definition
AuxLatitude.hpp:281
GeographicLib::AuxLatitude::Flattening
Math::real Flattening() const
Definition
AuxLatitude.hpp:289
GeographicLib::AuxLatitude::aux
aux
Definition
AuxLatitude.hpp:83
GeographicLib::AuxLatitude::PolarSemiAxis
Math::real PolarSemiAxis() const
Definition
AuxLatitude.hpp:285
GeographicLib::AuxLatitude::axes
static AuxLatitude axes(real a, real b)
Definition
AuxLatitude.hpp:192
GeographicLib::AuxLatitude::Clenshaw
static Math::real Clenshaw(bool sinp, real szeta, real czeta, const real c[], int K)
GeographicLib::Geocentric
Geocentric coordinates
Definition
Geocentric.hpp:67
GeographicLib::Math::real
double real
Definition
Math.hpp:109
GeographicLib
Namespace for GeographicLib.
Definition
Accumulator.cpp:12
include
GeographicLib
AuxLatitude.hpp
Generated by
1.9.8