GeographicLib  1.35
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 namespace GeographicLib {
13 
14  using namespace std;
15 
16  const Math::real EllipticFunction::tol_ =
17  numeric_limits<real>::epsilon() * real(0.01);
18  const Math::real EllipticFunction::tolRF_ = pow(3 * tol_, 1/real(8));
19  const Math::real EllipticFunction::tolRD_ = pow(real(0.2) * tol_, 1/real(8));
20  const Math::real EllipticFunction::tolRG0_ = real(2.7) * sqrt(tol_);
21  const Math::real EllipticFunction::tolJAC_ = sqrt(tol_);
22 
23  /*
24  * Implementation of methods given in
25  *
26  * B. C. Carlson
27  * Computation of elliptic integrals
28  * Numerical Algorithms 10, 13-26 (1995)
29  */
30 
31  Math::real EllipticFunction::RF(real x, real y, real z) throw() {
32  // Carlson, eqs 2.2 - 2.7
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) throw() {
67  // Carlson, eqs 2.36 - 2.38
68  real xn = sqrt(x), yn = sqrt(y);
69  if (xn < yn) swap(xn, yn);
70  while (abs(xn-yn) > tolRG0_ * xn) {
71  // Max 4 trips
72  real t = (xn + yn) /2;
73  yn = sqrt(xn * yn);
74  xn = t;
75  }
76  return Math::pi<real>() / (xn + yn);
77  }
78 
79  Math::real EllipticFunction::RC(real x, real y) throw() {
80  return ( !(x >= y) ? // x < y and catch nans
81  // http://dlmf.nist.gov/19.2.E18
82  atan(sqrt((y - x) / x)) / sqrt(y - x) :
83  ( x == y && y > 0 ? 1 / sqrt(y) :
84  Math::atanh( y > 0 ?
85  // http://dlmf.nist.gov/19.2.E19
86  sqrt((x - y) / x) :
87  // http://dlmf.nist.gov/19.2.E20
88  sqrt(x / (x - y)) ) / sqrt(x - y) ) );
89  }
90 
91  Math::real EllipticFunction::RG(real x, real y, real z) throw() {
92  if (z == 0)
93  swap(y, z);
94  // Carlson, eq 1.7
95  return (z * RF(x, y, z) - (x-z) * (y-z) * RD(x, y, z) / 3
96  + sqrt(x * y / z)) / 2;
97  }
98 
99  Math::real EllipticFunction::RG(real x, real y) throw() {
100  // Carlson, eqs 2.36 - 2.39
101  real
102  x0 = sqrt(max(x, y)),
103  y0 = sqrt(min(x, y)),
104  xn = x0,
105  yn = y0,
106  s = 0,
107  mul = real(0.25);
108  while (abs(xn-yn) > tolRG0_ * xn) {
109  // Max 4 trips
110  real t = (xn + yn) /2;
111  yn = sqrt(xn * yn);
112  xn = t;
113  mul *= 2;
114  t = xn - yn;
115  s += mul * t * t;
116  }
117  return (Math::sq( (x0 + y0)/2 ) - s) * Math::pi<real>() / (2 * (xn + yn));
118  }
119 
120  Math::real EllipticFunction::RJ(real x, real y, real z, real p) throw() {
121  // Carlson, eqs 2.17 - 2.25
122  real
123  A0 = (x + y + z + 2*p)/5,
124  An = A0,
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_,
127  x0 = x,
128  y0 = y,
129  z0 = z,
130  p0 = p,
131  mul = 1,
132  mul3 = 1,
133  s = 0;
134  while (Q >= mul * abs(An)) {
135  // Max 7 trips
136  real
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)),
139  e0 = delta/(mul3 * Math::sq(d0));
140  s += RC(1, 1 + e0)/(mul * d0);
141  An = (An + lam)/4;
142  x0 = (x0 + lam)/4;
143  y0 = (y0 + lam)/4;
144  z0 = (z0 + lam)/4;
145  p0 = (p0 + lam)/4;
146  mul *= 4;
147  mul3 *= 64;
148  }
149  real
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,
157  E5 = X*Y*Z*P*P;
158  // http://dlmf.nist.gov/19.36.E2
159  // Polynomial is
160  // (1 - 3*E2/14 + E3/6 + 9*E2^2/88 - 3*E4/22 - 9*E2*E3/52 + 3*E5/26
161  // - E2^3/16 + 3*E3^2/40 + 3*E2*E4/20 + 45*E2^2*E3/272
162  // - 9*(E3*E4+E2*E5)/68)
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;
168  }
169 
170  Math::real EllipticFunction::RD(real x, real y, real z) throw() {
171  // Carlson, eqs 2.28 - 2.34
172  real
173  A0 = (x + y + 3*z)/5,
174  An = A0,
175  Q = max(max(abs(A0-x), abs(A0-y)), abs(A0-z)) / tolRD_,
176  x0 = x,
177  y0 = y,
178  z0 = z,
179  mul = 1,
180  s = 0;
181  while (Q >= mul * abs(An)) {
182  // Max 7 trips
183  real lam = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0);
184  s += 1/(mul * sqrt(z0) * (z0 + lam));
185  An = (An + lam)/4;
186  x0 = (x0 + lam)/4;
187  y0 = (y0 + lam)/4;
188  z0 = (z0 + lam)/4;
189  mul *= 4;
190  }
191  real
192  X = (A0 - x) / (mul * An),
193  Y = (A0 - y) / (mul * An),
194  Z = -(X + Y) / 3,
195  E2 = X*Y - 6*Z*Z,
196  E3 = (3*X*Y - 8*Z*Z)*Z,
197  E4 = 3 * (X*Y - Z*Z) * Z*Z,
198  E5 = X*Y*Z*Z*Z;
199  // http://dlmf.nist.gov/19.36.E2
200  // Polynomial is
201  // (1 - 3*E2/14 + E3/6 + 9*E2^2/88 - 3*E4/22 - 9*E2*E3/52 + 3*E5/26
202  // - E2^3/16 + 3*E3^2/40 + 3*E2*E4/20 + 45*E2^2*E3/272
203  // - 9*(E3*E4+E2*E5)/68)
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;
209  }
210 
211  EllipticFunction::EllipticFunction(real k2, real alpha2) throw()
212  : _k2(k2)
213  , _kp2(1 - k2)
214  , _alpha2(alpha2)
215  , _alphap2(1 - alpha2)
216  , _eps(_k2/Math::sq(sqrt(_kp2) + 1))
217  // Don't initialize _Kc, _Ec, _Dc since this constructor might be called
218  // before the static real constants tolRF_, etc., are initialized.
219  , _init(false)
220  {}
221 
223  real kp2, real alphap2) throw()
224  : _k2(k2)
225  , _kp2(kp2)
226  , _alpha2(alpha2)
227  , _alphap2(alphap2)
228  , _eps(_k2/Math::sq(sqrt(_kp2) + 1))
229  , _init(false)
230  {}
231 
232  void EllipticFunction::Reset(real k2, real alpha2,
233  real kp2, real alphap2) throw() {
234  _k2 = k2;
235  _kp2 = kp2;
236  _alpha2 = alpha2;
237  _alphap2 = alphap2;
238  _eps = _k2/Math::sq(sqrt(_kp2) + 1);
239  _init = false;
240  }
241 
242  bool EllipticFunction::Init() const throw() {
243  // Complete elliptic integral K(k), Carlson eq. 4.1
244  // http://dlmf.nist.gov/19.25.E1
245  _Kc = _kp2 ? RF(_kp2, 1) : Math::infinity<real>();
246  // Complete elliptic integral E(k), Carlson eq. 4.2
247  // http://dlmf.nist.gov/19.25.E1
248  _Ec = _kp2 ? 2 * RG(_kp2, 1) : 1;
249  // D(k) = (K(k) - E(k))/m, Carlson eq.4.3
250  // http://dlmf.nist.gov/19.25.E1
251  _Dc = _kp2 ? RD(real(0), _kp2, 1) / 3 : Math::infinity<real>();
252  if (_alpha2) {
253  // http://dlmf.nist.gov/19.25.E2
254  real rj = _kp2 ? RJ(0, _kp2, 1, _alphap2) : Math::infinity<real>();
255  // Pi(alpha^2, k)
256  _Pic = _Kc + _alpha2 * rj / 3;
257  // G(alpha^2, k)
258  _Gc = _kp2 ? _Kc + (_alpha2 - _k2) * rj / 3 : RC(1, _alphap2);
259  // H(alpha^2, k)
260  _Hc = _kp2 ? _Kc - _alphap2 * rj / 3 : RC(1, _alphap2);
261  } else {
262  _Pic = _Kc; _Gc = _Ec; _Hc = _Kc - _Dc;
263  }
264  return _init = true;
265  }
266 
267  /*
268  * Implementation of methods given in
269  *
270  * R. Bulirsch
271  * Numerical Calculation of Elliptic Integrals and Elliptic Functions
272  * Numericshe Mathematik 7, 78-90 (1965)
273  */
274 
275  void EllipticFunction::sncndn(real x, real& sn, real& cn, real& dn)
276  const throw() {
277  // Bulirsch's sncndn routine, p 89.
278  if (_kp2 != 0) {
279  real mc = _kp2, d = 0;
280  if (_kp2 < 0) {
281  d = 1 - mc;
282  mc /= -d;
283  d = sqrt(d);
284  x *= d;
285  }
286  real c = 0; // To suppress warning about uninitialized variable
287  real m[num_], n[num_];
288  unsigned l = 0;
289  for (real a = 1; l < num_; ++l) {
290  // Max 5 trips
291  m[l] = a;
292  n[l] = mc = sqrt(mc);
293  c = (a + mc) / 2;
294  if (!(abs(a - mc) > tolJAC_ * a)) {
295  ++l;
296  break;
297  }
298  mc *= a;
299  a = c;
300  }
301  x *= c;
302  sn = sin(x);
303  cn = cos(x);
304  dn = 1;
305  if (sn != 0) {
306  real a = cn / sn;
307  c *= a;
308  while (l--) {
309  real b = m[l];
310  a *= c;
311  c *= dn;
312  dn = (n[l] + a) / (b + a);
313  a = c / b;
314  }
315  a = 1 / sqrt(c*c + 1);
316  sn = sn < 0 ? -a : a;
317  cn = c * sn;
318  if (_kp2 < 0) {
319  swap(cn, dn);
320  sn /= d;
321  }
322  }
323  } else {
324  sn = tanh(x);
325  dn = cn = 1 / cosh(x);
326  }
327  }
328 
329  Math::real EllipticFunction::F(real sn, real cn, real dn) const throw() {
330  // Carlson, eq. 4.5 and
331  // http://dlmf.nist.gov/19.25.E5
332  real fi = abs(sn) * RF(cn*cn, dn*dn, 1);
333  // Enforce usual trig-like symmetries
334  if (cn < 0)
335  fi = 2 * K() - fi;
336  if (sn < 0)
337  fi = -fi;
338  return fi;
339  }
340 
341  Math::real EllipticFunction::E(real sn, real cn, real dn) const throw() {
342  real
343  cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
344  ei = ( _k2 <= 0 ?
345  // Carlson, eq. 4.6 and
346  // http://dlmf.nist.gov/19.25.E9
347  RF(cn2, dn2, 1) - _k2 * sn2 * RD(cn2, dn2, 1) / 3 :
348  ( _kp2 >= 0 ?
349  // http://dlmf.nist.gov/19.25.E10
350  _kp2 * RF(cn2, dn2, 1) +
351  _k2 * _kp2 * sn2 * RD(cn2, 1, dn2) / 3 +
352  _k2 * abs(cn) / dn :
353  // http://dlmf.nist.gov/19.25.E11
354  - _kp2 * sn2 * RD(dn2, 1, cn2) / 3 + dn / abs(cn) ) );
355  ei *= abs(sn);
356  // Enforce usual trig-like symmetries
357  if (cn < 0)
358  ei = 2 * E() - ei;
359  if (sn < 0)
360  ei = -ei;
361  return ei;
362  }
363 
364  Math::real EllipticFunction::D(real sn, real cn, real dn) const throw() {
365  // Carlson, eq. 4.8 and
366  // http://dlmf.nist.gov/19.25.E5
367  real di = abs(sn) * sn*sn * RD(cn*cn, dn*dn, 1) / 3;
368  // Enforce usual trig-like symmetries
369  if (cn < 0)
370  di = 2 * D() - di;
371  if (sn < 0)
372  di = -di;
373  return di;
374  }
375 
376  Math::real EllipticFunction::Pi(real sn, real cn, real dn) const throw() {
377  // Carlson, eq. 4.5 and
378  // http://dlmf.nist.gov/19.25.E5
379  real
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);
383  // Enforce usual trig-like symmetries
384  if (cn < 0)
385  pii = 2 * Pi() - pii;
386  if (sn < 0)
387  pii = -pii;
388  return pii;
389  }
390 
391  Math::real EllipticFunction::G(real sn, real cn, real dn) const throw() {
392  real
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);
397  // Enforce usual trig-like symmetries
398  if (cn < 0)
399  gi = 2 * G() - gi;
400  if (sn < 0)
401  gi = -gi;
402  return gi;
403  }
404 
405  Math::real EllipticFunction::H(real sn, real cn, real dn) const throw() {
406  real
407  cn2 = cn*cn, dn2 = dn*dn, sn2 = sn*sn,
408  hi = abs(sn) * (RF(cn2, dn2, 1) -
409  _alphap2 * sn2 *
410  RJ(cn2, dn2, 1, cn2 + _alphap2 * sn2) / 3);
411  // Enforce usual trig-like symmetries
412  if (cn < 0)
413  hi = 2 * H() - hi;
414  if (sn < 0)
415  hi = -hi;
416  return hi;
417  }
418 
419  Math::real EllipticFunction::deltaF(real sn, real cn, real dn) const throw() {
420  // Function is periodic with period pi
421  if (cn < 0) { cn = -cn; sn = -sn; }
422  return F(sn, cn, dn) * (Math::pi<real>()/2) / K() - atan2(sn, cn);
423  }
424 
425  Math::real EllipticFunction::deltaE(real sn, real cn, real dn) const throw() {
426  // Function is periodic with period pi
427  if (cn < 0) { cn = -cn; sn = -sn; }
428  return E(sn, cn, dn) * (Math::pi<real>()/2) / E() - atan2(sn, cn);
429  }
430 
431  Math::real EllipticFunction::deltaPi(real sn, real cn, real dn)
432  const throw() {
433  // Function is periodic with period pi
434  if (cn < 0) { cn = -cn; sn = -sn; }
435  return Pi(sn, cn, dn) * (Math::pi<real>()/2) / Pi() - atan2(sn, cn);
436  }
437 
438  Math::real EllipticFunction::deltaD(real sn, real cn, real dn) const throw() {
439  // Function is periodic with period pi
440  if (cn < 0) { cn = -cn; sn = -sn; }
441  return D(sn, cn, dn) * (Math::pi<real>()/2) / D() - atan2(sn, cn);
442  }
443 
444  Math::real EllipticFunction::deltaG(real sn, real cn, real dn) const throw() {
445  // Function is periodic with period pi
446  if (cn < 0) { cn = -cn; sn = -sn; }
447  return G(sn, cn, dn) * (Math::pi<real>()/2) / G() - atan2(sn, cn);
448  }
449 
450  Math::real EllipticFunction::deltaH(real sn, real cn, real dn) const throw() {
451  // Function is periodic with period pi
452  if (cn < 0) { cn = -cn; sn = -sn; }
453  return H(sn, cn, dn) * (Math::pi<real>()/2) / H() - atan2(sn, cn);
454  }
455 
456  Math::real EllipticFunction::F(real phi) const throw() {
457  real sn = sin(phi), cn = cos(phi);
458  return (deltaF(sn, cn, Delta(sn, cn)) + phi) * K() / (Math::pi<real>()/2);
459  }
460 
461  Math::real EllipticFunction::E(real phi) const throw() {
462  real sn = sin(phi), cn = cos(phi);
463  return (deltaE(sn, cn, Delta(sn, cn)) + phi) * E() / (Math::pi<real>()/2);
464  }
465 
466  Math::real EllipticFunction::Ed(real ang) const throw() {
467  real n = ceil(ang/360 - real(0.5));
468  ang -= 360 * n;
469  real
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;
474  }
475 
476  Math::real EllipticFunction::Pi(real phi) const throw() {
477  real sn = sin(phi), cn = cos(phi);
478  return (deltaPi(sn, cn, Delta(sn, cn)) + phi) * Pi() / (Math::pi<real>()/2);
479  }
480 
481  Math::real EllipticFunction::D(real phi) const throw() {
482  real sn = sin(phi), cn = cos(phi);
483  return (deltaD(sn, cn, Delta(sn, cn)) + phi) * D() / (Math::pi<real>()/2);
484  }
485 
486  Math::real EllipticFunction::G(real phi) const throw() {
487  real sn = sin(phi), cn = cos(phi);
488  return (deltaG(sn, cn, Delta(sn, cn)) + phi) * G() / (Math::pi<real>()/2);
489  }
490 
491  Math::real EllipticFunction::H(real phi) const throw() {
492  real sn = sin(phi), cn = cos(phi);
493  return (deltaH(sn, cn, Delta(sn, cn)) + phi) * H() / (Math::pi<real>()/2);
494  }
495 
496  Math::real EllipticFunction::Einv(real x) const throw() {
497  _init || Init();
498  real n = floor(x / (2 * _Ec) + 0.5);
499  x -= 2 * _Ec * n; // x now in [-ec, ec)
500  // Linear approximation
501  real phi = Math::pi<real>() * x / (2 * _Ec); // phi in [-pi/2, pi/2)
502  // First order correction
503  phi -= _eps * sin(2 * phi) / 2;
504  for (int i = 0; i < num_; ++i) {
505  real
506  sn = sin(phi),
507  cn = cos(phi),
508  dn = Delta(sn, cn),
509  err = (E(sn, cn, dn) - x)/dn;
510  phi = phi - err;
511  if (abs(err) < tolJAC_)
512  break;
513  }
514  return n * Math::pi<real>() + phi;
515  }
516 
517  Math::real EllipticFunction::deltaEinv(real stau, real ctau) const throw() {
518  // Function is periodic with period pi
519  if (ctau < 0) { ctau = -ctau; stau = -stau; }
520  real tau = atan2(stau, ctau);
521  return Einv( tau * E() / (Math::pi<real>()/2) ) - tau;
522  }
523 
524 } // namespace GeographicLib
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
Definition: GeodSolve.cpp:40
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:73
Math::real deltaPi(real sn, real cn, real dn) const
static T atanh(T x)
Definition: Math.hpp:315
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:153
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