GeographicLib 2.1.2
AuxLatitude.cpp
Go to the documentation of this file.
1/**
2 * \file AuxLatitude.cpp
3 * \brief Implementation for the GeographicLib::AuxLatitude and
4 * GeographicLib::AuxAngle 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#include "AuxLatitude.hpp"
20
21/// \cond SKIP
22
23#if !defined(AUXLATITUDE_UNOPT)
24// If 1 use the unoptimized formulas in the paper. Otherwise use the
25// optimizations to improve roundoff for high eccentricity.
26#define AUXLATITUDE_UNOPT 0
27#endif
28
29#if !defined(AUXLATITUDE_NATIVE_ELLIPTIC)
30// If 1 use the templated implementation of RD and RF in this file. Otherwise
31// use the non-templated versions in the EllipticFunction class.
32#define AUXLATITUDE_NATIVE_ELLIPTIC 1
33#endif
34
35/// \endcond
36
37#if !AUXLATITUDE_NATIVE_ELLIPTIC
38// Use GeographicLib::EllipticFunction class
39#include <GeographicLib/EllipticFunction.hpp>
40#endif
41
42#if defined(_MSC_VER)
43// Squelch warnings about constant conditional expressions
44# pragma warning (disable: 4127)
45#endif
46
47namespace GeographicLib {
48
49 using namespace std;
50
51 template<typename T>
52 AuxAngle<T> AuxAngle<T>::NaN() {
53 return AuxAngle(std::numeric_limits<real>::quiet_NaN(),
54 std::numeric_limits<real>::quiet_NaN());
55 }
56
57 template<typename T>
58 AuxAngle<T> AuxAngle<T>::normalized() const {
59 if ( isnan( tan() ) ||
60 (fabs(_y) > std::numeric_limits<real>::max()/2 &&
61 fabs(_x) > std::numeric_limits<real>::max()/2) )
62 // deal with
63 // (0,0), (inf,inf), (nan,nan), (nan,x), (y,nan), (toobig,toobig)
64 return NaN();
65 real r = hypot(_y, _x),
66 y = _y/r, x = _x/r;
67 // deal with r = inf, then one of y,x
68 if (isnan(y)) y = copysign(real(1), _y);
69 if (isnan(x)) x = copysign(real(1), _x);
70 return AuxAngle(y, x);
71 }
72
73 template<typename T>
74 AuxAngle<T> AuxAngle<T>::copyquadrant(const AuxAngle& p) const {
75 return AuxAngle(copysign(y(), p.y()), copysign(x(), p.x()));
76 }
77
78 template<typename T>
79 AuxAngle<T>& AuxAngle<T>::operator+=(const AuxAngle& p) {
80 // Do nothing if p.tan() == 0 to preserve signs of y() and x()
81 if (p.tan() != 0) {
82 real x = _x * p._x - _y * p._y;
83 _y = _y * p._x + _x * p._y;
84 _x = x;
85 }
86 return *this;
87 }
88
89 template<typename T>
90 AuxLatitude<T>::AuxLatitude(real f)
91 : tol_( sqrt(numeric_limits<real>::epsilon()) )
92#if AUXLATITUDE_UNOPT
93 , bmin_( numeric_limits<real>::min() )
94 , bmax_( numeric_limits<real>::max() )
95#else
96 , bmin_( log2(numeric_limits<real>::min()) )
97 , bmax_( log2(numeric_limits<real>::max()) )
98#endif
99 , _f( f )
100 , _fm1( 1 - _f )
101 , _e2( _f * (2 - _f) )
102 , _e2m1( _fm1 * _fm1 )
103 , _e12( _e2/(1 - _e2) )
104 , _e12p1( 1 / _e2m1 )
105 , _n( _f/(2 - _f) )
106 , _e( sqrt(fabs(_e2)) )
107 , _e1( sqrt(fabs(_e12)) )
108 , _n2( _n * _n )
109#if AUXLATITUDE_UNOPT
110 , _q( _e12p1 + (_f == 0 ? 1 : (_f > 0 ? asinh(_e1) : atan(_e)) / _e) )
111#else
112 , _q( _e12p1 + (_f == 0 ? 1 : (_f > 0 ? atanh(_e) : atan(_e)) / _e) )
113#endif
114 {
115 if (!(isfinite(_f) && _f < 1 &&
116 isfinite(_n) && fabs(_n) < 1))
117 throw GeographicErr("Bad ellipsoid parameters");
118 fill(_c, _c + Lmax * AUXNUMBER * AUXNUMBER,
119 numeric_limits<real>::quiet_NaN());
120 }
121
122 template<typename T>
123 AuxLatitude<T>::AuxLatitude(real a, real b)
124 : tol_( sqrt(numeric_limits<real>::epsilon()) )
125#if AUXLATITUDE_UNOPT
126 , bmin_( numeric_limits<real>::min() )
127 , bmax_( numeric_limits<real>::max() )
128#else
129 , bmin_( log2(numeric_limits<real>::min()) )
130 , bmax_( log2(numeric_limits<real>::max()) )
131#endif
132 , _f( (a - b) / a )
133 , _fm1( b / a )
134 , _e2( ((a - b) * (a + b)) / (a * a) )
135 , _e2m1( (b * b) / (a * a) )
136 , _e12( ((a - b) * (a + b)) / (b * b) )
137 , _e12p1( (a * a) / (b * b) )
138 , _n( (a - b) / (a + b) )
139 , _e( sqrt(fabs(a - b) * (a + b)) / a )
140 , _e1( sqrt(fabs(a - b) * (a + b)) / b )
141 , _n2( _n * _n )
142#if AUXLATITUDE_UNOPT
143 , _q( _e12p1 + (_f == 0 ? 1 : (_f > 0 ? asinh(_e1) : atan(_e)) / _e) )
144#else
145 , _q( _e12p1 + (_f == 0 ? 1 : (_f > 0 ? atanh(_e) : atan(_e)) / _e) )
146#endif
147 {
148 if (!(isfinite(_f) && _f < 1 &&
149 isfinite(_n) && fabs(_n) < 1))
150 throw GeographicErr("Bad ellipsoid parameters");
151 fill(_c, _c + Lmax * AUXNUMBER * AUXNUMBER,
152 numeric_limits<real>::quiet_NaN());
153 }
154
155 template<typename T>
156 AuxAngle<T> AuxLatitude<T>::Parametric(const angle& phi, real* diff) const {
157 if (diff) *diff = _fm1;
158 return angle(phi.y() * _fm1, phi.x());
159 }
160
161 template<typename T>
162 AuxAngle<T> AuxLatitude<T>::Geocentric(const angle& phi, real* diff) const {
163 if (diff) *diff = _e2m1;
164 return angle(phi.y() * _e2m1, phi.x());
165 }
166
167 template<typename T>
168 AuxAngle<T> AuxLatitude<T>::Rectifying(const angle& phi, real* diff) const {
169 angle beta(Parametric(phi).normalized());
170 real sbeta = fabs(beta.y()), cbeta = fabs(beta.x());
171 real a = 1, b = _fm1, ka = _e2, kb = -_e12, ka1 = _e2m1, kb1 = _e12p1,
172 smu, cmu, mr;
173 if (_f < 0) {
174 swap(a, b); swap(ka, kb); swap(ka1, kb1); swap(sbeta, cbeta);
175 }
176 // now a,b = larger/smaller semiaxis
177 // beta now measured from larger semiaxis
178 // kb,ka = modulus-squared for distance from beta = 0,pi/2
179 // NB kb <= 0; 0 <= ka <= 1
180 // sa = b*E(beta,sqrt(kb)), sb = a*E(beta',sqrt(ka))
181 // 1 - ka * (1 - sb2) = 1 -ka + ka*sb2
182 real
183 sb2 = sbeta * sbeta,
184 cb2 = cbeta * cbeta,
185 db2 = 1 - kb * sb2,
186 da2 = ka1 + ka * sb2,
187 // DLMF Eq. 19.25.9
188 sa = b * sbeta * ( RF(cb2, db2, 1) - kb * sb2 * RD(cb2, db2, 1) / 3 ),
189 // DLMF Eq. 19.25.10 with complementary angles
190 sb = a * cbeta * ( ka1 * RF(sb2, da2, 1)
191 + ka * ka1 * cb2 * RD(sb2, 1, da2) / 3
192 + ka * sbeta / sqrt(da2) );
193 // sa + sb = 2*EllipticFunction::RG(a*a, b*b) = a*E(e) = b*E(i*e')
194 // mr = a*E(e)*(2/pi) = b*E(i*e')*(2/pi)
195 mr = (2 * (sa + sb)) / Math::pi<real>();
196 smu = sin(sa / mr);
197 cmu = sin(sb / mr);
198 if (_f < 0) { swap(smu, cmu); swap(a, b); }
199 // mu is normalized
200 angle mu(angle(smu, cmu).copyquadrant(phi));
201 if (diff) {
202 real cphi = phi.normalized().x(), tphi = phi.tan();
203 if (isfinite(tphi) || isnan(tphi)) {
204 cmu = mu.x(); cbeta = beta.x();
205 *diff = _fm1 * b/mr * Math::sq(cbeta / cmu) * (cbeta / cphi);
206 } else
207 *diff = _fm1 * mr/a;
208 }
209 return mu;
210 }
211
212 template<typename T>
213 AuxAngle<T> AuxLatitude<T>::Conformal(const angle& phi, real* diff) const {
214 real tphi = fabs(phi.tan()), tchi = tphi;
215 if ( !( !isfinite(tphi) || tphi == 0 || _f == 0 ) ) {
216 real scphi = sc(tphi),
217 sig = sinh(_e2 * atanhee(tphi) ),
218 scsig = sc(sig);
219#if AUXLATITUDE_UNOPT
220 tchi = tphi * scsig - sig * scphi;
221#else
222 if (_f <= 0) {
223 tchi = tphi * scsig - sig * scphi;
224 } else {
225 // The general expression for tchi is
226 // tphi * scsig - sig * scphi
227 // This involves cancellation if f > 0, so change to
228 // (tphi - sig) * (tphi + sig) / (tphi * scsig + sig * scphi)
229 // To control overflow, write as (sigtphi = sig / tphi)
230 // (tphi - sig) * (1 + sigtphi) / (scsig + sigtphi * scphi)
231 real sigtphi = sig / tphi, tphimsig;
232 if (sig < tphi / 2)
233 tphimsig = tphi - sig;
234 else {
235 // Still have possibly dangerous cancellation in tphi - sig.
236 //
237 // Write tphi - sig = (1 - e) * Dg(1, e)
238 // Dg(x, y) = (g(x) - g(y)) / (x - y)
239 // g(x) = sinh(x * atanh(sphi * x))
240 // Note sinh(atanh(sphi)) = tphi
241 // Turn the crank on divided differences, substitute
242 // sphi = tphi/sc(tphi)
243 // atanh(x) = asinh(x/sqrt(1-x^2))
244 real
245 em1 = _e2m1 / (1 + _e), // 1 - e
246 atanhs = asinh(tphi), // atanh(sphi)
247 scbeta = sc(_fm1 * tphi), // sec(beta)
248 scphibeta = sc(tphi) / scbeta, // sec(phi)/sec(beta)
249 atanhes = asinh(_e * tphi / scbeta), // atanh(e * sphi)
250 t1 = (atanhs - _e * atanhes)/2,
251 t2 = asinh(em1 * (tphi * scphibeta)) / em1,
252 Dg = cosh((atanhs + _e * atanhes)/2) * (sinh(t1) / t1)
253 * ((atanhs + atanhes)/2 + (1 + _e)/2 * t2);
254 tphimsig = em1 * Dg; // tphi - sig
255 }
256 tchi = tphimsig * (1 + sigtphi) / (scsig + sigtphi * scphi);
257 }
258#endif
259 }
260 angle chi(angle(tchi).copyquadrant(phi));
261 if (diff) {
262 if (isfinite(tphi) || isnan(tphi)) {
263 real cchi = chi.normalized().x(),
264 cphi = phi.normalized().x(),
265 cbeta = Parametric(phi).normalized().x();
266 *diff = _e2m1 * (cbeta / cchi) * (cbeta / cphi);
267 } else {
268#if AUXLATITUDE_UNOPT
269 real ss = _f > 0 ? sinh(_e * asinh(_e1)) : sinh(-_e * atan(_e));
270 *diff = _f > 0 ? 1/( sc(ss) + ss ) : sc(ss) - ss;
271#else
272 real ss = _f > 0 ? sinh(_e * atanh(_e)) : sinh(-_e * atan(_e));
273 *diff = sc(ss) - ss;
274#endif
275 }
276 }
277 return chi;
278 }
279
280 template<typename T>
281 AuxAngle<T> AuxLatitude<T>::Authalic(const angle& phi, real* diff) const {
282 real tphi = fabs(phi.tan());
283 angle xi(phi), phin(phi.normalized());
284 if ( !( !isfinite(tphi) || tphi == 0 || _f == 0 ) ) {
285 real qv = q(tphi),
286 Dqp = Dq(tphi),
287 Dqm = (_q + qv) / (1 + fabs(phin.y())); // Dq(-tphi)
288 xi = angle( copysign(qv, phi.y()), phin.x() * sqrt(Dqp * Dqm) );
289 }
290 if (diff) {
291 if (isfinite(tphi) || isnan(tphi)) {
292 real
293 cbeta = Parametric(phi).normalized().x(),
294 cxi = xi.normalized().x();
295 *diff =
296 (2/_q) * Math::sq(cbeta / cxi) * (cbeta / cxi) * (cbeta / phin.x());
297 } else
298 *diff = _e2m1 * sqrt(_q/2);
299 }
300 return xi;
301 }
302
303 template<typename T>
304 AuxAngle<T> AuxLatitude<T>::ToAuxiliary(int auxout, const angle& phi,
305 real* diff) const {
306 switch (auxout) {
307 case GEOGRAPHIC: if (diff) *diff = 1; return phi; break;
308 case PARAMETRIC: return Parametric(phi, diff); break;
309 case GEOCENTRIC: return Geocentric(phi, diff); break;
310 case RECTIFYING: return Rectifying(phi, diff); break;
311 case CONFORMAL : return Conformal (phi, diff); break;
312 case AUTHALIC : return Authalic (phi, diff); break;
313 default:
314 if (diff) *diff = numeric_limits<real>::quiet_NaN();
315 return angle::NaN();
316 break;
317 }
318 }
319
320 template<typename T>
321 AuxAngle<T> AuxLatitude<T>::FromAuxiliary(int auxin, const angle& zeta,
322 int* niter) const {
323 int n = 0; if (niter) *niter = n;
324 real tphi = _fm1;
325 switch (auxin) {
326 case GEOGRAPHIC: return zeta; break;
327 // case PARAMETRIC: break;
328 case PARAMETRIC: return angle(zeta.y() / _fm1, zeta.x()); break;
329 // case GEOCENTRIC: tphi *= _fm1 ; break;
330 case GEOCENTRIC: return angle(zeta.y() / _e2m1, zeta.x()); break;
331 case RECTIFYING: tphi *= sqrt(_fm1); break;
332 case CONFORMAL : tphi *= _fm1 ; break;
333 case AUTHALIC : tphi *= cbrt(_fm1); break;
334 default: return angle::NaN(); break;
335 }
336
337 // Drop through to solution by Newton's method
338#if AUXLATITUDE_UNOPT
339 real tzeta = fabs(zeta.tan());
340 if (!isfinite(tzeta)) return zeta;
341 tphi = tzeta / tphi;
342 real bmin = fmin(tphi, bmin_), bmax = fmax(tphi, bmax_);
343 for (int sign = 0, osign = 0, ntrip = 0; n < numit_;) {
344 ++n;
345 real diff;
346 angle zeta1(ToAuxiliary(auxin, angle(tphi), &diff));
347 real tzeta1 = zeta1.tan();
348 osign = sign;
349 if (tzeta1 == tzeta)
350 break;
351 else if (tzeta1 > tzeta) {
352 sign = 1;
353 bmax = tphi;
354 } else {
355 sign = -1;
356 bmin = tphi;
357 }
358 real dtphi = -(tzeta1 - tzeta) / diff;
359 tphi += dtphi;
360 if (!(fabs(dtphi) >= tol_ * tphi))
361 break;
362 // Skip the Newton guards if UNOPT
363 if (false && ((sign * osign < 0 && n - ntrip > 2) ||
364 tphi >= bmax || tphi <= bmin)) {
365 sign = 0; ntrip = n;
366 tphi = sqrt(bmin * bmax);
367 }
368 }
369#else
370 real tzeta = fabs(zeta.tan()), ltzeta = log2(tzeta);
371 if (!isfinite(ltzeta)) return zeta;
372 tphi = tzeta / tphi;
373 real ltphi = log2(tphi),
374 bmin = fmin(ltphi, bmin_), bmax = fmax(ltphi, bmax_);
375 for (int sign = 0, osign = 0, ntrip = 0; n < numit_;) {
376 ++n;
377 real diff;
378 angle zeta1(ToAuxiliary(auxin, angle(tphi), &diff));
379 real tzeta1 = zeta1.tan(), ltzeta1 = log2(tzeta1);
380 // Convert derivative from dtan(zeta)/dtan(phi) to
381 // dlog(tan(zeta))/dlog(tan(phi))
382 diff *= tphi/tzeta1;
383 osign = sign;
384 if (tzeta1 == tzeta)
385 break;
386 else if (tzeta1 > tzeta) {
387 sign = 1;
388 bmax = ltphi;
389 } else {
390 sign = -1;
391 bmin = ltphi;
392 }
393 real dltphi = -(ltzeta1 - ltzeta) / diff;
394 ltphi += dltphi;
395 tphi = exp2(ltphi);
396 if (!(fabs(dltphi) >= tol_)) {
397 ++n;
398 // Final Newton iteration without the logs
399 zeta1 = ToAuxiliary(auxin, angle(tphi), &diff);
400 tphi -= (zeta1.tan() - tzeta) / diff;
401 break;
402 }
403 if ((sign * osign < 0 && n - ntrip > 2) ||
404 ltphi >= bmax || ltphi <= bmin) {
405 sign = 0; ntrip = n;
406 ltphi = (bmin + bmax) / 2;
407 tphi = exp2(ltphi);
408 }
409 }
410#endif
411 if (niter) *niter = n;
412 return angle(tphi).copyquadrant(zeta);
413 }
414
415 template<typename T>
416 AuxAngle<T> AuxLatitude<T>::Convert(int auxin, int auxout, const angle& zeta,
417 bool series) const {
418 int k = ind(auxout, auxin);
419 if (k < 0) return angle::NaN();
420 if (auxin == auxout) return zeta;
421 if (series) {
422 if ( isnan(_c[Lmax * k]) ) fillcoeff(auxin, auxout, k);
423 angle zetan(zeta.normalized());
424 real d = Clenshaw(zetan.y(), zetan.x(), _c + Lmax * k, Lmax);
425 zetan += angle::radians(d);
426 return zetan;
427 } else {
428 if (auxin < 3 && auxout < 3)
429 return angle(zeta.y() * pow(_fm1, auxout - auxin), zeta.x());
430 else
431 return ToAuxiliary(auxout, FromAuxiliary(auxin, zeta));
432 }
433 }
434
435 template<typename T>
436 typename AuxLatitude<T>::real AuxLatitude<T>::RD(real x, real y, real z) {
437#if AUXLATITUDE_NATIVE_ELLIPTIC
438 // Carlson, eqs 2.28 - 2.34
439 static const real
440 tolRD = pow(real(0.2) * (numeric_limits<real>::epsilon() * real(0.01)),
441 1/real(8));
442 real
443 A0 = (x + y + 3*z)/5,
444 An = A0,
445 Q = fmax(fmax(fabs(A0-x), fabs(A0-y)), fabs(A0-z)) / tolRD,
446 x0 = x,
447 y0 = y,
448 z0 = z,
449 mul = 1,
450 s = 0;
451 while (Q >= mul * fabs(An)) {
452 // Max 7 trips
453 real lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0);
454 s += 1/(mul * sqrt(z0) * (z0 + lam));
455 An = (An + lam)/4;
456 x0 = (x0 + lam)/4;
457 y0 = (y0 + lam)/4;
458 z0 = (z0 + lam)/4;
459 mul *= 4;
460 }
461 real
462 X = (A0 - x) / (mul * An),
463 Y = (A0 - y) / (mul * An),
464 Z = -(X + Y) / 3,
465 E2 = X*Y - 6*Z*Z,
466 E3 = (3*X*Y - 8*Z*Z)*Z,
467 E4 = 3 * (X*Y - Z*Z) * Z*Z,
468 E5 = X*Y*Z*Z*Z;
469 // https://dlmf.nist.gov/19.36.E2
470 // Polynomial is
471 // (1 - 3*E2/14 + E3/6 + 9*E2^2/88 - 3*E4/22 - 9*E2*E3/52 + 3*E5/26
472 // - E2^3/16 + 3*E3^2/40 + 3*E2*E4/20 + 45*E2^2*E3/272
473 // - 9*(E3*E4+E2*E5)/68)
474 return ((471240 - 540540 * E2) * E5 +
475 (612612 * E2 - 540540 * E3 - 556920) * E4 +
476 E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
477 E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
478 (4084080 * mul * An * sqrt(An)) + 3 * s;
479#else
480 return real(EllipticFunction::RD(Math::real(x), Math::real(y),
481 Math::real(z)));
482#endif
483 }
484
485 template<typename T>
486 typename AuxLatitude<T>::real AuxLatitude<T>::RF(real x, real y, real z) {
487#if AUXLATITUDE_NATIVE_ELLIPTIC
488 // Carlson, eqs 2.2 - 2.7
489 static const real tolRF =
490 pow(3 * numeric_limits<real>::epsilon() * real(0.01), 1/real(8));
491 real
492 A0 = (x + y + z)/3,
493 An = A0,
494 Q = fmax(fmax(fabs(A0-x), fabs(A0-y)), fabs(A0-z)) / tolRF,
495 x0 = x,
496 y0 = y,
497 z0 = z,
498 mul = 1;
499 while (Q >= mul * fabs(An)) {
500 // Max 6 trips
501 real lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0);
502 An = (An + lam)/4;
503 x0 = (x0 + lam)/4;
504 y0 = (y0 + lam)/4;
505 z0 = (z0 + lam)/4;
506 mul *= 4;
507 }
508 real
509 X = (A0 - x) / (mul * An),
510 Y = (A0 - y) / (mul * An),
511 Z = - (X + Y),
512 E2 = X*Y - Z*Z,
513 E3 = X*Y*Z;
514 // https://dlmf.nist.gov/19.36.E1
515 // Polynomial is
516 // (1 - E2/10 + E3/14 + E2^2/24 - 3*E2*E3/44
517 // - 5*E2^3/208 + 3*E3^2/104 + E2^2*E3/16)
518 // convert to Horner form...
519 return (E3 * (6930 * E3 + E2 * (15015 * E2 - 16380) + 17160) +
520 E2 * ((10010 - 5775 * E2) * E2 - 24024) + 240240) /
521 (240240 * sqrt(An));
522#else
523 return real(EllipticFunction::RF(Math::real(x), Math::real(y),
524 Math::real(z)));
525#endif
526 }
527
528 template<typename T>
529 typename AuxLatitude<T>::real AuxLatitude<T>::atanhee(real tphi) const {
530#if AUXLATITUDE_UNOPT
531 real sphi = sn(tphi);
532 return _f == 0 ? sphi :
533 (_f < 0 ? atan( _e * sphi ) : atanh( _e * sphi )) / _e;
534#else
535 real s = _f <= 0 ? sn(tphi) : sn(_fm1 * tphi);
536 return _f == 0 ? s :
537 // atanh(e * sphi) = asinh(e' * sbeta)
538 (_f < 0 ? atan( _e * s ) : asinh( _e1 * s )) / _e;
539#endif
540 }
541
542 template<typename T>
543 typename AuxLatitude<T>::real AuxLatitude<T>::q(real tphi) const {
544#if AUXLATITUDE_UNOPT
545 real sphi = sn(tphi);
546 return atanhee(tphi) + sphi / (1 - _e2 * sphi*sphi);
547#else
548 real scbeta = sc(_fm1 * tphi);
549 return atanhee(tphi) + (tphi / scbeta) * (sc(tphi) / scbeta);
550#endif
551 }
552
553 template<typename T>
554 typename AuxLatitude<T>::real AuxLatitude<T>::Dq(real tphi) const {
555#if AUXLATITUDE_UNOPT
556 real sphi = sn(tphi), d = 1 - sphi;
557 if (tphi <= 0)
558 // This branch is not reached; this case is open-coded in Authalic.
559 return (_q - q(tphi)) / d;
560 else if (d == 0)
561 return 2 / Math::sq(_e2m1);
562 else
563 return (_f == 0 ? 1 / (1 - _e2 * sphi) :
564 (_f < 0 ? atan(_e * d / (1 - _e2 * sphi)) / (_e * d) :
565 atanh(_e * d / (1 - _e2 * sphi)) / (_e * d) )) +
566 (1 + _e2 * sphi) / ((1 - _e2 * sphi*sphi) * _e2m1);
567#else
568 real scphi = sc(tphi), sphi = sn(tphi),
569 // d = (1 - sphi) can underflow to zero for large tphi
570 d = tphi > 0 ? 1 / (scphi * scphi * (1 + sphi)) : 1 - sphi;
571 if (tphi <= 0)
572 // This branch is not reached; this case is open-coded in Authalic.
573 return (_q - q(tphi)) / d;
574 else if (d == 0)
575 return 2 / Math::sq(_e2m1);
576 else {
577 // General expression for Dq(1, sphi) is
578 // atanh(e * d / (1 - e2 * sphi)) / (e * d) +
579 // (1 + e2 * sphi) / ((1 - e2 * sphi * sphi) * e2m1);
580 // atanh( e * d / (1 - e2 * sphi))
581 // = atanh( e * d * scphi/(scphi - e2 * tphi))
582 // =
583 real scbeta = sc(_fm1 * tphi);
584 return (_f == 0 ? 1 :
585 (_f > 0 ? asinh(_e1 * d * scphi / scbeta) :
586 atan(_e * d / (1 - _e2 * sphi))) / (_e * d) ) +
587 (_f > 0 ?
588 ((scphi + _e2 * tphi) / (_e2m1 * scbeta)) * (scphi / scbeta) :
589 (1 + _e2 * sphi) / ((1 - _e2 * sphi*sphi) * _e2m1) );
590 }
591#endif
592 }
593
594 template<typename T>
595 void AuxLatitude<T>::fillcoeff(int auxin, int auxout, int k) const {
596#if GEOGRAPHICLIB_AUXLATITUDE_ORDER == 4
597 static const real coeffs[] = {
598 // C[phi,phi] skipped
599 // C[phi,beta]; even coeffs only
600 0, 1,
601 0, 1/real(2),
602 1/real(3),
603 1/real(4),
604 // C[phi,theta]; even coeffs only
605 -2, 2,
606 -4, 2,
607 8/real(3),
608 4,
609 // C[phi,mu]; even coeffs only
610 -27/real(32), 3/real(2),
611 -55/real(32), 21/real(16),
612 151/real(96),
613 1097/real(512),
614 // C[phi,chi]
615 116/real(45), -2, -2/real(3), 2,
616 -227/real(45), -8/real(5), 7/real(3),
617 -136/real(35), 56/real(15),
618 4279/real(630),
619 // C[phi,xi]
620 -2582/real(14175), -16/real(35), 4/real(45), 4/real(3),
621 -11966/real(14175), 152/real(945), 46/real(45),
622 3802/real(14175), 3044/real(2835),
623 6059/real(4725),
624 // C[beta,phi]; even coeffs only
625 0, -1,
626 0, 1/real(2),
627 -1/real(3),
628 1/real(4),
629 // C[beta,beta] skipped
630 // C[beta,theta]; even coeffs only
631 0, 1,
632 0, 1/real(2),
633 1/real(3),
634 1/real(4),
635 // C[beta,mu]; even coeffs only
636 -9/real(32), 1/real(2),
637 -37/real(96), 5/real(16),
638 29/real(96),
639 539/real(1536),
640 // C[beta,chi]
641 38/real(45), -1/real(3), -2/real(3), 1,
642 -7/real(9), -14/real(15), 5/real(6),
643 -34/real(21), 16/real(15),
644 2069/real(1260),
645 // C[beta,xi]
646 -1082/real(14175), -46/real(315), 4/real(45), 1/real(3),
647 -338/real(2025), 68/real(945), 17/real(90),
648 1102/real(14175), 461/real(2835),
649 3161/real(18900),
650 // C[theta,phi]; even coeffs only
651 2, -2,
652 -4, 2,
653 -8/real(3),
654 4,
655 // C[theta,beta]; even coeffs only
656 0, -1,
657 0, 1/real(2),
658 -1/real(3),
659 1/real(4),
660 // C[theta,theta] skipped
661 // C[theta,mu]; even coeffs only
662 -23/real(32), -1/real(2),
663 -5/real(96), 5/real(16),
664 1/real(32),
665 283/real(1536),
666 // C[theta,chi]
667 4/real(9), -2/real(3), -2/real(3), 0,
668 -23/real(45), -4/real(15), 1/real(3),
669 -24/real(35), 2/real(5),
670 83/real(126),
671 // C[theta,xi]
672 -2102/real(14175), -158/real(315), 4/real(45), -2/real(3),
673 934/real(14175), -16/real(945), 16/real(45),
674 922/real(14175), -232/real(2835),
675 719/real(4725),
676 // C[mu,phi]; even coeffs only
677 9/real(16), -3/real(2),
678 -15/real(32), 15/real(16),
679 -35/real(48),
680 315/real(512),
681 // C[mu,beta]; even coeffs only
682 3/real(16), -1/real(2),
683 1/real(32), -1/real(16),
684 -1/real(48),
685 -5/real(512),
686 // C[mu,theta]; even coeffs only
687 13/real(16), 1/real(2),
688 33/real(32), -1/real(16),
689 -5/real(16),
690 -261/real(512),
691 // C[mu,mu] skipped
692 // C[mu,chi]
693 41/real(180), 5/real(16), -2/real(3), 1/real(2),
694 557/real(1440), -3/real(5), 13/real(48),
695 -103/real(140), 61/real(240),
696 49561/real(161280),
697 // C[mu,xi]
698 -1609/real(28350), 121/real(1680), 4/real(45), -1/real(6),
699 16463/real(453600), 26/real(945), -29/real(720),
700 449/real(28350), -1003/real(45360),
701 -40457/real(2419200),
702 // C[chi,phi]
703 -82/real(45), 4/real(3), 2/real(3), -2,
704 -13/real(9), -16/real(15), 5/real(3),
705 34/real(21), -26/real(15),
706 1237/real(630),
707 // C[chi,beta]
708 -16/real(45), 0, 2/real(3), -1,
709 19/real(45), -2/real(5), 1/real(6),
710 16/real(105), -1/real(15),
711 17/real(1260),
712 // C[chi,theta]
713 -2/real(9), 2/real(3), 2/real(3), 0,
714 43/real(45), 4/real(15), -1/real(3),
715 2/real(105), -2/real(5),
716 -55/real(126),
717 // C[chi,mu]
718 1/real(360), -37/real(96), 2/real(3), -1/real(2),
719 437/real(1440), -1/real(15), -1/real(48),
720 37/real(840), -17/real(480),
721 -4397/real(161280),
722 // C[chi,chi] skipped
723 // C[chi,xi]
724 -2312/real(14175), -88/real(315), 34/real(45), -2/real(3),
725 6079/real(14175), -184/real(945), 1/real(45),
726 772/real(14175), -106/real(2835),
727 -167/real(9450),
728 // C[xi,phi]
729 538/real(4725), 88/real(315), -4/real(45), -4/real(3),
730 -2482/real(14175), 8/real(105), 34/real(45),
731 -898/real(14175), -1532/real(2835),
732 6007/real(14175),
733 // C[xi,beta]
734 34/real(675), 32/real(315), -4/real(45), -1/real(3),
735 74/real(2025), -4/real(315), -7/real(90),
736 2/real(14175), -83/real(2835),
737 -797/real(56700),
738 // C[xi,theta]
739 778/real(4725), 62/real(105), -4/real(45), 2/real(3),
740 12338/real(14175), -32/real(315), 4/real(45),
741 -1618/real(14175), -524/real(2835),
742 -5933/real(14175),
743 // C[xi,mu]
744 1297/real(18900), -817/real(10080), -4/real(45), 1/real(6),
745 -29609/real(453600), -2/real(35), 49/real(720),
746 -2917/real(56700), 4463/real(90720),
747 331799/real(7257600),
748 // C[xi,chi]
749 2458/real(4725), 46/real(315), -34/real(45), 2/real(3),
750 3413/real(14175), -256/real(315), 19/real(45),
751 -15958/real(14175), 248/real(567),
752 16049/real(28350),
753 // C[xi,xi] skipped
754 };
755 static const int ptrs[] = {
756 0, 0, 6, 12, 18, 28, 38, 44, 44, 50, 56, 66, 76, 82, 88, 88, 94, 104,
757 114, 120, 126, 132, 132, 142, 152, 162, 172, 182, 192, 192, 202, 212,
758 222, 232, 242, 252, 252,
759 };
760#elif GEOGRAPHICLIB_AUXLATITUDE_ORDER == 6
761 static const real coeffs[] = {
762 // C[phi,phi] skipped
763 // C[phi,beta]; even coeffs only
764 0, 0, 1,
765 0, 0, 1/real(2),
766 0, 1/real(3),
767 0, 1/real(4),
768 1/real(5),
769 1/real(6),
770 // C[phi,theta]; even coeffs only
771 2, -2, 2,
772 6, -4, 2,
773 -8, 8/real(3),
774 -16, 4,
775 32/real(5),
776 32/real(3),
777 // C[phi,mu]; even coeffs only
778 269/real(512), -27/real(32), 3/real(2),
779 6759/real(4096), -55/real(32), 21/real(16),
780 -417/real(128), 151/real(96),
781 -15543/real(2560), 1097/real(512),
782 8011/real(2560),
783 293393/real(61440),
784 // C[phi,chi]
785 -2854/real(675), 26/real(45), 116/real(45), -2, -2/real(3), 2,
786 2323/real(945), 2704/real(315), -227/real(45), -8/real(5), 7/real(3),
787 73814/real(2835), -1262/real(105), -136/real(35), 56/real(15),
788 -399572/real(14175), -332/real(35), 4279/real(630),
789 -144838/real(6237), 4174/real(315),
790 601676/real(22275),
791 // C[phi,xi]
792 28112932/real(212837625), 60136/real(467775), -2582/real(14175),
793 -16/real(35), 4/real(45), 4/real(3),
794 251310128/real(638512875), -21016/real(51975), -11966/real(14175),
795 152/real(945), 46/real(45),
796 -8797648/real(10945935), -94388/real(66825), 3802/real(14175),
797 3044/real(2835),
798 -1472637812/real(638512875), 41072/real(93555), 6059/real(4725),
799 455935736/real(638512875), 768272/real(467775),
800 4210684958LL/real(1915538625),
801 // C[beta,phi]; even coeffs only
802 0, 0, -1,
803 0, 0, 1/real(2),
804 0, -1/real(3),
805 0, 1/real(4),
806 -1/real(5),
807 1/real(6),
808 // C[beta,beta] skipped
809 // C[beta,theta]; even coeffs only
810 0, 0, 1,
811 0, 0, 1/real(2),
812 0, 1/real(3),
813 0, 1/real(4),
814 1/real(5),
815 1/real(6),
816 // C[beta,mu]; even coeffs only
817 205/real(1536), -9/real(32), 1/real(2),
818 1335/real(4096), -37/real(96), 5/real(16),
819 -75/real(128), 29/real(96),
820 -2391/real(2560), 539/real(1536),
821 3467/real(7680),
822 38081/real(61440),
823 // C[beta,chi]
824 -3118/real(4725), -1/real(3), 38/real(45), -1/real(3), -2/real(3), 1,
825 -247/real(270), 50/real(21), -7/real(9), -14/real(15), 5/real(6),
826 17564/real(2835), -5/real(3), -34/real(21), 16/real(15),
827 -49877/real(14175), -28/real(9), 2069/real(1260),
828 -28244/real(4455), 883/real(315),
829 797222/real(155925),
830 // C[beta,xi]
831 7947332/real(212837625), 11824/real(467775), -1082/real(14175),
832 -46/real(315), 4/real(45), 1/real(3),
833 39946703/real(638512875), -16672/real(155925), -338/real(2025),
834 68/real(945), 17/real(90),
835 -255454/real(1563705), -101069/real(467775), 1102/real(14175),
836 461/real(2835),
837 -189032762/real(638512875), 1786/real(18711), 3161/real(18900),
838 80274086/real(638512875), 88868/real(467775),
839 880980241/real(3831077250LL),
840 // C[theta,phi]; even coeffs only
841 -2, 2, -2,
842 6, -4, 2,
843 8, -8/real(3),
844 -16, 4,
845 -32/real(5),
846 32/real(3),
847 // C[theta,beta]; even coeffs only
848 0, 0, -1,
849 0, 0, 1/real(2),
850 0, -1/real(3),
851 0, 1/real(4),
852 -1/real(5),
853 1/real(6),
854 // C[theta,theta] skipped
855 // C[theta,mu]; even coeffs only
856 499/real(1536), -23/real(32), -1/real(2),
857 6565/real(12288), -5/real(96), 5/real(16),
858 -77/real(128), 1/real(32),
859 -4037/real(7680), 283/real(1536),
860 1301/real(7680),
861 17089/real(61440),
862 // C[theta,chi]
863 -3658/real(4725), 2/real(9), 4/real(9), -2/real(3), -2/real(3), 0,
864 61/real(135), 68/real(45), -23/real(45), -4/real(15), 1/real(3),
865 9446/real(2835), -46/real(35), -24/real(35), 2/real(5),
866 -34712/real(14175), -80/real(63), 83/real(126),
867 -2362/real(891), 52/real(45),
868 335882/real(155925),
869 // C[theta,xi]
870 216932/real(2627625), 109042/real(467775), -2102/real(14175),
871 -158/real(315), 4/real(45), -2/real(3),
872 117952358/real(638512875), -7256/real(155925), 934/real(14175),
873 -16/real(945), 16/real(45),
874 -7391576/real(54729675), -25286/real(66825), 922/real(14175),
875 -232/real(2835),
876 -67048172/real(638512875), 268/real(18711), 719/real(4725),
877 46774256/real(638512875), 14354/real(467775),
878 253129538/real(1915538625),
879 // C[mu,phi]; even coeffs only
880 -3/real(32), 9/real(16), -3/real(2),
881 135/real(2048), -15/real(32), 15/real(16),
882 105/real(256), -35/real(48),
883 -189/real(512), 315/real(512),
884 -693/real(1280),
885 1001/real(2048),
886 // C[mu,beta]; even coeffs only
887 -1/real(32), 3/real(16), -1/real(2),
888 -9/real(2048), 1/real(32), -1/real(16),
889 3/real(256), -1/real(48),
890 3/real(512), -5/real(512),
891 -7/real(1280),
892 -7/real(2048),
893 // C[mu,theta]; even coeffs only
894 -15/real(32), 13/real(16), 1/real(2),
895 -1673/real(2048), 33/real(32), -1/real(16),
896 349/real(256), -5/real(16),
897 963/real(512), -261/real(512),
898 -921/real(1280),
899 -6037/real(6144),
900 // C[mu,mu] skipped
901 // C[mu,chi]
902 7891/real(37800), -127/real(288), 41/real(180), 5/real(16), -2/real(3),
903 1/real(2),
904 -1983433/real(1935360), 281/real(630), 557/real(1440), -3/real(5),
905 13/real(48),
906 167603/real(181440), 15061/real(26880), -103/real(140), 61/real(240),
907 6601661/real(7257600), -179/real(168), 49561/real(161280),
908 -3418889/real(1995840), 34729/real(80640),
909 212378941/real(319334400),
910 // C[mu,xi]
911 12674323/real(851350500), -384229/real(14968800), -1609/real(28350),
912 121/real(1680), 4/real(45), -1/real(6),
913 -31621753811LL/real(1307674368000LL), -431/real(17325),
914 16463/real(453600), 26/real(945), -29/real(720),
915 -32844781/real(1751349600), 3746047/real(119750400), 449/real(28350),
916 -1003/real(45360),
917 10650637121LL/real(326918592000LL), 629/real(53460),
918 -40457/real(2419200),
919 205072597/real(20432412000LL), -1800439/real(119750400),
920 -59109051671LL/real(3923023104000LL),
921 // C[chi,phi]
922 4642/real(4725), 32/real(45), -82/real(45), 4/real(3), 2/real(3), -2,
923 -1522/real(945), 904/real(315), -13/real(9), -16/real(15), 5/real(3),
924 -12686/real(2835), 8/real(5), 34/real(21), -26/real(15),
925 -24832/real(14175), -12/real(5), 1237/real(630),
926 109598/real(31185), -734/real(315),
927 444337/real(155925),
928 // C[chi,beta]
929 -998/real(4725), 2/real(5), -16/real(45), 0, 2/real(3), -1,
930 -2/real(27), -22/real(105), 19/real(45), -2/real(5), 1/real(6),
931 116/real(567), -22/real(105), 16/real(105), -1/real(15),
932 2123/real(14175), -8/real(105), 17/real(1260),
933 128/real(4455), -1/real(105),
934 149/real(311850),
935 // C[chi,theta]
936 1042/real(4725), -14/real(45), -2/real(9), 2/real(3), 2/real(3), 0,
937 -712/real(945), -4/real(45), 43/real(45), 4/real(15), -1/real(3),
938 274/real(2835), 124/real(105), 2/real(105), -2/real(5),
939 21068/real(14175), -16/real(105), -55/real(126),
940 -9202/real(31185), -22/real(45),
941 -90263/real(155925),
942 // C[chi,mu]
943 -96199/real(604800), 81/real(512), 1/real(360), -37/real(96), 2/real(3),
944 -1/real(2),
945 1118711/real(3870720), -46/real(105), 437/real(1440), -1/real(15),
946 -1/real(48),
947 -5569/real(90720), 209/real(4480), 37/real(840), -17/real(480),
948 830251/real(7257600), 11/real(504), -4397/real(161280),
949 108847/real(3991680), -4583/real(161280),
950 -20648693/real(638668800),
951 // C[chi,chi] skipped
952 // C[chi,xi]
953 -55271278/real(212837625), 27128/real(93555), -2312/real(14175),
954 -88/real(315), 34/real(45), -2/real(3),
955 106691108/real(638512875), -65864/real(155925), 6079/real(14175),
956 -184/real(945), 1/real(45),
957 5921152/real(54729675), -14246/real(467775), 772/real(14175),
958 -106/real(2835),
959 75594328/real(638512875), -5312/real(467775), -167/real(9450),
960 2837636/real(638512875), -248/real(13365),
961 -34761247/real(1915538625),
962 // C[xi,phi]
963 -44732/real(2837835), 20824/real(467775), 538/real(4725), 88/real(315),
964 -4/real(45), -4/real(3),
965 -12467764/real(212837625), -37192/real(467775), -2482/real(14175),
966 8/real(105), 34/real(45),
967 100320856/real(1915538625), 54968/real(467775), -898/real(14175),
968 -1532/real(2835),
969 -5884124/real(70945875), 24496/real(467775), 6007/real(14175),
970 -839792/real(19348875), -23356/real(66825),
971 570284222/real(1915538625),
972 // C[xi,beta]
973 -70496/real(8513505), 2476/real(467775), 34/real(675), 32/real(315),
974 -4/real(45), -1/real(3),
975 53836/real(212837625), 3992/real(467775), 74/real(2025), -4/real(315),
976 -7/real(90),
977 -661844/real(1915538625), 7052/real(467775), 2/real(14175),
978 -83/real(2835),
979 1425778/real(212837625), 934/real(467775), -797/real(56700),
980 390088/real(212837625), -3673/real(467775),
981 -18623681/real(3831077250LL),
982 // C[xi,theta]
983 -4286228/real(42567525), -193082/real(467775), 778/real(4725),
984 62/real(105), -4/real(45), 2/real(3),
985 -61623938/real(70945875), 92696/real(467775), 12338/real(14175),
986 -32/real(315), 4/real(45),
987 427003576/real(1915538625), 612536/real(467775), -1618/real(14175),
988 -524/real(2835),
989 427770788/real(212837625), -8324/real(66825), -5933/real(14175),
990 -9153184/real(70945875), -320044/real(467775),
991 -1978771378/real(1915538625),
992 // C[xi,mu]
993 -9292991/real(302702400), 7764059/real(239500800), 1297/real(18900),
994 -817/real(10080), -4/real(45), 1/real(6),
995 36019108271LL/real(871782912000LL), 35474/real(467775),
996 -29609/real(453600), -2/real(35), 49/real(720),
997 3026004511LL/real(30648618000LL), -4306823/real(59875200),
998 -2917/real(56700), 4463/real(90720),
999 -368661577/real(4036032000LL), -102293/real(1871100),
1000 331799/real(7257600),
1001 -875457073/real(13621608000LL), 11744233/real(239500800),
1002 453002260127LL/real(7846046208000LL),
1003 // C[xi,chi]
1004 2706758/real(42567525), -55222/real(93555), 2458/real(4725),
1005 46/real(315), -34/real(45), 2/real(3),
1006 -340492279/real(212837625), 516944/real(467775), 3413/real(14175),
1007 -256/real(315), 19/real(45),
1008 4430783356LL/real(1915538625), 206834/real(467775), -15958/real(14175),
1009 248/real(567),
1010 62016436/real(70945875), -832976/real(467775), 16049/real(28350),
1011 -651151712/real(212837625), 15602/real(18711),
1012 2561772812LL/real(1915538625),
1013 // C[xi,xi] skipped
1014 };
1015 static const int ptrs[] = {
1016 0, 0, 12, 24, 36, 57, 78, 90, 90, 102, 114, 135, 156, 168, 180, 180, 192,
1017 213, 234, 246, 258, 270, 270, 291, 312, 333, 354, 375, 396, 396, 417,
1018 438, 459, 480, 501, 522, 522,
1019 };
1020#elif GEOGRAPHICLIB_AUXLATITUDE_ORDER == 8
1021 static const real coeffs[] = {
1022 // C[phi,phi] skipped
1023 // C[phi,beta]; even coeffs only
1024 0, 0, 0, 1,
1025 0, 0, 0, 1/real(2),
1026 0, 0, 1/real(3),
1027 0, 0, 1/real(4),
1028 0, 1/real(5),
1029 0, 1/real(6),
1030 1/real(7),
1031 1/real(8),
1032 // C[phi,theta]; even coeffs only
1033 -2, 2, -2, 2,
1034 -8, 6, -4, 2,
1035 16, -8, 8/real(3),
1036 40, -16, 4,
1037 -32, 32/real(5),
1038 -64, 32/real(3),
1039 128/real(7),
1040 32,
1041 // C[phi,mu]; even coeffs only
1042 -6607/real(24576), 269/real(512), -27/real(32), 3/real(2),
1043 -155113/real(122880), 6759/real(4096), -55/real(32), 21/real(16),
1044 87963/real(20480), -417/real(128), 151/real(96),
1045 2514467/real(245760), -15543/real(2560), 1097/real(512),
1046 -69119/real(6144), 8011/real(2560),
1047 -5962461/real(286720), 293393/real(61440),
1048 6459601/real(860160),
1049 332287993/real(27525120),
1050 // C[phi,chi]
1051 189416/real(99225), 16822/real(4725), -2854/real(675), 26/real(45),
1052 116/real(45), -2, -2/real(3), 2,
1053 141514/real(8505), -31256/real(1575), 2323/real(945), 2704/real(315),
1054 -227/real(45), -8/real(5), 7/real(3),
1055 -2363828/real(31185), 98738/real(14175), 73814/real(2835),
1056 -1262/real(105), -136/real(35), 56/real(15),
1057 14416399/real(935550), 11763988/real(155925), -399572/real(14175),
1058 -332/real(35), 4279/real(630),
1059 258316372/real(1216215), -2046082/real(31185), -144838/real(6237),
1060 4174/real(315),
1061 -2155215124LL/real(14189175), -115444544/real(2027025),
1062 601676/real(22275),
1063 -170079376/real(1216215), 38341552/real(675675),
1064 1383243703/real(11351340),
1065 // C[phi,xi]
1066 -1683291094/real(37574026875LL), 22947844/real(1915538625),
1067 28112932/real(212837625), 60136/real(467775), -2582/real(14175),
1068 -16/real(35), 4/real(45), 4/real(3),
1069 -14351220203LL/real(488462349375LL), 1228352/real(3007125),
1070 251310128/real(638512875), -21016/real(51975), -11966/real(14175),
1071 152/real(945), 46/real(45),
1072 505559334506LL/real(488462349375LL), 138128272/real(147349125),
1073 -8797648/real(10945935), -94388/real(66825), 3802/real(14175),
1074 3044/real(2835),
1075 973080708361LL/real(488462349375LL), -45079184/real(29469825),
1076 -1472637812/real(638512875), 41072/real(93555), 6059/real(4725),
1077 -1385645336626LL/real(488462349375LL), -550000184/real(147349125),
1078 455935736/real(638512875), 768272/real(467775),
1079 -2939205114427LL/real(488462349375LL), 443810768/real(383107725),
1080 4210684958LL/real(1915538625),
1081 101885255158LL/real(54273594375LL), 387227992/real(127702575),
1082 1392441148867LL/real(325641566250LL),
1083 // C[beta,phi]; even coeffs only
1084 0, 0, 0, -1,
1085 0, 0, 0, 1/real(2),
1086 0, 0, -1/real(3),
1087 0, 0, 1/real(4),
1088 0, -1/real(5),
1089 0, 1/real(6),
1090 -1/real(7),
1091 1/real(8),
1092 // C[beta,beta] skipped
1093 // C[beta,theta]; even coeffs only
1094 0, 0, 0, 1,
1095 0, 0, 0, 1/real(2),
1096 0, 0, 1/real(3),
1097 0, 0, 1/real(4),
1098 0, 1/real(5),
1099 0, 1/real(6),
1100 1/real(7),
1101 1/real(8),
1102 // C[beta,mu]; even coeffs only
1103 -4879/real(73728), 205/real(1536), -9/real(32), 1/real(2),
1104 -86171/real(368640), 1335/real(4096), -37/real(96), 5/real(16),
1105 2901/real(4096), -75/real(128), 29/real(96),
1106 1082857/real(737280), -2391/real(2560), 539/real(1536),
1107 -28223/real(18432), 3467/real(7680),
1108 -733437/real(286720), 38081/real(61440),
1109 459485/real(516096),
1110 109167851/real(82575360),
1111 // C[beta,chi]
1112 -25666/real(99225), 4769/real(4725), -3118/real(4725), -1/real(3),
1113 38/real(45), -1/real(3), -2/real(3), 1,
1114 193931/real(42525), -14404/real(4725), -247/real(270), 50/real(21),
1115 -7/real(9), -14/real(15), 5/real(6),
1116 -1709614/real(155925), -36521/real(14175), 17564/real(2835), -5/real(3),
1117 -34/real(21), 16/real(15),
1118 -637699/real(85050), 2454416/real(155925), -49877/real(14175),
1119 -28/real(9), 2069/real(1260),
1120 48124558/real(1216215), -20989/real(2835), -28244/real(4455),
1121 883/real(315),
1122 -16969807/real(1091475), -2471888/real(184275), 797222/real(155925),
1123 -1238578/real(42525), 2199332/real(225225),
1124 87600385/real(4540536),
1125 // C[beta,xi]
1126 -5946082372LL/real(488462349375LL), 9708931/real(1915538625),
1127 7947332/real(212837625), 11824/real(467775), -1082/real(14175),
1128 -46/real(315), 4/real(45), 1/real(3),
1129 190673521/real(69780335625LL), 164328266/real(1915538625),
1130 39946703/real(638512875), -16672/real(155925), -338/real(2025),
1131 68/real(945), 17/real(90),
1132 86402898356LL/real(488462349375LL), 236067184/real(1915538625),
1133 -255454/real(1563705), -101069/real(467775), 1102/real(14175),
1134 461/real(2835),
1135 110123070361LL/real(488462349375LL), -98401826/real(383107725),
1136 -189032762/real(638512875), 1786/real(18711), 3161/real(18900),
1137 -200020620676LL/real(488462349375LL), -802887278/real(1915538625),
1138 80274086/real(638512875), 88868/real(467775),
1139 -296107325077LL/real(488462349375LL), 66263486/real(383107725),
1140 880980241/real(3831077250LL),
1141 4433064236LL/real(18091198125LL), 37151038/real(127702575),
1142 495248998393LL/real(1302566265000LL),
1143 // C[theta,phi]; even coeffs only
1144 2, -2, 2, -2,
1145 -8, 6, -4, 2,
1146 -16, 8, -8/real(3),
1147 40, -16, 4,
1148 32, -32/real(5),
1149 -64, 32/real(3),
1150 -128/real(7),
1151 32,
1152 // C[theta,beta]; even coeffs only
1153 0, 0, 0, -1,
1154 0, 0, 0, 1/real(2),
1155 0, 0, -1/real(3),
1156 0, 0, 1/real(4),
1157 0, -1/real(5),
1158 0, 1/real(6),
1159 -1/real(7),
1160 1/real(8),
1161 // C[theta,theta] skipped
1162 // C[theta,mu]; even coeffs only
1163 -14321/real(73728), 499/real(1536), -23/real(32), -1/real(2),
1164 -201467/real(368640), 6565/real(12288), -5/real(96), 5/real(16),
1165 2939/real(4096), -77/real(128), 1/real(32),
1166 1155049/real(737280), -4037/real(7680), 283/real(1536),
1167 -19465/real(18432), 1301/real(7680),
1168 -442269/real(286720), 17089/real(61440),
1169 198115/real(516096),
1170 48689387/real(82575360),
1171 // C[theta,chi]
1172 64424/real(99225), 76/real(225), -3658/real(4725), 2/real(9), 4/real(9),
1173 -2/real(3), -2/real(3), 0,
1174 2146/real(1215), -2728/real(945), 61/real(135), 68/real(45),
1175 -23/real(45), -4/real(15), 1/real(3),
1176 -95948/real(10395), 428/real(945), 9446/real(2835), -46/real(35),
1177 -24/real(35), 2/real(5),
1178 29741/real(85050), 4472/real(525), -34712/real(14175), -80/real(63),
1179 83/real(126),
1180 280108/real(13365), -17432/real(3465), -2362/real(891), 52/real(45),
1181 -48965632/real(4729725), -548752/real(96525), 335882/real(155925),
1182 -197456/real(15795), 51368/real(12285),
1183 1461335/real(174636),
1184 // C[theta,xi]
1185 -230886326/real(6343666875LL), -189115382/real(1915538625),
1186 216932/real(2627625), 109042/real(467775), -2102/real(14175),
1187 -158/real(315), 4/real(45), -2/real(3),
1188 -11696145869LL/real(69780335625LL), 288456008/real(1915538625),
1189 117952358/real(638512875), -7256/real(155925), 934/real(14175),
1190 -16/real(945), 16/real(45),
1191 91546732346LL/real(488462349375LL), 478700902/real(1915538625),
1192 -7391576/real(54729675), -25286/real(66825), 922/real(14175),
1193 -232/real(2835),
1194 218929662961LL/real(488462349375LL), -67330724/real(383107725),
1195 -67048172/real(638512875), 268/real(18711), 719/real(4725),
1196 -129039188386LL/real(488462349375LL), -117954842/real(273648375),
1197 46774256/real(638512875), 14354/real(467775),
1198 -178084928947LL/real(488462349375LL), 2114368/real(34827975),
1199 253129538/real(1915538625),
1200 6489189398LL/real(54273594375LL), 13805944/real(127702575),
1201 59983985827LL/real(325641566250LL),
1202 // C[mu,phi]; even coeffs only
1203 57/real(2048), -3/real(32), 9/real(16), -3/real(2),
1204 -105/real(4096), 135/real(2048), -15/real(32), 15/real(16),
1205 -105/real(2048), 105/real(256), -35/real(48),
1206 693/real(16384), -189/real(512), 315/real(512),
1207 693/real(2048), -693/real(1280),
1208 -1287/real(4096), 1001/real(2048),
1209 -6435/real(14336),
1210 109395/real(262144),
1211 // C[mu,beta]; even coeffs only
1212 19/real(2048), -1/real(32), 3/real(16), -1/real(2),
1213 7/real(4096), -9/real(2048), 1/real(32), -1/real(16),
1214 -3/real(2048), 3/real(256), -1/real(48),
1215 -11/real(16384), 3/real(512), -5/real(512),
1216 7/real(2048), -7/real(1280),
1217 9/real(4096), -7/real(2048),
1218 -33/real(14336),
1219 -429/real(262144),
1220 // C[mu,theta]; even coeffs only
1221 509/real(2048), -15/real(32), 13/real(16), 1/real(2),
1222 2599/real(4096), -1673/real(2048), 33/real(32), -1/real(16),
1223 -2989/real(2048), 349/real(256), -5/real(16),
1224 -43531/real(16384), 963/real(512), -261/real(512),
1225 5545/real(2048), -921/real(1280),
1226 16617/real(4096), -6037/real(6144),
1227 -19279/real(14336),
1228 -490925/real(262144),
1229 // C[mu,mu] skipped
1230 // C[mu,chi]
1231 -18975107/real(50803200), 72161/real(387072), 7891/real(37800),
1232 -127/real(288), 41/real(180), 5/real(16), -2/real(3), 1/real(2),
1233 148003883/real(174182400), 13769/real(28800), -1983433/real(1935360),
1234 281/real(630), 557/real(1440), -3/real(5), 13/real(48),
1235 79682431/real(79833600), -67102379/real(29030400), 167603/real(181440),
1236 15061/real(26880), -103/real(140), 61/real(240),
1237 -40176129013LL/real(7664025600LL), 97445/real(49896),
1238 6601661/real(7257600), -179/real(168), 49561/real(161280),
1239 2605413599LL/real(622702080), 14644087/real(9123840),
1240 -3418889/real(1995840), 34729/real(80640),
1241 175214326799LL/real(58118860800LL), -30705481/real(10378368),
1242 212378941/real(319334400),
1243 -16759934899LL/real(3113510400LL), 1522256789/real(1383782400),
1244 1424729850961LL/real(743921418240LL),
1245 // C[mu,xi]
1246 -375027460897LL/real(125046361440000LL),
1247 7183403063LL/real(560431872000LL), 12674323/real(851350500),
1248 -384229/real(14968800), -1609/real(28350), 121/real(1680), 4/real(45),
1249 -1/real(6),
1250 30410873385097LL/real(2000741783040000LL),
1251 1117820213/real(122594472000LL), -31621753811LL/real(1307674368000LL),
1252 -431/real(17325), 16463/real(453600), 26/real(945), -29/real(720),
1253 151567502183LL/real(17863765920000LL),
1254 -116359346641LL/real(3923023104000LL), -32844781/real(1751349600),
1255 3746047/real(119750400), 449/real(28350), -1003/real(45360),
1256 -317251099510901LL/real(8002967132160000LL), -13060303/real(766215450),
1257 10650637121LL/real(326918592000LL), 629/real(53460),
1258 -40457/real(2419200),
1259 -2105440822861LL/real(125046361440000LL),
1260 146875240637LL/real(3923023104000LL), 205072597/real(20432412000LL),
1261 -1800439/real(119750400),
1262 91496147778023LL/real(2000741783040000LL), 228253559/real(24518894400LL),
1263 -59109051671LL/real(3923023104000LL),
1264 126430355893LL/real(13894040160000LL),
1265 -4255034947LL/real(261534873600LL),
1266 -791820407649841LL/real(42682491371520000LL),
1267 // C[chi,phi]
1268 1514/real(1323), -8384/real(4725), 4642/real(4725), 32/real(45),
1269 -82/real(45), 4/real(3), 2/real(3), -2,
1270 142607/real(42525), -2288/real(1575), -1522/real(945), 904/real(315),
1271 -13/real(9), -16/real(15), 5/real(3),
1272 120202/real(51975), 44644/real(14175), -12686/real(2835), 8/real(5),
1273 34/real(21), -26/real(15),
1274 -1097407/real(187110), 1077964/real(155925), -24832/real(14175),
1275 -12/real(5), 1237/real(630),
1276 -12870194/real(1216215), 1040/real(567), 109598/real(31185),
1277 -734/real(315),
1278 -126463/real(72765), -941912/real(184275), 444337/real(155925),
1279 3463678/real(467775), -2405834/real(675675),
1280 256663081/real(56756700),
1281 // C[chi,beta]
1282 1384/real(11025), -34/real(4725), -998/real(4725), 2/real(5),
1283 -16/real(45), 0, 2/real(3), -1,
1284 -12616/real(42525), 1268/real(4725), -2/real(27), -22/real(105),
1285 19/real(45), -2/real(5), 1/real(6),
1286 1724/real(51975), -1858/real(14175), 116/real(567), -22/real(105),
1287 16/real(105), -1/real(15),
1288 115249/real(935550), -26836/real(155925), 2123/real(14175), -8/real(105),
1289 17/real(1260),
1290 140836/real(1216215), -424/real(6237), 128/real(4455), -1/real(105),
1291 210152/real(4729725), -31232/real(2027025), 149/real(311850),
1292 30208/real(6081075), -499/real(225225),
1293 -68251/real(113513400),
1294 // C[chi,theta]
1295 -1738/real(11025), 18/real(175), 1042/real(4725), -14/real(45),
1296 -2/real(9), 2/real(3), 2/real(3), 0,
1297 23159/real(42525), 332/real(945), -712/real(945), -4/real(45),
1298 43/real(45), 4/real(15), -1/real(3),
1299 13102/real(31185), -1352/real(945), 274/real(2835), 124/real(105),
1300 2/real(105), -2/real(5),
1301 -2414843/real(935550), 1528/real(4725), 21068/real(14175), -16/real(105),
1302 -55/real(126),
1303 60334/real(93555), 20704/real(10395), -9202/real(31185), -22/real(45),
1304 40458083/real(14189175), -299444/real(675675), -90263/real(155925),
1305 -3818498/real(6081075), -8962/real(12285),
1306 -4259027/real(4365900),
1307 // C[chi,mu]
1308 -7944359/real(67737600), 5406467/real(38707200), -96199/real(604800),
1309 81/real(512), 1/real(360), -37/real(96), 2/real(3), -1/real(2),
1310 -24749483/real(348364800), -51841/real(1209600), 1118711/real(3870720),
1311 -46/real(105), 437/real(1440), -1/real(15), -1/real(48),
1312 6457463/real(17740800), -9261899/real(58060800), -5569/real(90720),
1313 209/real(4480), 37/real(840), -17/real(480),
1314 -324154477/real(7664025600LL), -466511/real(2494800),
1315 830251/real(7257600), 11/real(504), -4397/real(161280),
1316 -22894433/real(124540416), 8005831/real(63866880), 108847/real(3991680),
1317 -4583/real(161280),
1318 2204645983LL/real(12915302400LL), 16363163/real(518918400),
1319 -20648693/real(638668800),
1320 497323811/real(12454041600LL), -219941297/real(5535129600LL),
1321 -191773887257LL/real(3719607091200LL),
1322 // C[chi,chi] skipped
1323 // C[chi,xi]
1324 -17451293242LL/real(488462349375LL), 308365186/real(1915538625),
1325 -55271278/real(212837625), 27128/real(93555), -2312/real(14175),
1326 -88/real(315), 34/real(45), -2/real(3),
1327 -101520127208LL/real(488462349375LL), 149984636/real(1915538625),
1328 106691108/real(638512875), -65864/real(155925), 6079/real(14175),
1329 -184/real(945), 1/real(45),
1330 10010741462LL/real(37574026875LL), -99534832/real(383107725),
1331 5921152/real(54729675), -14246/real(467775), 772/real(14175),
1332 -106/real(2835),
1333 1615002539/real(75148053750LL), -35573728/real(273648375),
1334 75594328/real(638512875), -5312/real(467775), -167/real(9450),
1335 -3358119706LL/real(488462349375LL), 130601488/real(1915538625),
1336 2837636/real(638512875), -248/real(13365),
1337 46771947158LL/real(488462349375LL), -3196/real(3553875),
1338 -34761247/real(1915538625),
1339 -18696014/real(18091198125LL), -2530364/real(127702575),
1340 -14744861191LL/real(651283132500LL),
1341 // C[xi,phi]
1342 -88002076/real(13956067125LL), -86728/real(16372125),
1343 -44732/real(2837835), 20824/real(467775), 538/real(4725), 88/real(315),
1344 -4/real(45), -4/real(3),
1345 -2641983469LL/real(488462349375LL), -895712/real(147349125),
1346 -12467764/real(212837625), -37192/real(467775), -2482/real(14175),
1347 8/real(105), 34/real(45),
1348 8457703444LL/real(488462349375LL), 240616/real(4209975),
1349 100320856/real(1915538625), 54968/real(467775), -898/real(14175),
1350 -1532/real(2835),
1351 -4910552477LL/real(97692469875LL), -4832848/real(147349125),
1352 -5884124/real(70945875), 24496/real(467775), 6007/real(14175),
1353 9393713176LL/real(488462349375LL), 816824/real(13395375),
1354 -839792/real(19348875), -23356/real(66825),
1355 -4532926649LL/real(97692469875LL), 1980656/real(54729675),
1356 570284222/real(1915538625),
1357 -14848113968LL/real(488462349375LL), -496894276/real(1915538625),
1358 224557742191LL/real(976924698750LL),
1359 // C[xi,beta]
1360 29232878/real(97692469875LL), -18484/real(4343625), -70496/real(8513505),
1361 2476/real(467775), 34/real(675), 32/real(315), -4/real(45), -1/real(3),
1362 -324943819/real(488462349375LL), -4160804/real(1915538625),
1363 53836/real(212837625), 3992/real(467775), 74/real(2025), -4/real(315),
1364 -7/real(90),
1365 -168643106/real(488462349375LL), 237052/real(383107725),
1366 -661844/real(1915538625), 7052/real(467775), 2/real(14175),
1367 -83/real(2835),
1368 113042383/real(97692469875LL), -2915326/real(1915538625),
1369 1425778/real(212837625), 934/real(467775), -797/real(56700),
1370 -558526274/real(488462349375LL), 6064888/real(1915538625),
1371 390088/real(212837625), -3673/real(467775),
1372 155665021/real(97692469875LL), 41288/real(29469825),
1373 -18623681/real(3831077250LL),
1374 504234982/real(488462349375LL), -6205669/real(1915538625),
1375 -8913001661LL/real(3907698795000LL),
1376 // C[xi,theta]
1377 182466964/real(8881133625LL), 53702182/real(212837625),
1378 -4286228/real(42567525), -193082/real(467775), 778/real(4725),
1379 62/real(105), -4/real(45), 2/real(3),
1380 367082779691LL/real(488462349375LL), -32500616/real(273648375),
1381 -61623938/real(70945875), 92696/real(467775), 12338/real(14175),
1382 -32/real(315), 4/real(45),
1383 -42668482796LL/real(488462349375LL), -663111728/real(383107725),
1384 427003576/real(1915538625), 612536/real(467775), -1618/real(14175),
1385 -524/real(2835),
1386 -327791986997LL/real(97692469875LL), 421877252/real(1915538625),
1387 427770788/real(212837625), -8324/real(66825), -5933/real(14175),
1388 74612072536LL/real(488462349375LL), 6024982024LL/real(1915538625),
1389 -9153184/real(70945875), -320044/real(467775),
1390 489898512247LL/real(97692469875LL), -46140784/real(383107725),
1391 -1978771378/real(1915538625),
1392 -42056042768LL/real(488462349375LL), -2926201612LL/real(1915538625),
1393 -2209250801969LL/real(976924698750LL),
1394 // C[xi,mu]
1395 39534358147LL/real(2858202547200LL),
1396 -25359310709LL/real(1743565824000LL), -9292991/real(302702400),
1397 7764059/real(239500800), 1297/real(18900), -817/real(10080), -4/real(45),
1398 1/real(6),
1399 -13216941177599LL/real(571640509440000LL),
1400 -14814966289LL/real(245188944000LL), 36019108271LL/real(871782912000LL),
1401 35474/real(467775), -29609/real(453600), -2/real(35), 49/real(720),
1402 -27782109847927LL/real(250092722880000LL),
1403 99871724539LL/real(1569209241600LL), 3026004511LL/real(30648618000LL),
1404 -4306823/real(59875200), -2917/real(56700), 4463/real(90720),
1405 168979300892599LL/real(1600593426432000LL),
1406 2123926699/real(15324309000LL), -368661577/real(4036032000LL),
1407 -102293/real(1871100), 331799/real(7257600),
1408 1959350112697LL/real(9618950880000LL),
1409 -493031379277LL/real(3923023104000LL), -875457073/real(13621608000LL),
1410 11744233/real(239500800),
1411 -145659994071373LL/real(800296713216000LL),
1412 -793693009/real(9807557760LL), 453002260127LL/real(7846046208000LL),
1413 -53583096419057LL/real(500185445760000LL),
1414 103558761539LL/real(1426553856000LL),
1415 12272105438887727LL/real(128047474114560000LL),
1416 // C[xi,chi]
1417 -64724382148LL/real(97692469875LL), 16676974/real(30405375),
1418 2706758/real(42567525), -55222/real(93555), 2458/real(4725),
1419 46/real(315), -34/real(45), 2/real(3),
1420 85904355287LL/real(37574026875LL), 158999572/real(1915538625),
1421 -340492279/real(212837625), 516944/real(467775), 3413/real(14175),
1422 -256/real(315), 19/real(45),
1423 2986003168LL/real(37574026875LL), -7597644214LL/real(1915538625),
1424 4430783356LL/real(1915538625), 206834/real(467775), -15958/real(14175),
1425 248/real(567),
1426 -375566203/real(39037950), 851209552/real(174139875),
1427 62016436/real(70945875), -832976/real(467775), 16049/real(28350),
1428 5106181018156LL/real(488462349375LL), 3475643362LL/real(1915538625),
1429 -651151712/real(212837625), 15602/real(18711),
1430 34581190223LL/real(8881133625LL), -10656173804LL/real(1915538625),
1431 2561772812LL/real(1915538625),
1432 -5150169424688LL/real(488462349375LL), 873037408/real(383107725),
1433 7939103697617LL/real(1953849397500LL),
1434 // C[xi,xi] skipped
1435 };
1436 static const int ptrs[] = {
1437 0, 0, 20, 40, 60, 96, 132, 152, 152, 172, 192, 228, 264, 284, 304, 304,
1438 324, 360, 396, 416, 436, 456, 456, 492, 528, 564, 600, 636, 672, 672,
1439 708, 744, 780, 816, 852, 888, 888,
1440 };
1441#else
1442#error "Unsupported value for GEOGRAPHICLIB_AUXLATITUDE_ORDER"
1443#endif
1444
1445 static_assert(sizeof(ptrs) / sizeof(int) == AUXNUMBER*AUXNUMBER+1,
1446 "Mismatch in size of ptrs array");
1447 static_assert(sizeof(coeffs) / sizeof(real) ==
1448 (RECTIFYING+1)*RECTIFYING *
1449 (Lmax * (Lmax + 3) - 2*(Lmax/2))/4 // Even only arrays
1450 + (AUXNUMBER*(AUXNUMBER-1) - (RECTIFYING+1)*RECTIFYING) *
1451 (Lmax * (Lmax + 1))/2,
1452 "Mismatch in size of coeffs array");
1453
1454 if (k < 0) return; // auxout or auxin out of range
1455 if (auxout == auxin)
1456 fill(_c + Lmax * k, _c + Lmax * (k + 1), 0);
1457 else {
1458 int o = ptrs[k];
1459 real d = _n;
1460 if (auxin <= RECTIFYING && auxout <= RECTIFYING) {
1461 for (int l = 0; l < Lmax; ++l) {
1462 int m = (Lmax - l - 1) / 2; // order of polynomial in n^2
1463 _c[Lmax * k + l] = d * Math::polyval(m, coeffs + o, _n2);
1464 o += m + 1;
1465 d *= _n;
1466 }
1467 } else {
1468 for (int l = 0; l < Lmax; ++l) {
1469 int m = (Lmax - l - 1); // order of polynomial in n
1470 _c[Lmax * k + l] = d * Math::polyval(m, coeffs + o, _n);
1471 o += m + 1;
1472 d *= _n;
1473 }
1474 }
1475 // assert (o == ptrs[AUXNUMBER * auxout + auxin + 1])
1476 }
1477 }
1478
1479 template<typename T>
1480 typename AuxLatitude<T>::real
1481 AuxLatitude<T>::Clenshaw(real szeta, real czeta,
1482 const real c[], int K, bool alt) {
1483 // Evaluate
1484 // y = sum(c[k] * sin( (2*k+2) * zeta), i, 0, K-1)
1485 // Approx operation count = (K + 5) mult and (2 * K + 2) add
1486 // If alt, use the Reinsch optimizations for small szeta and small czeta.
1487 int k = K;
1488 real u0 = 0; // accumulator for sum
1489 if (alt && 2*fabs(szeta) < 1) {
1490 real X = -4 * szeta * szeta,
1491 d0 = 0; // other accumulator for sum
1492 for (; k > 0;) {
1493 d0 = X * u0 + d0 + c[--k];
1494 u0 = d0 + u0;
1495 }
1496 } else if (alt && 2*fabs(czeta) < 1) {
1497 real X = 4 * czeta * czeta,
1498 d0 = 0; // other accumulator for sum
1499 for (; k > 0;) {
1500 d0 = X * u0 - d0 + c[--k];
1501 u0 = d0 - u0;
1502 }
1503 } else {
1504 real X = 2 * (czeta - szeta) * (czeta + szeta), // 2 * cos(2*zeta)
1505 u1 = 0; // other accumulator for sum
1506 for (; k > 0;) {
1507 real t = X * u0 - u1 + c[--k];
1508 u1 = u0; u0 = t;
1509 }
1510 }
1511 return 2 * szeta * czeta * u0; // sin(2*zeta) * u0
1512 }
1513
1514 template class AuxAngle<Math::real>;
1515 template class AuxLatitude<Math::real>;
1516#if GEOGRAPHICLIB_PRECISION != 2
1517 template class AuxAngle<double>;
1518 template class AuxLatitude<double>;
1519#endif
1520
1521} // namespace GeographicLib
AuxLatitude.hpp
Header for the GeographicLib::AuxLatitude and GeographicLib::AuxAngle classes.
EllipticFunction.hpp
Header for GeographicLib::EllipticFunction class.
real
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
GeographicLib::AuxAngle
An accurate representation of angles.
Definition: AuxLatitude.hpp:55
GeographicLib::AuxAngle::operator+=
AuxAngle & operator+=(const AuxAngle &p)
Definition: AuxLatitude.cpp:79
GeographicLib::AuxAngle::copyquadrant
AuxAngle copyquadrant(const AuxAngle &p) const
Definition: AuxLatitude.cpp:74
GeographicLib::AuxAngle::NaN
static AuxAngle NaN()
Definition: AuxLatitude.cpp:52
GeographicLib::AuxAngle::normalized
AuxAngle normalized() const
Definition: AuxLatitude.cpp:58
GeographicLib::AuxLatitude
Conversions between auxiliary latitudes.
Definition: AuxLatitude.hpp:271
GeographicLib::AuxLatitude::FromAuxiliary
angle FromAuxiliary(int auxin, const angle &zeta, int *niter=nullptr) const
Definition: AuxLatitude.cpp:321
GeographicLib::AuxLatitude::ToAuxiliary
angle ToAuxiliary(int auxout, const angle &phi, real *diff=nullptr) const
Definition: AuxLatitude.cpp:304
GeographicLib::AuxLatitude::Convert
angle Convert(int auxin, int auxout, const angle &zeta, bool series=false) const
Definition: AuxLatitude.cpp:416
GeographicLib::EllipticFunction::RD
static real RD(real x, real y, real z)
Definition: EllipticFunction.cpp:182
GeographicLib::EllipticFunction::RF
static real RF(real x, real y, real z)
Definition: EllipticFunction.cpp:29
GeographicLib::Geocentric
Geocentric coordinates
Definition: Geocentric.hpp:67
GeographicLib::GeographicErr
Exception handling for GeographicLib.
Definition: Constants.hpp:316
GeographicLib::Math::sq
static T sq(T x)
Definition: Math.hpp:212
GeographicLib::Math::polyval
static T polyval(int N, const T p[], T x)
Definition: Math.hpp:271
GeographicLib
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
std::swap
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
Definition: NearestNeighbor.hpp:826