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