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