GeographicLib  1.37
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-2012) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
8  **********************************************************************/
9 
11 
12 #if defined(_MSC_VER)
13 // Squelch warnings about constant conditional expressions
14 # pragma warning (disable: 4127)
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  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 = max(max(abs(A0-x), abs(A0-y)), abs(A0-z)) / tolRF,
37  x0 = x,
38  y0 = y,
39  z0 = z,
40  mul = 1;
41  while (Q >= mul * abs(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  // http://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  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) {
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  return ( !(x >= y) ? // x < y and catch nans
83  // http://dlmf.nist.gov/19.2.E18
84  atan(sqrt((y - x) / x)) / sqrt(y - x) :
85  ( x == y && y > 0 ? 1 / sqrt(y) :
86  Math::atanh( y > 0 ?
87  // http://dlmf.nist.gov/19.2.E19
88  sqrt((x - y) / x) :
89  // http://dlmf.nist.gov/19.2.E20
90  sqrt(x / (x - y)) ) / sqrt(x - y) ) );
91  }
92 
93  Math::real EllipticFunction::RG(real x, real y, real z) {
94  if (z == 0)
95  swap(y, z);
96  // Carlson, eq 1.7
97  return (z * RF(x, y, z) - (x-z) * (y-z) * RD(x, y, z) / 3
98  + sqrt(x * y / z)) / 2;
99  }
100 
102  // Carlson, eqs 2.36 - 2.39
103  real tolRG0 =
104  real(2.7) * sqrt((numeric_limits<real>::epsilon() * real(0.01)));
105  real
106  x0 = sqrt(max(x, y)),
107  y0 = sqrt(min(x, y)),
108  xn = x0,
109  yn = y0,
110  s = 0,
111  mul = real(0.25);
112  while (abs(xn-yn) > tolRG0 * xn) {
113  // Max 4 trips
114  real t = (xn + yn) /2;
115  yn = sqrt(xn * yn);
116  xn = t;
117  mul *= 2;
118  t = xn - yn;
119  s += mul * t * t;
120  }
121  return (Math::sq( (x0 + y0)/2 ) - s) * Math::pi() / (2 * (xn + yn));
122  }
123 
124  Math::real EllipticFunction::RJ(real x, real y, real z, real p) {
125  // Carlson, eqs 2.17 - 2.25
126  real tolRD = pow(real(0.2) * (numeric_limits<real>::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 = max(max(abs(A0-x), abs(A0-y)), max(abs(A0-z), abs(A0-p))) / tolRD,
133  x0 = x,
134  y0 = y,
135  z0 = z,
136  p0 = p,
137  mul = 1,
138  mul3 = 1,
139  s = 0;
140  while (Q >= mul * abs(An)) {
141  // Max 7 trips
142  real
143  lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0),
144  d0 = (sqrt(p0)+sqrt(x0)) * (sqrt(p0)+sqrt(y0)) * (sqrt(p0)+sqrt(z0)),
145  e0 = delta/(mul3 * Math::sq(d0));
146  s += RC(1, 1 + e0)/(mul * d0);
147  An = (An + lam)/4;
148  x0 = (x0 + lam)/4;
149  y0 = (y0 + lam)/4;
150  z0 = (z0 + lam)/4;
151  p0 = (p0 + lam)/4;
152  mul *= 4;
153  mul3 *= 64;
154  }
155  real
156  X = (A0 - x) / (mul * An),
157  Y = (A0 - y) / (mul * An),
158  Z = (A0 - z) / (mul * An),
159  P = -(X + Y + Z) / 2,
160  E2 = X*Y + X*Z + Y*Z - 3*P*P,
161  E3 = X*Y*Z + 2*P * (E2 + 2*P*P),
162  E4 = (2*X*Y*Z + P * (E2 + 3*P*P)) * P,
163  E5 = X*Y*Z*P*P;
164  // http://dlmf.nist.gov/19.36.E2
165  // Polynomial is
166  // (1 - 3*E2/14 + E3/6 + 9*E2^2/88 - 3*E4/22 - 9*E2*E3/52 + 3*E5/26
167  // - E2^3/16 + 3*E3^2/40 + 3*E2*E4/20 + 45*E2^2*E3/272
168  // - 9*(E3*E4+E2*E5)/68)
169  return ((471240 - 540540 * E2) * E5 +
170  (612612 * E2 - 540540 * E3 - 556920) * E4 +
171  E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
172  E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
173  (4084080 * mul * An * sqrt(An)) + 6 * s;
174  }
175 
176  Math::real EllipticFunction::RD(real x, real y, real z) {
177  // Carlson, eqs 2.28 - 2.34
178  real tolRD = pow(real(0.2) * (numeric_limits<real>::epsilon() * real(0.01)),
179  1/real(8));
180  real
181  A0 = (x + y + 3*z)/5,
182  An = A0,
183  Q = max(max(abs(A0-x), abs(A0-y)), abs(A0-z)) / tolRD,
184  x0 = x,
185  y0 = y,
186  z0 = z,
187  mul = 1,
188  s = 0;
189  while (Q >= mul * abs(An)) {
190  // Max 7 trips
191  real lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0);
192  s += 1/(mul * sqrt(z0) * (z0 + lam));
193  An = (An + lam)/4;
194  x0 = (x0 + lam)/4;
195  y0 = (y0 + lam)/4;
196  z0 = (z0 + lam)/4;
197  mul *= 4;
198  }
199  real
200  X = (A0 - x) / (mul * An),
201  Y = (A0 - y) / (mul * An),
202  Z = -(X + Y) / 3,
203  E2 = X*Y - 6*Z*Z,
204  E3 = (3*X*Y - 8*Z*Z)*Z,
205  E4 = 3 * (X*Y - Z*Z) * Z*Z,
206  E5 = X*Y*Z*Z*Z;
207  // http://dlmf.nist.gov/19.36.E2
208  // Polynomial is
209  // (1 - 3*E2/14 + E3/6 + 9*E2^2/88 - 3*E4/22 - 9*E2*E3/52 + 3*E5/26
210  // - E2^3/16 + 3*E3^2/40 + 3*E2*E4/20 + 45*E2^2*E3/272
211  // - 9*(E3*E4+E2*E5)/68)
212  return ((471240 - 540540 * E2) * E5 +
213  (612612 * E2 - 540540 * E3 - 556920) * E4 +
214  E3 * (306306 * E3 + E2 * (675675 * E2 - 706860) + 680680) +
215  E2 * ((417690 - 255255 * E2) * E2 - 875160) + 4084080) /
216  (4084080 * mul * An * sqrt(An)) + 3 * s;
217  }
218 
219  void EllipticFunction::Reset(real k2, real alpha2,
220  real kp2, real alphap2) {
221  _k2 = k2;
222  _kp2 = kp2;
223  _alpha2 = alpha2;
224  _alphap2 = alphap2;
225  _eps = _k2/Math::sq(sqrt(_kp2) + 1);
226  if (_k2) {
227  // Complete elliptic integral K(k), Carlson eq. 4.1
228  // http://dlmf.nist.gov/19.25.E1
229  _Kc = _kp2 ? RF(_kp2, 1) : Math::infinity();
230  // Complete elliptic integral E(k), Carlson eq. 4.2
231  // http://dlmf.nist.gov/19.25.E1
232  _Ec = _kp2 ? 2 * RG(_kp2, 1) : 1;
233  // D(k) = (K(k) - E(k))/k^2, Carlson eq.4.3
234  // http://dlmf.nist.gov/19.25.E1
235  _Dc = _kp2 ? RD(real(0), _kp2, 1) / 3 : Math::infinity();
236  } else {
237  _Kc = _Ec = Math::pi()/2; _Dc = _Kc/2;
238  }
239  if (_alpha2) {
240  // http://dlmf.nist.gov/19.25.E2
241  real rj = _kp2 ? RJ(0, _kp2, 1, _alphap2) : Math::infinity();
242  // Pi(alpha^2, k)
243  _Pic = _Kc + _alpha2 * rj / 3;
244  // G(alpha^2, k)
245  _Gc = _kp2 ? _Kc + (_alpha2 - _k2) * rj / 3 : RC(1, _alphap2);
246  // H(alpha^2, k)
247  _Hc = _kp2 ? _Kc - _alphap2 * rj / 3 : RC(1, _alphap2);
248  } else {
249  _Pic = _Kc; _Gc = _Ec; _Hc = _Kc - _Dc;
250  }
251  }
252 
253  /*
254  * Implementation of methods given in
255  *
256  * R. Bulirsch
257  * Numerical Calculation of Elliptic Integrals and Elliptic Functions
258  * Numericshe Mathematik 7, 78-90 (1965)
259  */
260 
261  void EllipticFunction::sncndn(real x, real& sn, real& cn, real& dn)
262  const {
263  // Bulirsch's sncndn routine, p 89.
264  real tolJAC = sqrt(numeric_limits<real>::epsilon() * real(0.01));
265  if (_kp2 != 0) {
266  real mc = _kp2, d = 0;
267  if (_kp2 < 0) {
268  d = 1 - mc;
269  mc /= -d;
270  d = sqrt(d);
271  x *= d;
272  }
273  real c = 0; // To suppress warning about uninitialized variable
274  real m[num_], n[num_];
275  unsigned l = 0;
276  for (real a = 1; l < num_ || GEOGRAPHICLIB_PANIC; ++l) {
277  // This converges quadratically. Max 5 trips
278  m[l] = a;
279  n[l] = mc = sqrt(mc);
280  c = (a + mc) / 2;
281  if (!(abs(a - mc) > tolJAC * a)) {
282  ++l;
283  break;
284  }
285  mc *= a;
286  a = c;
287  }
288  x *= c;
289  sn = sin(x);
290  cn = cos(x);
291  dn = 1;
292  if (sn != 0) {
293  real a = cn / sn;
294  c *= a;
295  while (l--) {
296  real b = m[l];
297  a *= c;
298  c *= dn;
299  dn = (n[l] + a) / (b + a);
300  a = c / b;
301  }
302  a = 1 / sqrt(c*c + 1);
303  sn = sn < 0 ? -a : a;
304  cn = c * sn;
305  if (_kp2 < 0) {
306  swap(cn, dn);
307  sn /= d;
308  }
309  }
310  } else {
311  sn = tanh(x);
312  dn = cn = 1 / cosh(x);
313  }
314  }
315 
316  Math::real EllipticFunction::F(real sn, real cn, real dn) const {
317  // Carlson, eq. 4.5 and
318  // http://dlmf.nist.gov/19.25.E5
319  real fi = abs(sn) * RF(cn*cn, dn*dn, 1);
320  // Enforce usual trig-like symmetries
321  if (cn < 0)
322  fi = 2 * K() - fi;
323  if (sn < 0)
324  fi = -fi;
325  return fi;
326  }
327 
328  Math::real EllipticFunction::E(real sn, real cn, real dn) const {
329  real
330  cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
331  ei = ( _k2 <= 0 ?
332  // Carlson, eq. 4.6 and
333  // http://dlmf.nist.gov/19.25.E9
334  RF(cn2, dn2, 1) - _k2 * sn2 * RD(cn2, dn2, 1) / 3 :
335  ( _kp2 >= 0 ?
336  // http://dlmf.nist.gov/19.25.E10
337  _kp2 * RF(cn2, dn2, 1) +
338  _k2 * _kp2 * sn2 * RD(cn2, 1, dn2) / 3 +
339  _k2 * abs(cn) / dn :
340  // http://dlmf.nist.gov/19.25.E11
341  - _kp2 * sn2 * RD(dn2, 1, cn2) / 3 + dn / abs(cn) ) );
342  ei *= abs(sn);
343  // Enforce usual trig-like symmetries
344  if (cn < 0)
345  ei = 2 * E() - ei;
346  if (sn < 0)
347  ei = -ei;
348  return ei;
349  }
350 
351  Math::real EllipticFunction::D(real sn, real cn, real dn) const {
352  // Carlson, eq. 4.8 and
353  // http://dlmf.nist.gov/19.25.E5
354  real di = abs(sn) * sn*sn * RD(cn*cn, dn*dn, 1) / 3;
355  // Enforce usual trig-like symmetries
356  if (cn < 0)
357  di = 2 * D() - di;
358  if (sn < 0)
359  di = -di;
360  return di;
361  }
362 
363  Math::real EllipticFunction::Pi(real sn, real cn, real dn) const {
364  // Carlson, eq. 4.5 and
365  // http://dlmf.nist.gov/19.25.E5
366  real
367  cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
368  pii = abs(sn) * (RF(cn2, dn2, 1) +
369  _alpha2 * sn2 * RJ(cn2, dn2, 1, 1 - _alpha2 * sn2) / 3);
370  // Enforce usual trig-like symmetries
371  if (cn < 0)
372  pii = 2 * Pi() - pii;
373  if (sn < 0)
374  pii = -pii;
375  return pii;
376  }
377 
378  Math::real EllipticFunction::G(real sn, real cn, real dn) const {
379  real
380  cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
381  gi = abs(sn) * (RF(cn2, dn2, 1) +
382  (_alpha2 - _k2) * sn2 *
383  RJ(cn2, dn2, 1, cn2 + _alphap2 * sn2) / 3);
384  // Enforce usual trig-like symmetries
385  if (cn < 0)
386  gi = 2 * G() - gi;
387  if (sn < 0)
388  gi = -gi;
389  return gi;
390  }
391 
392  Math::real EllipticFunction::H(real sn, real cn, real dn) const {
393  real
394  cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
395  hi = abs(sn) * (RF(cn2, dn2, 1) -
396  _alphap2 * sn2 *
397  RJ(cn2, dn2, 1, cn2 + _alphap2 * sn2) / 3);
398  // Enforce usual trig-like symmetries
399  if (cn < 0)
400  hi = 2 * H() - hi;
401  if (sn < 0)
402  hi = -hi;
403  return hi;
404  }
405 
406  Math::real EllipticFunction::deltaF(real sn, real cn, real dn) const {
407  // Function is periodic with period pi
408  if (cn < 0) { cn = -cn; sn = -sn; }
409  return F(sn, cn, dn) * (Math::pi()/2) / K() - atan2(sn, cn);
410  }
411 
412  Math::real EllipticFunction::deltaE(real sn, real cn, real dn) const {
413  // Function is periodic with period pi
414  if (cn < 0) { cn = -cn; sn = -sn; }
415  return E(sn, cn, dn) * (Math::pi()/2) / E() - atan2(sn, cn);
416  }
417 
418  Math::real EllipticFunction::deltaPi(real sn, real cn, real dn)
419  const {
420  // Function is periodic with period pi
421  if (cn < 0) { cn = -cn; sn = -sn; }
422  return Pi(sn, cn, dn) * (Math::pi()/2) / Pi() - atan2(sn, cn);
423  }
424 
425  Math::real EllipticFunction::deltaD(real sn, real cn, real dn) const {
426  // Function is periodic with period pi
427  if (cn < 0) { cn = -cn; sn = -sn; }
428  return D(sn, cn, dn) * (Math::pi()/2) / D() - atan2(sn, cn);
429  }
430 
431  Math::real EllipticFunction::deltaG(real sn, real cn, real dn) const {
432  // Function is periodic with period pi
433  if (cn < 0) { cn = -cn; sn = -sn; }
434  return G(sn, cn, dn) * (Math::pi()/2) / G() - atan2(sn, cn);
435  }
436 
437  Math::real EllipticFunction::deltaH(real sn, real cn, real dn) const {
438  // Function is periodic with period pi
439  if (cn < 0) { cn = -cn; sn = -sn; }
440  return H(sn, cn, dn) * (Math::pi()/2) / H() - atan2(sn, cn);
441  }
442 
443  Math::real EllipticFunction::F(real phi) const {
444  real sn = sin(phi), cn = cos(phi);
445  return (deltaF(sn, cn, Delta(sn, cn)) + phi) * K() / (Math::pi()/2);
446  }
447 
448  Math::real EllipticFunction::E(real phi) const {
449  real sn = sin(phi), cn = cos(phi);
450  return (deltaE(sn, cn, Delta(sn, cn)) + phi) * E() / (Math::pi()/2);
451  }
452 
454  real n = ceil(ang/360 - real(0.5));
455  ang -= 360 * n;
456  real
457  phi = ang * Math::degree(),
458  sn = abs(ang) == 180 ? 0 : sin(phi),
459  cn = abs(ang) == 90 ? 0 : cos(phi);
460  return E(sn, cn, Delta(sn, cn)) + 4 * E() * n;
461  }
462 
464  real sn = sin(phi), cn = cos(phi);
465  return (deltaPi(sn, cn, Delta(sn, cn)) + phi) * Pi() / (Math::pi()/2);
466  }
467 
468  Math::real EllipticFunction::D(real phi) const {
469  real sn = sin(phi), cn = cos(phi);
470  return (deltaD(sn, cn, Delta(sn, cn)) + phi) * D() / (Math::pi()/2);
471  }
472 
473  Math::real EllipticFunction::G(real phi) const {
474  real sn = sin(phi), cn = cos(phi);
475  return (deltaG(sn, cn, Delta(sn, cn)) + phi) * G() / (Math::pi()/2);
476  }
477 
478  Math::real EllipticFunction::H(real phi) const {
479  real sn = sin(phi), cn = cos(phi);
480  return (deltaH(sn, cn, Delta(sn, cn)) + phi) * H() / (Math::pi()/2);
481  }
482 
484  real tolJAC = sqrt(numeric_limits<real>::epsilon() * real(0.01));
485  real n = floor(x / (2 * _Ec) + 0.5);
486  x -= 2 * _Ec * n; // x now in [-ec, ec)
487  // Linear approximation
488  real phi = Math::pi() * x / (2 * _Ec); // phi in [-pi/2, pi/2)
489  // First order correction
490  phi -= _eps * sin(2 * phi) / 2;
491  for (int i = 0; i < num_ || GEOGRAPHICLIB_PANIC; ++i) {
492  real
493  sn = sin(phi),
494  cn = cos(phi),
495  dn = Delta(sn, cn),
496  err = (E(sn, cn, dn) - x)/dn;
497  phi = phi - err;
498  if (abs(err) < tolJAC)
499  break;
500  }
501  return n * Math::pi() + phi;
502  }
503 
504  Math::real EllipticFunction::deltaEinv(real stau, real ctau) const {
505  // Function is periodic with period pi
506  if (ctau < 0) { ctau = -ctau; stau = -stau; }
507  real tau = atan2(stau, ctau);
508  return Einv( tau * E() / (Math::pi()/2) ) - tau;
509  }
510 
511 } // namespace GeographicLib
Math::real deltaEinv(real stau, real ctau) const
static T pi()
Definition: Math.hpp:213
void Reset(real k2=0, real alpha2=0)
static T infinity()
Definition: Math.hpp:491
Math::real deltaF(real sn, real cn, real dn) const
GeographicLib::Math::real real
Definition: GeodSolve.cpp:40
Math::real deltaPi(real sn, real cn, real dn) const
static T atanh(T x)
Definition: Math.hpp:339
static real RG(real x, real y, real z)
Math::real Ed(real ang) const
Math::real Einv(real x) const
static T sq(T x)
Definition: Math.hpp:243
static real RF(real x, real y, real z)
static real RC(real x, real y)
Math::real F(real phi) const
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
void sncndn(real x, real &sn, real &cn, real &dn) const
Header for GeographicLib::EllipticFunction class.
static T degree()
Definition: Math.hpp:227
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)
Math::real deltaD(real sn, real cn, real dn) const
Math::real deltaG(real sn, real cn, real dn) const
#define GEOGRAPHICLIB_PANIC
Definition: Math.hpp:86