GeographicLib: EllipticFunction.cpp Source File (original) (raw)
1
2
3
4
5
6
7
8
9
11
13
14 using namespace std;
15
16
17
18
19
20
21
22
23
25
26 static const real tolRF =
27 pow(3 * numeric_limits::epsilon() * real(0.01), 1/real(8));
28 real
29 A0 = (x + y + z)/3,
30 An = A0,
31 Q = fmax(fmax(fabs(A0-x), fabs(A0-y)), fabs(A0-z)) / tolRF,
32 x0 = x,
33 y0 = y,
34 z0 = z,
35 mul = 1;
36 while (Q >= mul * fabs(An)) {
37
38 real lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0);
39 An = (An + lam)/4;
40 x0 = (x0 + lam)/4;
41 y0 = (y0 + lam)/4;
42 z0 = (z0 + lam)/4;
43 mul *= 4;
44 }
45 real
46 X = (A0 - x) / (mul * An),
47 Y = (A0 - y) / (mul * An),
48 Z = - (X + Y),
49 E2 = X*Y - Z*Z,
50 E3 = X*Y*Z;
51
52
53
54
55
56 return (E3 * (6930 * E3 + E2 * (15015 * E2 - 16380) + 17160) +
57 E2 * ((10010 - 5775 * E2) * E2 - 24024) + 240240) /
58 (240240 * sqrt(An));
59 }
60
62
63 static const real tolRG0 =
64 real(2.7) * sqrt((numeric_limits::epsilon() * real(0.01)));
65 real xn = sqrt(x), yn = sqrt(y);
66 if (xn < yn) swap(xn, yn);
67 while (fabs(xn-yn) > tolRG0 * xn) {
68
69 real t = (xn + yn) /2;
70 yn = sqrt(xn * yn);
71 xn = t;
72 }
73 return Math::pi() / (xn + yn);
74 }
75
77
78 return ( !(x >= y) ?
79
80 atan(sqrt((y - x) / x)) / sqrt(y - x) :
81 ( x == y ? 1 / sqrt(y) :
82 asinh( y > 0 ?
83
84
85 sqrt((x - y) / y) :
86
87
88 sqrt(-x / y) ) / sqrt(x - y) ) );
89 }
90
92 return (x == 0 ? RG(y, z) :
93 (y == 0 ? RG(z, x) :
94 (z == 0 ? RG(x, y) :
95
96 (z * RF(x, y, z) - (x-z) * (y-z) * RD(x, y, z) / 3
97 + sqrt(x * y / z)) / 2 )));
98 }
99
101
102 static const real tolRG0 =
103 real(2.7) * sqrt((numeric_limits::epsilon() * real(0.01)));
104 real
105 x0 = sqrt(fmax(x, y)),
106 y0 = sqrt(fmin(x, y)),
107 xn = x0,
108 yn = y0,
109 s = 0,
110 mul = real(0.25);
111 while (fabs(xn-yn) > tolRG0 * xn) {
112
113 real t = (xn + yn) /2;
114 yn = sqrt(xn * yn);
115 xn = t;
116 mul *= 2;
117 t = xn - yn;
118 s += mul * t * t;
119 }
120 return (Math::sq( (x0 + y0)/2 ) - s) * Math::pi() / (2 * (xn + yn));
121 }
122
124
125 static const real
126 tolRD = pow(real(0.2) * (numeric_limits::epsilon() * real(0.01)),
127 1/real(8));
128 real
129 A0 = (x + y + z + 2*p)/5,
130 An = A0,
131 delta = (p-x) * (p-y) * (p-z),
132 Q = fmax(fmax(fabs(A0-x), fabs(A0-y)),
133 fmax(fabs(A0-z), fabs(A0-p))) / tolRD,
134 x0 = x,
135 y0 = y,
136 z0 = z,
137 p0 = p,
138 mul = 1,
139 mul3 = 1,
140 s = 0;
141 while (Q >= mul * fabs(An)) {
142
143 real
144 lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0),
145 d0 = (sqrt(p0)+sqrt(x0)) * (sqrt(p0)+sqrt(y0)) * (sqrt(p0)+sqrt(z0)),
146 e0 = delta/(mul3 * Math::sq(d0));
147 s += RC(1, 1 + e0)/(mul * d0);
148 An = (An + lam)/4;
149 x0 = (x0 + lam)/4;
150 y0 = (y0 + lam)/4;
151 z0 = (z0 + lam)/4;
152 p0 = (p0 + lam)/4;
153 mul *= 4;
154 mul3 *= 64;
155 }
156 real
157 X = (A0 - x) / (mul * An),
158 Y = (A0 - y) / (mul * An),
159 Z = (A0 - z) / (mul * An),
160 P = -(X + Y + Z) / 2,
161 E2 = X*Y + X*Z + Y*Z - 3*P*P,
162 E3 = X*Y*Z + 2*P * (E2 + 2*P*P),
163 E4 = (2*X*Y*Z + P * (E2 + 3*P*P)) * P,
164 E5 = X*Y*Z*P*P;
165
166
167
168
169
170 return ((471240 - 540540 * E2) * E5 +
171 (612612 * E2 - 540540 * E3 - 556920) * E4 +
172 E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
173 E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
174 (4084080 * mul * An * sqrt(An)) + 6 * s;
175 }
176
178
179 static const real
180 tolRD = pow(real(0.2) * (numeric_limits::epsilon() * real(0.01)),
181 1/real(8));
182 real
183 A0 = (x + y + 3*z)/5,
184 An = A0,
185 Q = fmax(fmax(fabs(A0-x), fabs(A0-y)), fabs(A0-z)) / tolRD,
186 x0 = x,
187 y0 = y,
188 z0 = z,
189 mul = 1,
190 s = 0;
191 while (Q >= mul * fabs(An)) {
192
193 real lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0);
194 s += 1/(mul * sqrt(z0) * (z0 + lam));
195 An = (An + lam)/4;
196 x0 = (x0 + lam)/4;
197 y0 = (y0 + lam)/4;
198 z0 = (z0 + lam)/4;
199 mul *= 4;
200 }
201 real
202 X = (A0 - x) / (mul * An),
203 Y = (A0 - y) / (mul * An),
204 Z = -(X + Y) / 3,
205 E2 = X*Y - 6*Z*Z,
206 E3 = (3*X*Y - 8*Z*Z)*Z,
207 E4 = 3 * (X*Y - Z*Z) * Z*Z,
208 E5 = X*Y*Z*Z*Z;
209
210
211
212
213
214 return ((471240 - 540540 * E2) * E5 +
215 (612612 * E2 - 540540 * E3 - 556920) * E4 +
216 E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
217 E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
218 (4084080 * mul * An * sqrt(An)) + 3 * s;
219 }
220
222 real kp2, real alphap2) {
223
224 if (k2 > 1)
225 throw GeographicErr("Parameter k2 is not in (-inf, 1]");
227 throw GeographicErr("Parameter alpha2 is not in (-inf, 1]");
228 if (kp2 < 0)
229 throw GeographicErr("Parameter kp2 is not in [0, inf)");
231 throw GeographicErr("Parameter alphap2 is not in [0, inf)");
232 _k2 = k2;
233 _kp2 = kp2;
236 _eps = _k2/Math::sq(sqrt(_kp2) + 1);
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259 if (_k2 != 0) {
260
261
263
264
265 _eEc = _kp2 != 0 ? 2 * RG(_kp2, 1) : 1;
266
267
269 } else {
270 _kKc = _eEc = Math::pi()/2; _dDc = _kKc/2;
271 }
272 if (_alpha2 != 0) {
273
274 real rj = (_kp2 != 0 && _alphap2 != 0) ? RJ(0, _kp2, 1, _alphap2) :
276
277 rc = _kp2 != 0 ? 0 :
279
280 _pPic = _kp2 != 0 ? _kKc + _alpha2 * rj / 3 : Math::infinity();
281
282 _gGc = _kp2 != 0 ? _kKc + (_alpha2 - _k2) * rj / 3 : rc;
283
284 _hHc = _kp2 != 0 ? _kKc - (_alphap2 != 0 ? _alphap2 * rj : 0) / 3 : rc;
285 } else {
286 _pPic = _kKc; _gGc = _eEc;
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302 _hHc = _kp2 == 1 ? Math::pi()/4 :
303 (_kp2 == 0 ? 1 : _kp2 * RD(0, 1, _kp2) / 3);
304 }
305 }
306
307
308
309
310
311
312
313
314
316
317 static const real tolJAC =
318 sqrt(numeric_limits::epsilon() * real(0.01));
319 if (_kp2 != 0) {
320 real mc = _kp2, d = 0;
321 if (signbit(_kp2)) {
322
323
324
325 d = 1 - mc;
326 mc /= -d;
327 d = sqrt(d);
328 x *= d;
329 }
330 real c = 0;
331 real m[num_], n[num_];
332 unsigned l = 0;
333 for (real a = 1;
334 l < num_ ||
336 ++l) {
337
338 m[l] = a;
339 n[l] = mc = sqrt(mc);
340 c = (a + mc) / 2;
341 if (!(fabs(a - mc) > tolJAC * a)) {
342 ++l;
343 break;
344 }
345 mc *= a;
346 a = c;
347 }
348 x *= c;
349 sn = sin(x);
350 cn = cos(x);
351 dn = 1;
352 if (sn != 0) {
353 real a = cn / sn;
354 c *= a;
355 while (l--) {
356 real b = m[l];
357 a *= c;
358 c *= dn;
359 dn = (n[l] + a) / (b + a);
360 a = c / b;
361 }
362 a = 1 / sqrt(c*c + 1);
363 sn = signbit(sn) ? -a : a;
364 cn = c * sn;
365 if (signbit(_kp2)) {
366
367 swap(cn, dn);
368 sn /= d;
369 }
370 }
371 } else {
372 sn = tanh(x);
373 dn = cn = 1 / cosh(x);
374 }
375 }
376
378
379
380 static const real tolJAC =
381 pow(numeric_limits::epsilon(), real(0.75));
382 real k2 = _k2, kp2 = _kp2;
383 if (_k2 == 0)
384 return x;
385 else if (_kp2 == 0) {
386 return atan(sinh(x));
387 } else if (_k2 < 0) {
388
389 k2 = -_k2 / _kp2; kp2 = 1 / _kp2;
390 x *= sqrt(_kp2);
391 }
392 real a[num_], b, c[num_];
393 a[0] = 1; b = sqrt(kp2); c[0] = sqrt(k2);
394 int l = 1;
395 for (; l < num_ ||
397 a[l] = (a[l-1] + b) / 2;
398 c[l] = (a[l-1] - b) / 2;
399 b = sqrt(a[l-1] * b);
400 if (!(c[l] > tolJAC * a[l])) break;
401 ++l;
402 }
403
404
405 real phi = a[l] * x * real(1 << l), phi1 = 0;
406 for (; l > 0; --l) {
407 phi1 = phi;
408 phi = (phi + asin(c[l] * sin(phi) / a[l])) / 2;
409 }
410
411 return _k2 < 0 ? phi1 - phi : phi;
412 }
413
415 real phi = am(x);
416 if (_kp2 == 0) {
417
418
419 sn = tanh(x); cn = dn = 1 / cosh(x);
420 } else {
421 sn = sin(phi); cn = cos(phi);
422
423
424 dn = Delta(sn, cn);
425 }
426 return phi;
427 }
428
430
431
432 real cn2 = cn*cn, dn2 = dn*dn,
433 fi = cn2 != 0 ? fabs(sn) * RF(cn2, dn2, 1) : K();
434
435 if (signbit(cn))
436 fi = 2 * K() - fi;
437 return copysign(fi, sn);
438 }
439
441 real
442 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
443 ei = cn2 != 0 ?
444 fabs(sn) * ( _k2 <= 0 ?
445
446
447 RF(cn2, dn2, 1) - _k2 * sn2 * RD(cn2, dn2, 1) / 3 :
448 ( _kp2 >= 0 ?
449
450 _kp2 * RF(cn2, dn2, 1) +
451 _k2 * _kp2 * sn2 * RD(cn2, 1, dn2) / 3 +
452 _k2 * fabs(cn) / dn :
453
454 - _kp2 * sn2 * RD(dn2, 1, cn2) / 3 +
455 dn / fabs(cn) ) ) :
456 E();
457
458 if (signbit(cn))
459 ei = 2 * E() - ei;
460 return copysign(ei, sn);
461 }
462
464
465
466 real
467 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
468 di = cn2 != 0 ? fabs(sn) * sn2 * RD(cn2, dn2, 1) / 3 : D();
469
470 if (signbit(cn))
471 di = 2 * D() - di;
472 return copysign(di, sn);
473 }
474
476
477
478 real
479 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
480 pii = cn2 != 0 ? fabs(sn) * (RF(cn2, dn2, 1) +
481 _alpha2 * sn2 *
482 RJ(cn2, dn2, 1, cn2 + _alphap2 * sn2) / 3) :
483 Pi();
484
485 if (signbit(cn))
486 pii = 2 * Pi() - pii;
487 return copysign(pii, sn);
488 }
489
491 real
492 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
493 gi = cn2 != 0 ? fabs(sn) * (RF(cn2, dn2, 1) +
494 (_alpha2 - _k2) * sn2 *
495 RJ(cn2, dn2, 1, cn2 + _alphap2 * sn2) / 3) :
496 G();
497
498 if (signbit(cn))
499 gi = 2 * G() - gi;
500 return copysign(gi, sn);
501 }
502
504 real
505 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
506
507 hi = cn2 != 0 ? fabs(sn) * (RF(cn2, dn2, 1) -
508 _alphap2 * sn2 *
509 RJ(cn2, dn2, 1, cn2 + _alphap2 * sn2) / 3) :
510 H();
511
512 if (signbit(cn))
513 hi = 2 * H() - hi;
514 return copysign(hi, sn);
515 }
516
518
519 if (signbit(cn)) { cn = -cn; sn = -sn; }
520 return F(sn, cn, dn) * (Math::pi()/2) / K() - atan2(sn, cn);
521 }
522
524
525 if (signbit(cn)) { cn = -cn; sn = -sn; }
526 return E(sn, cn, dn) * (Math::pi()/2) / E() - atan2(sn, cn);
527 }
528
530
531 if (signbit(cn)) { cn = -cn; sn = -sn; }
532 return Pi(sn, cn, dn) * (Math::pi()/2) / Pi() - atan2(sn, cn);
533 }
534
536
537 if (signbit(cn)) { cn = -cn; sn = -sn; }
538 return D(sn, cn, dn) * (Math::pi()/2) / D() - atan2(sn, cn);
539 }
540
542
543 if (signbit(cn)) { cn = -cn; sn = -sn; }
544 return G(sn, cn, dn) * (Math::pi()/2) / G() - atan2(sn, cn);
545 }
546
548
549 if (signbit(cn)) { cn = -cn; sn = -sn; }
550 return H(sn, cn, dn) * (Math::pi()/2) / H() - atan2(sn, cn);
551 }
552
554 if (_k2 == 0)
555 return phi;
556 else if (_kp2 == 0)
557 return asinh(tan(phi));
558 real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
559 return fabs(phi) < Math::pi() ? F(sn, cn, dn) :
561 }
562
564 if (_k2 == 0)
565 return phi;
566
567
568
569
570 real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
571 return fabs(phi) < Math::pi() ? E(sn, cn, dn) :
573 }
574
576
578 real sn, cn;
580 return E(sn, cn, Delta(sn, cn)) + 4 * E() * n;
581 }
582
584 real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
585 return fabs(phi) < Math::pi() ? Pi(sn, cn, dn) :
587 }
588
590 real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
591 return fabs(phi) < Math::pi() ? D(sn, cn, dn) :
593 }
594
596 real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
597 return fabs(phi) < Math::pi() ? G(sn, cn, dn) :
599 }
600
602 real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
603 return fabs(phi) < Math::pi() ? H(sn, cn, dn) :
605 }
606
608 static const real tolJAC =
609 sqrt(numeric_limits::epsilon() * real(0.01));
610 real n = floor(x / (2 * _eEc) + real(0.5));
611 x -= 2 * _eEc * n;
612
613 real phi = Math::pi() * x / (2 * _eEc);
614
615 phi -= _eps * sin(2 * phi) / 2;
616
617
618
619 for (int i = 0;
620 i < num_ ||
622 ++i) {
623 real
624 sn = sin(phi),
625 cn = cos(phi),
626 dn = Delta(sn, cn),
627 err = (E(sn, cn, dn) - x)/dn;
628 phi -= err;
629 if (!(fabs(err) > tolJAC))
630 break;
631 }
633 }
634
636
637 if (signbit(ctau)) { ctau = -ctau; stau = -stau; }
638 real tau = atan2(stau, ctau);
640 }
641
642}
Header for GeographicLib::EllipticFunction class.
#define GEOGRAPHICLIB_PANIC(msg)
void sncndn(real x, real &sn, real &cn, real &dn) const
Definition EllipticFunction.cpp:315
static real RJ(real x, real y, real z, real p)
Definition EllipticFunction.cpp:123
Math::real deltaG(real sn, real cn, real dn) const
Definition EllipticFunction.cpp:541
static real RG(real x, real y, real z)
Definition EllipticFunction.cpp:91
Math::real deltaE(real sn, real cn, real dn) const
Definition EllipticFunction.cpp:523
Math::real F(real phi) const
Definition EllipticFunction.cpp:553
static real RC(real x, real y)
Definition EllipticFunction.cpp:76
Math::real Einv(real x) const
Definition EllipticFunction.cpp:607
static real RD(real x, real y, real z)
Definition EllipticFunction.cpp:177
Math::real alphap2() const
void Reset(real k2=0, real alpha2=0)
Math::real am(real x) const
Definition EllipticFunction.cpp:377
Math::real Delta(real sn, real cn) const
Math::real deltaD(real sn, real cn, real dn) const
Definition EllipticFunction.cpp:535
Math::real Ed(real ang) const
Definition EllipticFunction.cpp:575
Math::real deltaH(real sn, real cn, real dn) const
Definition EllipticFunction.cpp:547
Math::real deltaF(real sn, real cn, real dn) const
Definition EllipticFunction.cpp:517
static real RF(real x, real y, real z)
Definition EllipticFunction.cpp:24
Math::real deltaPi(real sn, real cn, real dn) const
Definition EllipticFunction.cpp:529
Math::real deltaEinv(real stau, real ctau) const
Definition EllipticFunction.cpp:635
Math::real alpha2() const
Exception handling for GeographicLib.
static void sincosd(T x, T &sinx, T &cosx)
static T AngNormalize(T x)
static constexpr int td
degrees per turn
Namespace for GeographicLib.
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)