14 # pragma warning (disable: 4127) 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) /
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;
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;
129 real tolRD = pow(
real(0.2) * (numeric_limits<real>::epsilon() *
real(0.01)),
132 A0 = (x + y + z + 2*p)/5,
134 delta = (p-x) * (p-y) * (p-z),
135 Q = max(max(abs(A0-x), abs(A0-y)), max(abs(A0-z), abs(A0-p))) / tolRD,
143 while (Q >= mul * abs(An)) {
146 lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0),
147 d0 = (sqrt(p0)+sqrt(x0)) * (sqrt(p0)+sqrt(y0)) * (sqrt(p0)+sqrt(z0)),
149 s += RC(1, 1 + e0)/(mul * d0);
159 X = (A0 - x) / (mul * An),
160 Y = (A0 - y) / (mul * An),
161 Z = (A0 - z) / (mul * An),
162 P = -(X + Y + Z) / 2,
163 E2 = X*Y + X*Z + Y*Z - 3*P*P,
164 E3 = X*Y*Z + 2*P * (E2 + 2*P*P),
165 E4 = (2*X*Y*Z + P * (E2 + 3*P*P)) * P,
172 return ((471240 - 540540 * E2) * E5 +
173 (612612 * E2 - 540540 * E3 - 556920) * E4 +
174 E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
175 E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
176 (4084080 * mul * An * sqrt(An)) + 6 * s;
181 real tolRD = pow(
real(0.2) * (numeric_limits<real>::epsilon() *
real(0.01)),
184 A0 = (x + y + 3*z)/5,
186 Q = max(max(abs(A0-x), abs(A0-y)), abs(A0-z)) / tolRD,
192 while (Q >= mul * abs(An)) {
194 real lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0);
195 s += 1/(mul * sqrt(z0) * (z0 + lam));
203 X = (A0 - x) / (mul * An),
204 Y = (A0 - y) / (mul * An),
207 E3 = (3*X*Y - 8*Z*Z)*Z,
208 E4 = 3 * (X*Y - Z*Z) * Z*Z,
215 return ((471240 - 540540 * E2) * E5 +
216 (612612 * E2 - 540540 * E3 - 556920) * E4 +
217 E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
218 E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
219 (4084080 * mul * An * sqrt(An)) + 3 * s;
223 real kp2, real alphap2) {
228 _eps = _k2/
Math::sq(sqrt(_kp2) + 1);
235 _Ec = _kp2 ? 2 * RG(_kp2, 1) : 1;
240 _Kc = _Ec =
Math::pi()/2; _Dc = _Kc/2;
246 _Pic = _Kc + _alpha2 * rj / 3;
248 _Gc = _kp2 ? _Kc + (_alpha2 - _k2) * rj / 3 : RC(1, _alphap2);
250 _Hc = _kp2 ? _Kc - _alphap2 * rj / 3 : RC(1, _alphap2);
252 _Pic = _Kc; _Gc = _Ec; _Hc = _Kc - _Dc;
267 real tolJAC = sqrt(numeric_limits<real>::epsilon() *
real(0.01));
269 real mc = _kp2, d = 0;
277 real m[num_], n[num_];
282 n[l] = mc = sqrt(mc);
284 if (!(abs(a - mc) > tolJAC * a)) {
302 dn = (n[l] + a) / (b + a);
305 a = 1 / sqrt(c*c + 1);
306 sn = sn < 0 ? -a : a;
315 dn = cn = 1 / cosh(x);
322 real fi = abs(sn) * RF(cn*cn, dn*dn, 1);
333 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
337 RF(cn2, dn2, 1) - _k2 * sn2 * RD(cn2, dn2, 1) / 3 :
340 _kp2 * RF(cn2, dn2, 1) +
341 _k2 * _kp2 * sn2 * RD(cn2, 1, dn2) / 3 +
344 - _kp2 * sn2 * RD(dn2, 1, cn2) / 3 + dn / abs(cn) ) );
357 real di = abs(sn) * sn*sn * RD(cn*cn, dn*dn, 1) / 3;
370 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
371 pii = abs(sn) * (RF(cn2, dn2, 1) +
372 _alpha2 * sn2 * RJ(cn2, dn2, 1, 1 - _alpha2 * sn2) / 3);
375 pii = 2 * Pi() - pii;
383 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
384 gi = abs(sn) * (RF(cn2, dn2, 1) +
385 (_alpha2 - _k2) * sn2 *
386 RJ(cn2, dn2, 1, cn2 + _alphap2 * sn2) / 3);
397 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
398 hi = abs(sn) * (RF(cn2, dn2, 1) -
400 RJ(cn2, dn2, 1, cn2 + _alphap2 * sn2) / 3);
411 if (cn < 0) { cn = -cn; sn = -sn; }
412 return F(sn, cn, dn) * (
Math::pi()/2) / K() - atan2(sn, cn);
417 if (cn < 0) { cn = -cn; sn = -sn; }
418 return E(sn, cn, dn) * (
Math::pi()/2) / E() - atan2(sn, cn);
424 if (cn < 0) { cn = -cn; sn = -sn; }
425 return Pi(sn, cn, dn) * (
Math::pi()/2) / Pi() - atan2(sn, cn);
430 if (cn < 0) { cn = -cn; sn = -sn; }
431 return D(sn, cn, dn) * (
Math::pi()/2) / D() - atan2(sn, cn);
436 if (cn < 0) { cn = -cn; sn = -sn; }
437 return G(sn, cn, dn) * (
Math::pi()/2) / G() - atan2(sn, cn);
442 if (cn < 0) { cn = -cn; sn = -sn; }
443 return H(sn, cn, dn) * (
Math::pi()/2) / H() - atan2(sn, cn);
447 real sn = sin(phi), cn = cos(phi);
448 return (deltaF(sn, cn, Delta(sn, cn)) + phi) * K() / (
Math::pi()/2);
452 real sn = sin(phi), cn = cos(phi);
453 return (deltaE(sn, cn, Delta(sn, cn)) + phi) * E() / (
Math::pi()/2);
457 real n = ceil(ang/360 -
real(0.5));
461 sn = abs(ang) == 180 ? 0 : sin(phi),
462 cn = abs(ang) == 90 ? 0 : cos(phi);
463 return E(sn, cn, Delta(sn, cn)) + 4 * E() * n;
467 real sn = sin(phi), cn = cos(phi);
468 return (deltaPi(sn, cn, Delta(sn, cn)) + phi) * Pi() / (
Math::pi()/2);
472 real sn = sin(phi), cn = cos(phi);
473 return (deltaD(sn, cn, Delta(sn, cn)) + phi) * D() / (
Math::pi()/2);
477 real sn = sin(phi), cn = cos(phi);
478 return (deltaG(sn, cn, Delta(sn, cn)) + phi) * G() / (
Math::pi()/2);
482 real sn = sin(phi), cn = cos(phi);
483 return (deltaH(sn, cn, Delta(sn, cn)) + phi) * H() / (
Math::pi()/2);
487 real tolJAC = sqrt(numeric_limits<real>::epsilon() *
real(0.01));
488 real n = floor(x / (2 * _Ec) + 0.5);
491 real phi =
Math::pi() * x / (2 * _Ec);
493 phi -= _eps * sin(2 * phi) / 2;
499 err = (E(sn, cn, dn) - x)/dn;
501 if (abs(err) < tolJAC)
509 if (ctau < 0) { ctau = -ctau; stau = -stau; }
510 real tau = atan2(stau, ctau);
511 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
Math::real deltaG(real sn, real cn, real dn) const
Math::real Ed(real ang) const
Math::real Einv(real x) const
Math::real deltaD(real sn, real cn, real dn) const
static real RD(real x, real y, real z)
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