14 # pragma warning (disable: 4127) 31 static const real tolRF =
32 pow(3 * numeric_limits<real>::epsilon() *
real(0.01), 1/
real(8));
36 Q = max(max(abs(A0-x), abs(A0-y)), abs(A0-z)) / tolRF,
41 while (Q >= mul * abs(An)) {
43 real lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0);
51 X = (A0 - x) / (mul * An),
52 Y = (A0 - y) / (mul * An),
61 return (E3 * (6930 * E3 + E2 * (15015 * E2 - 16380) + 17160) +
62 E2 * ((10010 - 5775 * E2) * E2 - 24024) + 240240) /
68 static const real tolRG0 =
69 real(2.7) * sqrt((numeric_limits<real>::epsilon() *
real(0.01)));
70 real xn = sqrt(x), yn = sqrt(y);
71 if (xn < yn)
swap(xn, yn);
72 while (abs(xn-yn) > tolRG0 * xn) {
74 real t = (xn + yn) /2;
85 atan(sqrt((y - x) / x)) / sqrt(y - x) :
86 ( x == y ? 1 / sqrt(y) :
93 sqrt(-x / y) ) / sqrt(x - y) ) );
100 return (z * RF(x, y, z) - (x-z) * (y-z) * RD(x, y, z) / 3
101 + sqrt(x * y / z)) / 2;
106 static const real tolRG0 =
107 real(2.7) * sqrt((numeric_limits<real>::epsilon() *
real(0.01)));
109 x0 = sqrt(max(x, y)),
110 y0 = sqrt(min(x, y)),
115 while (abs(xn-yn) > tolRG0 * xn) {
117 real t = (xn + yn) /2;
130 tolRD = pow(
real(0.2) * (numeric_limits<real>::epsilon() *
real(0.01)),
133 A0 = (x + y + z + 2*p)/5,
135 delta = (p-x) * (p-y) * (p-z),
136 Q = max(max(abs(A0-x), abs(A0-y)), max(abs(A0-z), abs(A0-p))) / tolRD,
144 while (Q >= mul * abs(An)) {
147 lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0),
148 d0 = (sqrt(p0)+sqrt(x0)) * (sqrt(p0)+sqrt(y0)) * (sqrt(p0)+sqrt(z0)),
150 s += RC(1, 1 + e0)/(mul * d0);
160 X = (A0 - x) / (mul * An),
161 Y = (A0 - y) / (mul * An),
162 Z = (A0 - z) / (mul * An),
163 P = -(X + Y + Z) / 2,
164 E2 = X*Y + X*Z + Y*Z - 3*P*P,
165 E3 = X*Y*Z + 2*P * (E2 + 2*P*P),
166 E4 = (2*X*Y*Z + P * (E2 + 3*P*P)) * P,
173 return ((471240 - 540540 * E2) * E5 +
174 (612612 * E2 - 540540 * E3 - 556920) * E4 +
175 E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
176 E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
177 (4084080 * mul * An * sqrt(An)) + 6 * s;
183 tolRD = pow(
real(0.2) * (numeric_limits<real>::epsilon() *
real(0.01)),
186 A0 = (x + y + 3*z)/5,
188 Q = max(max(abs(A0-x), abs(A0-y)), abs(A0-z)) / tolRD,
194 while (Q >= mul * abs(An)) {
196 real lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0);
197 s += 1/(mul * sqrt(z0) * (z0 + lam));
205 X = (A0 - x) / (mul * An),
206 Y = (A0 - y) / (mul * An),
209 E3 = (3*X*Y - 8*Z*Z)*Z,
210 E4 = 3 * (X*Y - Z*Z) * Z*Z,
217 return ((471240 - 540540 * E2) * E5 +
218 (612612 * E2 - 540540 * E3 - 556920) * E4 +
219 E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
220 E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
221 (4084080 * mul * An * sqrt(An)) + 3 * s;
225 real kp2, real alphap2) {
239 _eps = _k2/
Math::sq(sqrt(_kp2) + 1);
268 _Ec = _kp2 != 0 ? 2 * RG(_kp2, 1) : 1;
273 _Kc = _Ec =
Math::pi()/2; _Dc = _Kc/2;
277 real rj = (_kp2 != 0 && _alphap2 != 0) ? RJ(0, _kp2, 1, _alphap2) :
285 _Gc = _kp2 != 0 ? _Kc + (_alpha2 - _k2) * rj / 3 : rc;
287 _Hc = _kp2 != 0 ? _Kc - (_alphap2 != 0 ? _alphap2 * rj : 0) / 3 : rc;
289 _Pic = _Kc; _Gc = _Ec;
303 _Hc = _kp2 != 0 ? _kp2 * RD(0, 1, _kp2) / 3 : 1;
317 static const real tolJAC =
318 sqrt(numeric_limits<real>::epsilon() *
real(0.01));
320 real mc = _kp2, d = 0;
328 real m[num_], n[num_];
333 n[l] = mc = sqrt(mc);
335 if (!(abs(a - mc) > tolJAC * a)) {
353 dn = (n[l] + a) / (b + a);
356 a = 1 / sqrt(c*c + 1);
357 sn = sn < 0 ? -a : a;
366 dn = cn = 1 / cosh(x);
373 real cn2 = cn*cn, dn2 = dn*dn,
374 fi = cn2 != 0 ? abs(sn) * RF(cn2, dn2, 1) : K();
378 return copysign(fi, sn);
383 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
385 abs(sn) * ( _k2 <= 0 ?
388 RF(cn2, dn2, 1) - _k2 * sn2 * RD(cn2, dn2, 1) / 3 :
391 _kp2 * RF(cn2, dn2, 1) +
392 _k2 * _kp2 * sn2 * RD(cn2, 1, dn2) / 3 +
395 - _kp2 * sn2 * RD(dn2, 1, cn2) / 3 +
401 return copysign(ei, sn);
408 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
409 di = cn2 != 0 ? abs(sn) * sn2 * RD(cn2, dn2, 1) / 3 : D();
413 return copysign(di, sn);
420 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
421 pii = cn2 != 0 ? abs(sn) * (RF(cn2, dn2, 1) +
423 RJ(cn2, dn2, 1, cn2 + _alphap2 * sn2) / 3) :
427 pii = 2 * Pi() - pii;
428 return copysign(pii, sn);
433 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
434 gi = cn2 != 0 ? abs(sn) * (RF(cn2, dn2, 1) +
435 (_alpha2 - _k2) * sn2 *
436 RJ(cn2, dn2, 1, cn2 + _alphap2 * sn2) / 3) :
441 return copysign(gi, sn);
446 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
448 hi = cn2 != 0 ? abs(sn) * (RF(cn2, dn2, 1) -
450 RJ(cn2, dn2, 1, cn2 + _alphap2 * sn2) / 3) :
455 return copysign(hi, sn);
460 if (cn < 0) { cn = -cn; sn = -sn; }
461 return F(sn, cn, dn) * (
Math::pi()/2) / K() - atan2(sn, cn);
466 if (cn < 0) { cn = -cn; sn = -sn; }
467 return E(sn, cn, dn) * (
Math::pi()/2) / E() - atan2(sn, cn);
472 if (cn < 0) { cn = -cn; sn = -sn; }
473 return Pi(sn, cn, dn) * (
Math::pi()/2) / Pi() - atan2(sn, cn);
478 if (cn < 0) { cn = -cn; sn = -sn; }
479 return D(sn, cn, dn) * (
Math::pi()/2) / D() - atan2(sn, cn);
484 if (cn < 0) { cn = -cn; sn = -sn; }
485 return G(sn, cn, dn) * (
Math::pi()/2) / G() - atan2(sn, cn);
490 if (cn < 0) { cn = -cn; sn = -sn; }
491 return H(sn, cn, dn) * (
Math::pi()/2) / H() - atan2(sn, cn);
495 real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
496 return abs(phi) <
Math::pi() ? F(sn, cn, dn) :
497 (deltaF(sn, cn, dn) + phi) * K() / (
Math::pi()/2);
501 real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
502 return abs(phi) <
Math::pi() ? E(sn, cn, dn) :
503 (deltaE(sn, cn, dn) + phi) * E() / (
Math::pi()/2);
507 real n = ceil(ang/360 -
real(0.5));
511 return E(sn, cn, Delta(sn, cn)) + 4 * E() * n;
515 real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
516 return abs(phi) <
Math::pi() ? Pi(sn, cn, dn) :
517 (deltaPi(sn, cn, dn) + phi) * Pi() / (
Math::pi()/2);
521 real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
522 return abs(phi) <
Math::pi() ? D(sn, cn, dn) :
523 (deltaD(sn, cn, dn) + phi) * D() / (
Math::pi()/2);
527 real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
528 return abs(phi) <
Math::pi() ? G(sn, cn, dn) :
529 (deltaG(sn, cn, dn) + phi) * G() / (
Math::pi()/2);
533 real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
534 return abs(phi) <
Math::pi() ? H(sn, cn, dn) :
535 (deltaH(sn, cn, dn) + phi) * H() / (
Math::pi()/2);
539 static const real tolJAC =
540 sqrt(numeric_limits<real>::epsilon() *
real(0.01));
541 real n = floor(x / (2 * _Ec) +
real(0.5));
544 real phi =
Math::pi() * x / (2 * _Ec);
546 phi -= _eps * sin(2 * phi) / 2;
555 err = (E(sn, cn, dn) - x)/dn;
557 if (!(abs(err) > tolJAC))
565 if (ctau < 0) { ctau = -ctau; stau = -stau; }
566 real tau = atan2(stau, ctau);
567 return Einv( tau * E() / (
Math::pi()/2) ) - tau;
void sncndn(real x, real &sn, real &cn, real &dn) const
void Reset(real k2=0, real alpha2=0)
GeographicLib::Math::real real
Math::real deltaPi(real sn, real cn, real dn) const
Math::real deltaE(real sn, real cn, real dn) const
static real RG(real x, real y, real z)
static real RF(real x, real y, real z)
static real RC(real x, real y)
Namespace for GeographicLib.
Header for GeographicLib::EllipticFunction class.
Math::real F(real phi) const
Math::real deltaH(real sn, real cn, real dn) const
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
Math::real deltaG(real sn, real cn, real dn) const
Math::real Ed(real ang) const
Math::real Einv(real x) const
Exception handling for GeographicLib.
Math::real deltaD(real sn, real cn, real dn) const
static real RD(real x, real y, real z)
static void sincosd(T x, T &sinx, T &cosx)
Math::real deltaEinv(real stau, real ctau) const
static real RJ(real x, real y, real z, real p)
#define GEOGRAPHICLIB_PANIC
Math::real deltaF(real sn, real cn, real dn) const