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