GeographicLib  2.1.2
EllipticFunction.cpp
Go to the documentation of this file.
1 /**
2  * \file EllipticFunction.cpp
3  * \brief Implementation for GeographicLib::EllipticFunction class
4  *
5  * Copyright (c) Charles Karney (2008-2022) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * https://geographiclib.sourceforge.io/
8  **********************************************************************/
9 
11 
12 #if defined(_MSC_VER)
13 // Squelch warnings about constant conditional and enum-float expressions
14 # pragma warning (disable: 4127 5055)
15 #endif
16 
17 namespace GeographicLib {
18 
19  using namespace std;
20 
21  /*
22  * Implementation of methods given in
23  *
24  * B. C. Carlson
25  * Computation of elliptic integrals
26  * Numerical Algorithms 10, 13-26 (1995)
27  */
28 
29  Math::real EllipticFunction::RF(real x, real y, real z) {
30  // Carlson, eqs 2.2 - 2.7
31  static const real tolRF =
32  pow(3 * numeric_limits<real>::epsilon() * real(0.01), 1/real(8));
33  real
34  A0 = (x + y + z)/3,
35  An = A0,
36  Q = fmax(fmax(fabs(A0-x), fabs(A0-y)), fabs(A0-z)) / tolRF,
37  x0 = x,
38  y0 = y,
39  z0 = z,
40  mul = 1;
41  while (Q >= mul * fabs(An)) {
42  // Max 6 trips
43  real lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0);
44  An = (An + lam)/4;
45  x0 = (x0 + lam)/4;
46  y0 = (y0 + lam)/4;
47  z0 = (z0 + lam)/4;
48  mul *= 4;
49  }
50  real
51  X = (A0 - x) / (mul * An),
52  Y = (A0 - y) / (mul * An),
53  Z = - (X + Y),
54  E2 = X*Y - Z*Z,
55  E3 = X*Y*Z;
56  // https://dlmf.nist.gov/19.36.E1
57  // Polynomial is
58  // (1 - E2/10 + E3/14 + E2^2/24 - 3*E2*E3/44
59  // - 5*E2^3/208 + 3*E3^2/104 + E2^2*E3/16)
60  // convert to Horner form...
61  return (E3 * (6930 * E3 + E2 * (15015 * E2 - 16380) + 17160) +
62  E2 * ((10010 - 5775 * E2) * E2 - 24024) + 240240) /
63  (240240 * sqrt(An));
64  }
65 
66  Math::real EllipticFunction::RF(real x, real y) {
67  // Carlson, eqs 2.36 - 2.38
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 (fabs(xn-yn) > tolRG0 * xn) {
73  // Max 4 trips
74  real t = (xn + yn) /2;
75  yn = sqrt(xn * yn);
76  xn = t;
77  }
78  return Math::pi() / (xn + yn);
79  }
80 
81  Math::real EllipticFunction::RC(real x, real y) {
82  // Defined only for y != 0 and x >= 0.
83  return ( !(x >= y) ? // x < y and catch nans
84  // https://dlmf.nist.gov/19.2.E18
85  atan(sqrt((y - x) / x)) / sqrt(y - x) :
86  ( x == y ? 1 / sqrt(y) :
87  asinh( y > 0 ?
88  // https://dlmf.nist.gov/19.2.E19
89  // atanh(sqrt((x - y) / x))
90  sqrt((x - y) / y) :
91  // https://dlmf.nist.gov/19.2.E20
92  // atanh(sqrt(x / (x - y)))
93  sqrt(-x / y) ) / sqrt(x - y) ) );
94  }
95 
96  Math::real EllipticFunction::RG(real x, real y, real z) {
97  return (x == 0 ? RG(y, z) :
98  (y == 0 ? RG(z, x) :
99  (z == 0 ? RG(x, y) :
100  // Carlson, eq 1.7
101  (z * RF(x, y, z) - (x-z) * (y-z) * RD(x, y, z) / 3
102  + sqrt(x * y / z)) / 2 )));
103  }
104 
106  // Carlson, eqs 2.36 - 2.39
107  static const real tolRG0 =
108  real(2.7) * sqrt((numeric_limits<real>::epsilon() * real(0.01)));
109  real
110  x0 = sqrt(fmax(x, y)),
111  y0 = sqrt(fmin(x, y)),
112  xn = x0,
113  yn = y0,
114  s = 0,
115  mul = real(0.25);
116  while (fabs(xn-yn) > tolRG0 * xn) {
117  // Max 4 trips
118  real t = (xn + yn) /2;
119  yn = sqrt(xn * yn);
120  xn = t;
121  mul *= 2;
122  t = xn - yn;
123  s += mul * t * t;
124  }
125  return (Math::sq( (x0 + y0)/2 ) - s) * Math::pi() / (2 * (xn + yn));
126  }
127 
128  Math::real EllipticFunction::RJ(real x, real y, real z, real p) {
129  // Carlson, eqs 2.17 - 2.25
130  static const real
131  tolRD = pow(real(0.2) * (numeric_limits<real>::epsilon() * real(0.01)),
132  1/real(8));
133  real
134  A0 = (x + y + z + 2*p)/5,
135  An = A0,
136  delta = (p-x) * (p-y) * (p-z),
137  Q = fmax(fmax(fabs(A0-x), fabs(A0-y)),
138  fmax(fabs(A0-z), fabs(A0-p))) / tolRD,
139  x0 = x,
140  y0 = y,
141  z0 = z,
142  p0 = p,
143  mul = 1,
144  mul3 = 1,
145  s = 0;
146  while (Q >= mul * fabs(An)) {
147  // Max 7 trips
148  real
149  lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0),
150  d0 = (sqrt(p0)+sqrt(x0)) * (sqrt(p0)+sqrt(y0)) * (sqrt(p0)+sqrt(z0)),
151  e0 = delta/(mul3 * Math::sq(d0));
152  s += RC(1, 1 + e0)/(mul * d0);
153  An = (An + lam)/4;
154  x0 = (x0 + lam)/4;
155  y0 = (y0 + lam)/4;
156  z0 = (z0 + lam)/4;
157  p0 = (p0 + lam)/4;
158  mul *= 4;
159  mul3 *= 64;
160  }
161  real
162  X = (A0 - x) / (mul * An),
163  Y = (A0 - y) / (mul * An),
164  Z = (A0 - z) / (mul * An),
165  P = -(X + Y + Z) / 2,
166  E2 = X*Y + X*Z + Y*Z - 3*P*P,
167  E3 = X*Y*Z + 2*P * (E2 + 2*P*P),
168  E4 = (2*X*Y*Z + P * (E2 + 3*P*P)) * P,
169  E5 = X*Y*Z*P*P;
170  // https://dlmf.nist.gov/19.36.E2
171  // Polynomial is
172  // (1 - 3*E2/14 + E3/6 + 9*E2^2/88 - 3*E4/22 - 9*E2*E3/52 + 3*E5/26
173  // - E2^3/16 + 3*E3^2/40 + 3*E2*E4/20 + 45*E2^2*E3/272
174  // - 9*(E3*E4+E2*E5)/68)
175  return ((471240 - 540540 * E2) * E5 +
176  (612612 * E2 - 540540 * E3 - 556920) * E4 +
177  E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
178  E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
179  (4084080 * mul * An * sqrt(An)) + 6 * s;
180  }
181 
182  Math::real EllipticFunction::RD(real x, real y, real z) {
183  // Carlson, eqs 2.28 - 2.34
184  static const real
185  tolRD = pow(real(0.2) * (numeric_limits<real>::epsilon() * real(0.01)),
186  1/real(8));
187  real
188  A0 = (x + y + 3*z)/5,
189  An = A0,
190  Q = fmax(fmax(fabs(A0-x), fabs(A0-y)), fabs(A0-z)) / tolRD,
191  x0 = x,
192  y0 = y,
193  z0 = z,
194  mul = 1,
195  s = 0;
196  while (Q >= mul * fabs(An)) {
197  // Max 7 trips
198  real lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0);
199  s += 1/(mul * sqrt(z0) * (z0 + lam));
200  An = (An + lam)/4;
201  x0 = (x0 + lam)/4;
202  y0 = (y0 + lam)/4;
203  z0 = (z0 + lam)/4;
204  mul *= 4;
205  }
206  real
207  X = (A0 - x) / (mul * An),
208  Y = (A0 - y) / (mul * An),
209  Z = -(X + Y) / 3,
210  E2 = X*Y - 6*Z*Z,
211  E3 = (3*X*Y - 8*Z*Z)*Z,
212  E4 = 3 * (X*Y - Z*Z) * Z*Z,
213  E5 = X*Y*Z*Z*Z;
214  // https://dlmf.nist.gov/19.36.E2
215  // Polynomial is
216  // (1 - 3*E2/14 + E3/6 + 9*E2^2/88 - 3*E4/22 - 9*E2*E3/52 + 3*E5/26
217  // - E2^3/16 + 3*E3^2/40 + 3*E2*E4/20 + 45*E2^2*E3/272
218  // - 9*(E3*E4+E2*E5)/68)
219  return ((471240 - 540540 * E2) * E5 +
220  (612612 * E2 - 540540 * E3 - 556920) * E4 +
221  E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
222  E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
223  (4084080 * mul * An * sqrt(An)) + 3 * s;
224  }
225 
226  void EllipticFunction::Reset(real k2, real alpha2,
227  real kp2, real alphap2) {
228  // Accept nans here (needed for GeodesicExact)
229  if (k2 > 1)
230  throw GeographicErr("Parameter k2 is not in (-inf, 1]");
231  if (alpha2 > 1)
232  throw GeographicErr("Parameter alpha2 is not in (-inf, 1]");
233  if (kp2 < 0)
234  throw GeographicErr("Parameter kp2 is not in [0, inf)");
235  if (alphap2 < 0)
236  throw GeographicErr("Parameter alphap2 is not in [0, inf)");
237  _k2 = k2;
238  _kp2 = kp2;
239  _alpha2 = alpha2;
240  _alphap2 = alphap2;
241  _eps = _k2/Math::sq(sqrt(_kp2) + 1);
242  // Values of complete elliptic integrals for k = 0,1 and alpha = 0,1
243  // K E D
244  // k = 0: pi/2 pi/2 pi/4
245  // k = 1: inf 1 inf
246  // Pi G H
247  // k = 0, alpha = 0: pi/2 pi/2 pi/4
248  // k = 1, alpha = 0: inf 1 1
249  // k = 0, alpha = 1: inf inf pi/2
250  // k = 1, alpha = 1: inf inf inf
251  //
252  // Pi(0, k) = K(k)
253  // G(0, k) = E(k)
254  // H(0, k) = K(k) - D(k)
255  // Pi(0, k) = K(k)
256  // G(0, k) = E(k)
257  // H(0, k) = K(k) - D(k)
258  // Pi(alpha2, 0) = pi/(2*sqrt(1-alpha2))
259  // G(alpha2, 0) = pi/(2*sqrt(1-alpha2))
260  // H(alpha2, 0) = pi/(2*(1 + sqrt(1-alpha2)))
261  // Pi(alpha2, 1) = inf
262  // H(1, k) = K(k)
263  // G(alpha2, 1) = H(alpha2, 1) = RC(1, alphap2)
264  if (_k2 != 0) {
265  // Complete elliptic integral K(k), Carlson eq. 4.1
266  // https://dlmf.nist.gov/19.25.E1
267  _kKc = _kp2 != 0 ? RF(_kp2, 1) : Math::infinity();
268  // Complete elliptic integral E(k), Carlson eq. 4.2
269  // https://dlmf.nist.gov/19.25.E1
270  _eEc = _kp2 != 0 ? 2 * RG(_kp2, 1) : 1;
271  // D(k) = (K(k) - E(k))/k^2, Carlson eq.4.3
272  // https://dlmf.nist.gov/19.25.E1
273  _dDc = _kp2 != 0 ? RD(0, _kp2, 1) / 3 : Math::infinity();
274  } else {
275  _kKc = _eEc = Math::pi()/2; _dDc = _kKc/2;
276  }
277  if (_alpha2 != 0) {
278  // https://dlmf.nist.gov/19.25.E2
279  real rj = (_kp2 != 0 && _alphap2 != 0) ? RJ(0, _kp2, 1, _alphap2) :
280  Math::infinity(),
281  // Only use rc if _kp2 = 0.
282  rc = _kp2 != 0 ? 0 :
283  (_alphap2 != 0 ? RC(1, _alphap2) : Math::infinity());
284  // Pi(alpha^2, k)
285  _pPic = _kp2 != 0 ? _kKc + _alpha2 * rj / 3 : Math::infinity();
286  // G(alpha^2, k)
287  _gGc = _kp2 != 0 ? _kKc + (_alpha2 - _k2) * rj / 3 : rc;
288  // H(alpha^2, k)
289  _hHc = _kp2 != 0 ? _kKc - (_alphap2 != 0 ? _alphap2 * rj : 0) / 3 : rc;
290  } else {
291  _pPic = _kKc; _gGc = _eEc;
292  // Hc = Kc - Dc but this involves large cancellations if k2 is close to
293  // 1. So write (for alpha2 = 0)
294  // Hc = int(cos(phi)^2/sqrt(1-k2*sin(phi)^2),phi,0,pi/2)
295  // = 1/sqrt(1-k2) * int(sin(phi)^2/sqrt(1-k2/kp2*sin(phi)^2,...)
296  // = 1/kp * D(i*k/kp)
297  // and use D(k) = RD(0, kp2, 1) / 3
298  // so Hc = 1/kp * RD(0, 1/kp2, 1) / 3
299  // = kp2 * RD(0, 1, kp2) / 3
300  // using https://dlmf.nist.gov/19.20.E18
301  // Equivalently
302  // RF(x, 1) - RD(0, x, 1)/3 = x * RD(0, 1, x)/3 for x > 0
303  // For k2 = 1 and alpha2 = 0, we have
304  // Hc = int(cos(phi),...) = 1
305  _hHc = _kp2 != 0 ? _kp2 * RD(0, 1, _kp2) / 3 : 1;
306  }
307  }
308 
309  /*
310  * Implementation of methods given in
311  *
312  * R. Bulirsch
313  * Numerical Calculation of Elliptic Integrals and Elliptic Functions
314  * Numericshe Mathematik 7, 78-90 (1965)
315  */
316 
317  void EllipticFunction::sncndn(real x, real& sn, real& cn, real& dn) const {
318  // Bulirsch's sncndn routine, p 89.
319  static const real tolJAC =
320  sqrt(numeric_limits<real>::epsilon() * real(0.01));
321  if (_kp2 != 0) {
322  real mc = _kp2, d = 0;
323  if (signbit(_kp2)) {
324  d = 1 - mc;
325  mc /= -d;
326  d = sqrt(d);
327  x *= d;
328  }
329  real c = 0; // To suppress warning about uninitialized variable
330  real m[num_], n[num_];
331  unsigned l = 0;
332  for (real a = 1; l < num_ || GEOGRAPHICLIB_PANIC; ++l) {
333  // This converges quadratically. Max 5 trips
334  m[l] = a;
335  n[l] = mc = sqrt(mc);
336  c = (a + mc) / 2;
337  if (!(fabs(a - mc) > tolJAC * a)) {
338  ++l;
339  break;
340  }
341  mc *= a;
342  a = c;
343  }
344  x *= c;
345  sn = sin(x);
346  cn = cos(x);
347  dn = 1;
348  if (sn != 0) {
349  real a = cn / sn;
350  c *= a;
351  while (l--) {
352  real b = m[l];
353  a *= c;
354  c *= dn;
355  dn = (n[l] + a) / (b + a);
356  a = c / b;
357  }
358  a = 1 / sqrt(c*c + 1);
359  sn = signbit(sn) ? -a : a;
360  cn = c * sn;
361  if (signbit(_kp2)) {
362  swap(cn, dn);
363  sn /= d;
364  }
365  }
366  } else {
367  sn = tanh(x);
368  dn = cn = 1 / cosh(x);
369  }
370  }
371 
372  Math::real EllipticFunction::F(real sn, real cn, real dn) const {
373  // Carlson, eq. 4.5 and
374  // https://dlmf.nist.gov/19.25.E5
375  real cn2 = cn*cn, dn2 = dn*dn,
376  fi = cn2 != 0 ? fabs(sn) * RF(cn2, dn2, 1) : K();
377  // Enforce usual trig-like symmetries
378  if (signbit(cn))
379  fi = 2 * K() - fi;
380  return copysign(fi, sn);
381  }
382 
383  Math::real EllipticFunction::E(real sn, real cn, real dn) const {
384  real
385  cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
386  ei = cn2 != 0 ?
387  fabs(sn) * ( _k2 <= 0 ?
388  // Carlson, eq. 4.6 and
389  // https://dlmf.nist.gov/19.25.E9
390  RF(cn2, dn2, 1) - _k2 * sn2 * RD(cn2, dn2, 1) / 3 :
391  ( _kp2 >= 0 ?
392  // https://dlmf.nist.gov/19.25.E10
393  _kp2 * RF(cn2, dn2, 1) +
394  _k2 * _kp2 * sn2 * RD(cn2, 1, dn2) / 3 +
395  _k2 * fabs(cn) / dn :
396  // https://dlmf.nist.gov/19.25.E11
397  - _kp2 * sn2 * RD(dn2, 1, cn2) / 3 +
398  dn / fabs(cn) ) ) :
399  E();
400  // Enforce usual trig-like symmetries
401  if (signbit(cn))
402  ei = 2 * E() - ei;
403  return copysign(ei, sn);
404  }
405 
406  Math::real EllipticFunction::D(real sn, real cn, real dn) const {
407  // Carlson, eq. 4.8 and
408  // https://dlmf.nist.gov/19.25.E13
409  real
410  cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
411  di = cn2 != 0 ? fabs(sn) * sn2 * RD(cn2, dn2, 1) / 3 : D();
412  // Enforce usual trig-like symmetries
413  if (signbit(cn))
414  di = 2 * D() - di;
415  return copysign(di, sn);
416  }
417 
418  Math::real EllipticFunction::Pi(real sn, real cn, real dn) const {
419  // Carlson, eq. 4.7 and
420  // https://dlmf.nist.gov/19.25.E14
421  real
422  cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
423  pii = cn2 != 0 ? fabs(sn) * (RF(cn2, dn2, 1) +
424  _alpha2 * sn2 *
425  RJ(cn2, dn2, 1, cn2 + _alphap2 * sn2) / 3) :
426  Pi();
427  // Enforce usual trig-like symmetries
428  if (signbit(cn))
429  pii = 2 * Pi() - pii;
430  return copysign(pii, sn);
431  }
432 
433  Math::real EllipticFunction::G(real sn, real cn, real dn) const {
434  real
435  cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
436  gi = cn2 != 0 ? fabs(sn) * (RF(cn2, dn2, 1) +
437  (_alpha2 - _k2) * sn2 *
438  RJ(cn2, dn2, 1, cn2 + _alphap2 * sn2) / 3) :
439  G();
440  // Enforce usual trig-like symmetries
441  if (signbit(cn))
442  gi = 2 * G() - gi;
443  return copysign(gi, sn);
444  }
445 
446  Math::real EllipticFunction::H(real sn, real cn, real dn) const {
447  real
448  cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
449  // WARNING: large cancellation if k2 = 1, alpha2 = 0, and phi near pi/2
450  hi = cn2 != 0 ? fabs(sn) * (RF(cn2, dn2, 1) -
451  _alphap2 * sn2 *
452  RJ(cn2, dn2, 1, cn2 + _alphap2 * sn2) / 3) :
453  H();
454  // Enforce usual trig-like symmetries
455  if (signbit(cn))
456  hi = 2 * H() - hi;
457  return copysign(hi, sn);
458  }
459 
460  Math::real EllipticFunction::deltaF(real sn, real cn, real dn) const {
461  // Function is periodic with period pi
462  if (signbit(cn)) { cn = -cn; sn = -sn; }
463  return F(sn, cn, dn) * (Math::pi()/2) / K() - atan2(sn, cn);
464  }
465 
466  Math::real EllipticFunction::deltaE(real sn, real cn, real dn) const {
467  // Function is periodic with period pi
468  if (signbit(cn)) { cn = -cn; sn = -sn; }
469  return E(sn, cn, dn) * (Math::pi()/2) / E() - atan2(sn, cn);
470  }
471 
472  Math::real EllipticFunction::deltaPi(real sn, real cn, real dn) const {
473  // Function is periodic with period pi
474  if (signbit(cn)) { cn = -cn; sn = -sn; }
475  return Pi(sn, cn, dn) * (Math::pi()/2) / Pi() - atan2(sn, cn);
476  }
477 
478  Math::real EllipticFunction::deltaD(real sn, real cn, real dn) const {
479  // Function is periodic with period pi
480  if (signbit(cn)) { cn = -cn; sn = -sn; }
481  return D(sn, cn, dn) * (Math::pi()/2) / D() - atan2(sn, cn);
482  }
483 
484  Math::real EllipticFunction::deltaG(real sn, real cn, real dn) const {
485  // Function is periodic with period pi
486  if (signbit(cn)) { cn = -cn; sn = -sn; }
487  return G(sn, cn, dn) * (Math::pi()/2) / G() - atan2(sn, cn);
488  }
489 
490  Math::real EllipticFunction::deltaH(real sn, real cn, real dn) const {
491  // Function is periodic with period pi
492  if (signbit(cn)) { cn = -cn; sn = -sn; }
493  return H(sn, cn, dn) * (Math::pi()/2) / H() - atan2(sn, cn);
494  }
495 
496  Math::real EllipticFunction::F(real phi) const {
497  real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
498  return fabs(phi) < Math::pi() ? F(sn, cn, dn) :
499  (deltaF(sn, cn, dn) + phi) * K() / (Math::pi()/2);
500  }
501 
502  Math::real EllipticFunction::E(real phi) const {
503  real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
504  return fabs(phi) < Math::pi() ? E(sn, cn, dn) :
505  (deltaE(sn, cn, dn) + phi) * E() / (Math::pi()/2);
506  }
507 
509  // ang - Math::AngNormalize(ang) is (nearly) an exact multiple of 360
510  real n = round((ang - Math::AngNormalize(ang))/Math::td);
511  real sn, cn;
512  Math::sincosd(ang, sn, cn);
513  return E(sn, cn, Delta(sn, cn)) + 4 * E() * n;
514  }
515 
517  real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
518  return fabs(phi) < Math::pi() ? Pi(sn, cn, dn) :
519  (deltaPi(sn, cn, dn) + phi) * Pi() / (Math::pi()/2);
520  }
521 
522  Math::real EllipticFunction::D(real phi) const {
523  real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
524  return fabs(phi) < Math::pi() ? D(sn, cn, dn) :
525  (deltaD(sn, cn, dn) + phi) * D() / (Math::pi()/2);
526  }
527 
528  Math::real EllipticFunction::G(real phi) const {
529  real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
530  return fabs(phi) < Math::pi() ? G(sn, cn, dn) :
531  (deltaG(sn, cn, dn) + phi) * G() / (Math::pi()/2);
532  }
533 
534  Math::real EllipticFunction::H(real phi) const {
535  real sn = sin(phi), cn = cos(phi), dn = Delta(sn, cn);
536  return fabs(phi) < Math::pi() ? H(sn, cn, dn) :
537  (deltaH(sn, cn, dn) + phi) * H() / (Math::pi()/2);
538  }
539 
541  static const real tolJAC =
542  sqrt(numeric_limits<real>::epsilon() * real(0.01));
543  real n = floor(x / (2 * _eEc) + real(0.5));
544  x -= 2 * _eEc * n; // x now in [-ec, ec)
545  // Linear approximation
546  real phi = Math::pi() * x / (2 * _eEc); // phi in [-pi/2, pi/2)
547  // First order correction
548  phi -= _eps * sin(2 * phi) / 2;
549  // For kp2 close to zero use asin(x/_eEc) or
550  // J. P. Boyd, Applied Math. and Computation 218, 7005-7013 (2012)
551  // https://doi.org/10.1016/j.amc.2011.12.021
552  for (int i = 0; i < num_ || GEOGRAPHICLIB_PANIC; ++i) {
553  real
554  sn = sin(phi),
555  cn = cos(phi),
556  dn = Delta(sn, cn),
557  err = (E(sn, cn, dn) - x)/dn;
558  phi -= err;
559  if (!(fabs(err) > tolJAC))
560  break;
561  }
562  return n * Math::pi() + phi;
563  }
564 
565  Math::real EllipticFunction::deltaEinv(real stau, real ctau) const {
566  // Function is periodic with period pi
567  if (signbit(ctau)) { ctau = -ctau; stau = -stau; }
568  real tau = atan2(stau, ctau);
569  return Einv( tau * E() / (Math::pi()/2) ) - tau;
570  }
571 
572 } // namespace GeographicLib
Header for GeographicLib::EllipticFunction class.
GeographicLib::Math::real real
Definition: GeodSolve.cpp:31
#define GEOGRAPHICLIB_PANIC
Definition: Math.hpp:61
void sncndn(real x, real &sn, real &cn, real &dn) const
static real RJ(real x, real y, real z, real p)
Math::real deltaG(real sn, real cn, real dn) const
static real RG(real x, real y, real z)
Math::real deltaE(real sn, real cn, real dn) const
Math::real F(real phi) const
static real RC(real x, real y)
Math::real Einv(real x) const
static real RD(real x, real y, real z)
void Reset(real k2=0, real alpha2=0)
Math::real deltaD(real sn, real cn, real dn) const
Math::real Ed(real ang) const
Math::real deltaH(real sn, real cn, real dn) const
Math::real deltaF(real sn, real cn, real dn) const
static real RF(real x, real y, real z)
Math::real deltaPi(real sn, real cn, real dn) const
Math::real deltaEinv(real stau, real ctau) const
Exception handling for GeographicLib.
Definition: Constants.hpp:316
static void sincosd(T x, T &sinx, T &cosx)
Definition: Math.cpp:106
static T sq(T x)
Definition: Math.hpp:212
static T AngNormalize(T x)
Definition: Math.cpp:71
static T infinity()
Definition: Math.cpp:262
static T pi()
Definition: Math.hpp:190
@ td
degrees per turn
Definition: Math.hpp:145
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)