12 namespace GeographicLib {
17 numeric_limits<real>::epsilon() *
real(0.01);
18 const Math::real EllipticFunction::tolRF_ = pow(3 * tol_, 1/
real(8));
20 const Math::real EllipticFunction::tolRG0_ =
real(2.7) * sqrt(tol_);
21 const Math::real EllipticFunction::tolJAC_ = sqrt(tol_);
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 real xn = sqrt(x), yn = sqrt(y);
69 if (xn < yn) swap(xn, yn);
70 while (abs(xn-yn) > tolRG0_ * xn) {
72 real t = (xn + yn) /2;
76 return Math::pi<real>() / (xn + yn);
82 atan(sqrt((y - x) / x)) / sqrt(y - x) :
83 ( x == y && y > 0 ? 1 / sqrt(y) :
88 sqrt(x / (x - y)) ) / sqrt(x - y) ) );
95 return (z * RF(x, y, z) - (x-z) * (y-z) * RD(x, y, z) / 3
96 + sqrt(x * y / z)) / 2;
102 x0 = sqrt(max(x, y)),
103 y0 = sqrt(min(x, y)),
108 while (abs(xn-yn) > tolRG0_ * xn) {
110 real t = (xn + yn) /2;
117 return (
Math::sq( (x0 + y0)/2 ) - s) * Math::pi<real>() / (2 * (xn + yn));
123 A0 = (x + y + z + 2*p)/5,
125 delta = (p-x) * (p-y) * (p-z),
126 Q = max(max(abs(A0-x), abs(A0-y)), max(abs(A0-z), abs(A0-p))) / tolRD_,
134 while (Q >= mul * abs(An)) {
137 lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0),
138 d0 = (sqrt(p0)+sqrt(x0)) * (sqrt(p0)+sqrt(y0)) * (sqrt(p0)+sqrt(z0)),
140 s += RC(1, 1 + e0)/(mul * d0);
150 X = (A0 - x) / (mul * An),
151 Y = (A0 - y) / (mul * An),
152 Z = (A0 - z) / (mul * An),
153 P = -(X + Y + Z) / 2,
154 E2 = X*Y + X*Z + Y*Z - 3*P*P,
155 E3 = X*Y*Z + 2*P * (E2 + 2*P*P),
156 E4 = (2*X*Y*Z + P * (E2 + 3*P*P)) * P,
163 return ((471240 - 540540 * E2) * E5 +
164 (612612 * E2 - 540540 * E3 - 556920) * E4 +
165 E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
166 E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
167 (4084080 * mul * An * sqrt(An)) + 6 * s;
173 A0 = (x + y + 3*z)/5,
175 Q = max(max(abs(A0-x), abs(A0-y)), abs(A0-z)) / tolRD_,
181 while (Q >= mul * abs(An)) {
183 real lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0);
184 s += 1/(mul * sqrt(z0) * (z0 + lam));
192 X = (A0 - x) / (mul * An),
193 Y = (A0 - y) / (mul * An),
196 E3 = (3*X*Y - 8*Z*Z)*Z,
197 E4 = 3 * (X*Y - Z*Z) * Z*Z,
204 return ((471240 - 540540 * E2) * E5 +
205 (612612 * E2 - 540540 * E3 - 556920) * E4 +
206 E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
207 E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
208 (4084080 * mul * An * sqrt(An)) + 3 * s;
215 , _alphap2(1 - alpha2)
216 , _eps(_k2/
Math::sq(sqrt(_kp2) + 1))
223 real kp2, real alphap2)
throw()
228 , _eps(_k2/
Math::sq(sqrt(_kp2) + 1))
233 real kp2, real alphap2)
throw() {
238 _eps = _k2/
Math::sq(sqrt(_kp2) + 1);
242 bool EllipticFunction::Init()
const throw() {
245 _Kc = _kp2 ? RF(_kp2, 1) :
Math::infinity<
real>();
248 _Ec = _kp2 ? 2 * RG(_kp2, 1) : 1;
251 _Dc = _kp2 ? RD(
real(0), _kp2, 1) / 3 : Math::infinity<real>();
254 real rj = _kp2 ? RJ(0, _kp2, 1, _alphap2) :
Math::infinity<
real>();
256 _Pic = _Kc + _alpha2 * rj / 3;
258 _Gc = _kp2 ? _Kc + (_alpha2 - _k2) * rj / 3 : RC(1, _alphap2);
260 _Hc = _kp2 ? _Kc - _alphap2 * rj / 3 : RC(1, _alphap2);
262 _Pic = _Kc; _Gc = _Ec; _Hc = _Kc - _Dc;
279 real mc = _kp2, d = 0;
287 real m[num_], n[num_];
289 for (real a = 1; l < num_; ++l) {
292 n[l] = mc = sqrt(mc);
294 if (!(abs(a - mc) > tolJAC_ * a)) {
312 dn = (n[l] + a) / (b + a);
315 a = 1 / sqrt(c*c + 1);
316 sn = sn < 0 ? -a : a;
325 dn = cn = 1 / cosh(x);
332 real fi = abs(sn) * RF(cn*cn, dn*dn, 1);
343 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
347 RF(cn2, dn2, 1) - _k2 * sn2 * RD(cn2, dn2, 1) / 3 :
350 _kp2 * RF(cn2, dn2, 1) +
351 _k2 * _kp2 * sn2 * RD(cn2, 1, dn2) / 3 +
354 - _kp2 * sn2 * RD(dn2, 1, cn2) / 3 + dn / abs(cn) ) );
367 real di = abs(sn) * sn*sn * RD(cn*cn, dn*dn, 1) / 3;
380 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
381 pii = abs(sn) * (RF(cn2, dn2, 1) +
382 _alpha2 * sn2 * RJ(cn2, dn2, 1, 1 - _alpha2 * sn2) / 3);
385 pii = 2 * Pi() - pii;
393 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
394 gi = abs(sn) * (RF(cn2, dn2, 1) +
395 (_alpha2 - _k2) * sn2 *
396 RJ(cn2, dn2, 1, cn2 + _alphap2 * sn2) / 3);
407 cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
408 hi = abs(sn) * (RF(cn2, dn2, 1) -
410 RJ(cn2, dn2, 1, cn2 + _alphap2 * sn2) / 3);
421 if (cn < 0) { cn = -cn; sn = -sn; }
422 return F(sn, cn, dn) * (Math::pi<real>()/2) / K() - atan2(sn, cn);
427 if (cn < 0) { cn = -cn; sn = -sn; }
428 return E(sn, cn, dn) * (Math::pi<real>()/2) / E() - atan2(sn, cn);
434 if (cn < 0) { cn = -cn; sn = -sn; }
435 return Pi(sn, cn, dn) * (Math::pi<real>()/2) / Pi() - atan2(sn, cn);
440 if (cn < 0) { cn = -cn; sn = -sn; }
441 return D(sn, cn, dn) * (Math::pi<real>()/2) / D() - atan2(sn, cn);
446 if (cn < 0) { cn = -cn; sn = -sn; }
447 return G(sn, cn, dn) * (Math::pi<real>()/2) / G() - atan2(sn, cn);
452 if (cn < 0) { cn = -cn; sn = -sn; }
453 return H(sn, cn, dn) * (Math::pi<real>()/2) / H() - atan2(sn, cn);
457 real sn = sin(phi), cn = cos(phi);
458 return (deltaF(sn, cn, Delta(sn, cn)) + phi) * K() / (Math::pi<real>()/2);
462 real sn = sin(phi), cn = cos(phi);
463 return (deltaE(sn, cn, Delta(sn, cn)) + phi) * E() / (Math::pi<real>()/2);
467 real n = ceil(ang/360 -
real(0.5));
470 phi = ang * Math::degree<real>(),
471 sn = abs(ang) == 180 ? 0 : sin(phi),
472 cn = abs(ang) == 90 ? 0 : cos(phi);
473 return E(sn, cn, Delta(sn, cn)) + 4 * E() * n;
477 real sn = sin(phi), cn = cos(phi);
478 return (deltaPi(sn, cn, Delta(sn, cn)) + phi) * Pi() / (Math::pi<real>()/2);
482 real sn = sin(phi), cn = cos(phi);
483 return (deltaD(sn, cn, Delta(sn, cn)) + phi) * D() / (Math::pi<real>()/2);
487 real sn = sin(phi), cn = cos(phi);
488 return (deltaG(sn, cn, Delta(sn, cn)) + phi) * G() / (Math::pi<real>()/2);
492 real sn = sin(phi), cn = cos(phi);
493 return (deltaH(sn, cn, Delta(sn, cn)) + phi) * H() / (Math::pi<real>()/2);
498 real n = floor(x / (2 * _Ec) + 0.5);
501 real phi = Math::pi<real>() * x / (2 * _Ec);
503 phi -= _eps * sin(2 * phi) / 2;
504 for (
int i = 0; i < num_; ++i) {
509 err = (E(sn, cn, dn) - x)/dn;
511 if (abs(err) < tolJAC_)
514 return n * Math::pi<real>() + phi;
519 if (ctau < 0) { ctau = -ctau; stau = -stau; }
520 real tau = atan2(stau, ctau);
521 return Einv( tau * E() / (Math::pi<real>()/2) ) - tau;
Math::real deltaEinv(real stau, real ctau) const
void Reset(real k2=0, real alpha2=0)
Math::real deltaF(real sn, real cn, real dn) const
GeographicLib::Math::real real
Mathematical functions needed by GeographicLib.
Math::real deltaPi(real sn, real cn, real dn) const
static real RG(real x, real y, real z)
Math::real Ed(real ang) const
Math::real Einv(real x) const
static real RF(real x, real y, real z)
static real RC(real x, real y)
Math::real F(real phi) const
void sncndn(real x, real &sn, real &cn, real &dn) const
Header for GeographicLib::EllipticFunction class.
Math::real deltaE(real sn, real cn, real dn) const
Math::real deltaH(real sn, real cn, real dn) const
static real RD(real x, real y, real z)
static real RJ(real x, real y, real z, real p)
EllipticFunction(real k2=0, real alpha2=0)
Math::real deltaD(real sn, real cn, real dn) const
Math::real deltaG(real sn, real cn, real dn) const