GeographicLib  1.35
EllipticFunction.hpp
Go to the documentation of this file.
1 /**
2  * \file EllipticFunction.hpp
3  * \brief Header 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 
10 #if !defined(GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP)
11 #define GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP 1
12 
14 
15 namespace GeographicLib {
16 
17  /**
18  * \brief Elliptic integrals and functions
19  *
20  * This provides the elliptic functions and integrals needed for Ellipsoid,
21  * GeodesicExact, and TransverseMercatorExact. Two categories of function
22  * are provided:
23  * - \e static functions to compute symmetric elliptic integrals
24  * (http://dlmf.nist.gov/19.16.i)
25  * - \e member functions to compute Legrendre's elliptic
26  * integrals (http://dlmf.nist.gov/19.2.ii) and the
27  * Jacobi elliptic functions (http://dlmf.nist.gov/22.2).
28  * .
29  * In the latter case, an object is constructed giving the modulus \e k (and
30  * optionally the parameter &alpha;<sup>2</sup>). The modulus is always
31  * passed as its square <i>k</i><sup>2</sup> which allows \e k to be pure
32  * imaginary (<i>k</i><sup>2</sup> &lt; 0). (Confusingly, Abramowitz and
33  * Stegun call \e m = <i>k</i><sup>2</sup> the "parameter" and \e n =
34  * &alpha;<sup>2</sup> the "characteristic".)
35  *
36  * In geodesic applications, it is convenient to separate the incomplete
37  * integrals into secular and periodic components, e.g.,
38  * \f[
39  * E(\phi, k) = (2 E(\phi) / \pi) [ \phi + \delta E(\phi, k) ]
40  * \f]
41  * where &delta;\e E(&phi;, \e k) is an odd periodic function with period
42  * &pi;.
43  *
44  * The computation of the elliptic integrals uses the algorithms given in
45  * - B. C. Carlson,
46  * <a href="http://dx.doi.org/10.1007/BF02198293"> Computation of real or
47  * complex elliptic integrals</a>, Numerical Algorithms 10, 13--26 (1995)
48  * .
49  * with the additional optimizations given in http://dlmf.nist.gov/19.36.i.
50  * The computation of the Jacobi elliptic functions uses the algorithm given
51  * in
52  * - R. Bulirsch,
53  * <a href="http://dx.doi.org/10.1007/BF01397975"> Numerical Calculation of
54  * Elliptic Integrals and Elliptic Functions</a>, Numericshe Mathematik 7,
55  * 78--90 (1965).
56  * .
57  * The notation follows http://dlmf.nist.gov/19 and http://dlmf.nist.gov/22
58  *
59  * Example of use:
60  * \include example-EllipticFunction.cpp
61  **********************************************************************/
63  private:
64  typedef Math::real real;
65  static const real tol_;
66  static const real tolRF_;
67  static const real tolRD_;
68  static const real tolRG0_;
69  static const real tolJAC_;
70  enum { num_ = 13 }; // Max depth required for sncndn. Probably 5 is enough.
71  real _k2, _kp2, _alpha2, _alphap2, _eps;
72  mutable bool _init;
73  mutable real _Kc, _Ec, _Dc, _Pic, _Gc, _Hc;
74  bool Init() const throw();
75  public:
76  /** \name Constructor
77  **********************************************************************/
78  ///@{
79  /**
80  * Constructor specifying the modulus and parameter.
81  *
82  * @param[in] k2 the square of the modulus <i>k</i><sup>2</sup>.
83  * <i>k</i><sup>2</sup> must lie in (-&infin;, 1). (No checking is
84  * done.)
85  * @param[in] alpha2 the parameter &alpha;<sup>2</sup>.
86  * &alpha;<sup>2</sup> must lie in (-&infin;, 1). (No checking is done.)
87  *
88  * If only elliptic integrals of the first and second kinds are needed,
89  * then set &alpha;<sup>2</sup> = 0 (the default value); in this case, we
90  * have &Pi;(&phi;, 0, \e k) = \e F(&phi;, \e k), \e G(&phi;, 0, \e k) = \e
91  * E(&phi;, \e k), and \e H(&phi;, 0, \e k) = \e F(&phi;, \e k) - \e
92  * D(&phi;, \e k).
93  **********************************************************************/
94  EllipticFunction(real k2 = 0, real alpha2 = 0) throw();
95 
96  /**
97  * Constructor specifying the modulus and parameter and their complements.
98  *
99  * @param[in] k2 the square of the modulus <i>k</i><sup>2</sup>.
100  * <i>k</i><sup>2</sup> must lie in (-&infin;, 1). (No checking is
101  * done.)
102  * @param[in] alpha2 the parameter &alpha;<sup>2</sup>.
103  * &alpha;<sup>2</sup> must lie in (-&infin;, 1). (No checking is done.)
104  * @param[in] kp2 the complementary modulus squared <i>k'</i><sup>2</sup> =
105  * 1 &minus; <i>k</i><sup>2</sup>.
106  * @param[in] alphap2 the complementary parameter &alpha;'<sup>2</sup> = 1
107  * &minus; &alpha;<sup>2</sup>.
108  *
109  * The arguments must satisfy \e k2 + \e kp2 = 1 and \e alpha2 + \e alphap2
110  * = 1. (No checking is done that these conditions are met.) This
111  * constructor is provided to enable accuracy to be maintained, e.g., when
112  * \e k is very close to unity.
113  **********************************************************************/
114  EllipticFunction(real k2, real alpha2, real kp2, real alphap2) throw();
115 
116  /**
117  * Reset the modulus and parameter.
118  *
119  * @param[in] k2 the new value of square of the modulus
120  * <i>k</i><sup>2</sup> which must lie in (-&infin;, 1). (No checking is
121  * done.)
122  * @param[in] alpha2 the new value of parameter &alpha;<sup>2</sup>.
123  * &alpha;<sup>2</sup> must lie in (-&infin;, 1). (No checking is done.)
124  **********************************************************************/
125  void Reset(real k2 = 0, real alpha2 = 0) throw()
126  { Reset(k2, alpha2, 1 - k2, 1 - alpha2); }
127 
128  /**
129  * Reset the modulus and parameter supplying also their complements.
130  *
131  * @param[in] k2 the square of the modulus <i>k</i><sup>2</sup>.
132  * <i>k</i><sup>2</sup> must lie in (-&infin;, 1). (No checking is
133  * done.)
134  * @param[in] alpha2 the parameter &alpha;<sup>2</sup>.
135  * &alpha;<sup>2</sup> must lie in (-&infin;, 1). (No checking is done.)
136  * @param[in] kp2 the complementary modulus squared <i>k'</i><sup>2</sup> =
137  * 1 &minus; <i>k</i><sup>2</sup>.
138  * @param[in] alphap2 the complementary parameter &alpha;'<sup>2</sup> = 1
139  * &minus; &alpha;<sup>2</sup>.
140  *
141  * The arguments must satisfy \e k2 + \e kp2 = 1 and \e alpha2 + \e alphap2
142  * = 1. (No checking is done that these conditions are met.) This
143  * constructor is provided to enable accuracy to be maintained, e.g., when
144  * is very small.
145  **********************************************************************/
146  void Reset(real k2, real alpha2, real kp2, real alphap2) throw();
147 
148  ///@}
149 
150  /** \name Inspector functions.
151  **********************************************************************/
152  ///@{
153  /**
154  * @return the square of the modulus <i>k</i><sup>2</sup>.
155  **********************************************************************/
156  Math::real k2() const throw() { return _k2; }
157 
158  /**
159  * @return the square of the complementary modulus <i>k'</i><sup>2</sup> =
160  * 1 &minus; <i>k</i><sup>2</sup>.
161  **********************************************************************/
162  Math::real kp2() const throw() { return _kp2; }
163 
164  /**
165  * @return the parameter &alpha;<sup>2</sup>.
166  **********************************************************************/
167  Math::real alpha2() const throw() { return _alpha2; }
168 
169  /**
170  * @return the complementary parameter &alpha;'<sup>2</sup> = 1 &minus;
171  * &alpha;<sup>2</sup>.
172  **********************************************************************/
173  Math::real alphap2() const throw() { return _alphap2; }
174  ///@}
175 
176  /// \cond SKIP
177  /**
178  * @return the square of the modulus <i>k</i><sup>2</sup>.
179  *
180  * <b>DEPRECATED</b>, use k2() instead.
181  **********************************************************************/
182  Math::real m() const throw() { return _k2; }
183  /**
184  * @return the square of the complementary modulus <i>k'</i><sup>2</sup> =
185  * 1 &minus; <i>k</i><sup>2</sup>.
186  *
187  * <b>DEPRECATED</b>, use kp2() instead.
188  **********************************************************************/
189  Math::real m1() const throw() { return _kp2; }
190  /// \endcond
191 
192  /** \name Complete elliptic integrals.
193  **********************************************************************/
194  ///@{
195  /**
196  * The complete integral of the first kind.
197  *
198  * @return \e K(\e k).
199  *
200  * \e K(\e k) is defined in http://dlmf.nist.gov/19.2.E4
201  * \f[
202  * K(k) = \int_0^{\pi/2} \frac1{\sqrt{1-k^2\sin^2\phi}}\,d\phi.
203  * \f]
204  **********************************************************************/
205  Math::real K() const throw() { _init || Init(); return _Kc; }
206 
207  /**
208  * The complete integral of the second kind.
209  *
210  * @return \e E(\e k)
211  *
212  * \e E(\e k) is defined in http://dlmf.nist.gov/19.2.E5
213  * \f[
214  * E(k) = \int_0^{\pi/2} \sqrt{1-k^2\sin^2\phi}\,d\phi.
215  * \f]
216  **********************************************************************/
217  Math::real E() const throw() { _init || Init(); return _Ec; }
218 
219  /**
220  * Jahnke's complete integral.
221  *
222  * @return \e D(\e k).
223  *
224  * \e D(\e k) is defined in http://dlmf.nist.gov/19.2.E6
225  * \f[
226  * D(k) = \int_0^{\pi/2} \frac{\sin^2\phi}{\sqrt{1-k^2\sin^2\phi}}\,d\phi.
227  * \f]
228  **********************************************************************/
229  Math::real D() const throw() { _init || Init(); return _Dc; }
230 
231  /**
232  * The difference between the complete integrals of the first and second
233  * kinds.
234  *
235  * @return \e K(\e k) &minus; \e E(\e k).
236  **********************************************************************/
237  Math::real KE() const throw() { _init || Init(); return _k2 * _Dc; }
238 
239  /**
240  * The complete integral of the third kind.
241  *
242  * @return &Pi;(&alpha;<sup>2</sup>, \e k)
243  *
244  * &Pi;(&alpha;<sup>2</sup>, \e k) is defined in
245  * http://dlmf.nist.gov/19.2.E7
246  * \f[
247  * \Pi(\alpha^2, k) = \int_0^{\pi/2}
248  * \frac1{\sqrt{1-k^2\sin^2\phi}(1 - \alpha^2\sin^2\phi_)}\,d\phi.
249  * \f]
250  **********************************************************************/
251  Math::real Pi() const throw() { _init || Init(); return _Pic; }
252 
253  /**
254  * Legendre's complete geodesic longitude integral.
255  *
256  * @return \e G(&alpha;<sup>2</sup>, \e k)
257  *
258  * \e G(&alpha;<sup>2</sup>, \e k) is given by
259  * \f[
260  * G(\alpha^2, k) = \int_0^{\pi/2}
261  * \frac{\sqrt{1-k^2\sin^2\phi}}{1 - \alpha^2\sin^2\phi}\,d\phi.
262  * \f]
263  **********************************************************************/
264  Math::real G() const throw() { _init || Init(); return _Gc; }
265 
266  /**
267  * Cayley's complete geodesic longitude difference integral.
268  *
269  * @return \e H(&alpha;<sup>2</sup>, \e k)
270  *
271  * \e H(&alpha;<sup>2</sup>, \e k) is given by
272  * \f[
273  * H(\alpha^2, k) = \int_0^{\pi/2}
274  * \frac{\cos^2\phi}{(1-\alpha^2\sin^2\phi)\sqrt{1-k^2\sin^2\phi}}
275  * \,d\phi.
276  * \f]
277  **********************************************************************/
278  Math::real H() const throw() { _init || Init(); return _Hc; }
279  ///@}
280 
281  /** \name Incomplete elliptic integrals.
282  **********************************************************************/
283  ///@{
284  /**
285  * The incomplete integral of the first kind.
286  *
287  * @param[in] phi
288  * @return \e F(&phi;, \e k).
289  *
290  * \e F(&phi;, \e k) is defined in http://dlmf.nist.gov/19.2.E4
291  * \f[
292  * F(\phi, k) = \int_0^\phi \frac1{\sqrt{1-k^2\sin^2\theta}}\,d\theta.
293  * \f]
294  **********************************************************************/
295  Math::real F(real phi) const throw();
296 
297  /**
298  * The incomplete integral of the second kind.
299  *
300  * @param[in] phi
301  * @return \e E(&phi;, \e k).
302  *
303  * \e E(&phi;, \e k) is defined in http://dlmf.nist.gov/19.2.E5
304  * \f[
305  * E(\phi, k) = \int_0^\phi \sqrt{1-k^2\sin^2\theta}\,d\theta.
306  * \f]
307  **********************************************************************/
308  Math::real E(real phi) const throw();
309 
310  /**
311  * The incomplete integral of the second kind with the argument given in
312  * degrees.
313  *
314  * @param[in] ang in <i>degrees</i>.
315  * @return \e E(&pi; <i>ang</i>/180, \e k).
316  **********************************************************************/
317  Math::real Ed(real ang) const throw();
318 
319  /**
320  * The inverse of the incomplete integral of the second kind.
321  *
322  * @param[in] x
323  * @return &phi; = <i>E</i><sup>&minus;1</sup>(\e x, \e k); i.e., the
324  * solution of such that \e E(&phi;, \e k) = \e x.
325  **********************************************************************/
326  Math::real Einv(real x) const throw();
327 
328  /**
329  * The incomplete integral of the third kind.
330  *
331  * @param[in] phi
332  * @return &Pi;(&phi;, &alpha;<sup>2</sup>, \e k).
333  *
334  * &Pi;(&phi;, &alpha;<sup>2</sup>, \e k) is defined in
335  * http://dlmf.nist.gov/19.2.E7
336  * \f[
337  * \Pi(\phi, \alpha^2, k) = \int_0^\phi
338  * \frac1{\sqrt{1-k^2\sin^2\theta}(1 - \alpha^2\sin^2\theta_)}\,d\theta.
339  * \f]
340  **********************************************************************/
341  Math::real Pi(real phi) const throw();
342 
343  /**
344  * Jahnke's incomplete elliptic integral.
345  *
346  * @param[in] phi
347  * @return \e D(&phi;, \e k).
348  *
349  * \e D(&phi;, \e k) is defined in http://dlmf.nist.gov/19.2.E4
350  * \f[
351  * D(\phi, k) = \int_0^\phi
352  * \frac{\sin^2\theta}{\sqrt{1-k^2\sin^2\theta}}\,d\theta.
353  * \f]
354  **********************************************************************/
355  Math::real D(real phi) const throw();
356 
357  /**
358  * Legendre's geodesic longitude integral.
359  *
360  * @param[in] phi
361  * @return \e G(&phi;, &alpha;<sup>2</sup>, \e k).
362  *
363  * \e G(&phi;, &alpha;<sup>2</sup>, \e k) is defined by
364  * \f[
365  * \begin{aligned}
366  * G(\phi, \alpha^2, k) &=
367  * \frac{k^2}{\alpha^2} F(\phi, k) +
368  * \biggl(1 - \frac{k^2}{\alpha^2}\biggr) \Pi(\phi, \alpha^2, k) \\
369  * &= \int_0^\phi
370  * \frac{\sqrt{1-k^2\sin^2\theta}}{1 - \alpha^2\sin^2\theta}\,d\theta.
371  * \end{aligned}
372  * \f]
373  *
374  * Legendre expresses the longitude of a point on the geodesic in terms of
375  * this combination of elliptic integrals in Exercices de Calcul
376  * Int&eacute;gral, Vol. 1 (1811), p. 181,
377  * http://books.google.com/books?id=riIOAAAAQAAJ&pg=PA181.
378  *
379  * See \ref geodellip for the expression for the longitude in terms of this
380  * function.
381  **********************************************************************/
382  Math::real G(real phi) const throw();
383 
384  /**
385  * Cayley's geodesic longitude difference integral.
386  *
387  * @param[in] phi
388  * @return \e H(&phi;, &alpha;<sup>2</sup>, \e k).
389  *
390  * \e H(&phi;, &alpha;<sup>2</sup>, \e k) is defined by
391  * \f[
392  * \begin{aligned}
393  * H(\phi, \alpha^2, k) &=
394  * \frac1{\alpha^2} F(\phi, k) +
395  * \biggl(1 - \frac1{\alpha^2}\biggr) \Pi(\phi, \alpha^2, k) \\
396  * &= \int_0^\phi
397  * \frac{\cos^2\theta}{(1-\alpha^2\sin^2\theta)\sqrt{1-k^2\sin^2\theta}}
398  * \,d\theta.
399  * \end{aligned}
400  * \f]
401  *
402  * Cayley expresses the longitude difference of a point on the geodesic in
403  * terms of this combination of elliptic integrals in Phil. Mag. <b>40</b>
404  * (1870), p. 333, http://books.google.com/books?id=Zk0wAAAAIAAJ&pg=PA333.
405  *
406  * See \ref geodellip for the expression for the longitude in terms of this
407  * function.
408  **********************************************************************/
409  Math::real H(real phi) const throw();
410  ///@}
411 
412  /** \name Incomplete integrals in terms of Jacobi elliptic functions.
413  **********************************************************************/
414  /**
415  * The incomplete integral of the first kind in terms of Jacobi elliptic
416  * functions.
417  *
418  * @param[in] sn = sin&phi;
419  * @param[in] cn = cos&phi;
420  * @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
421  * sin<sup>2</sup>&phi;)
422  * @return \e F(&phi;, \e k) as though &phi; &isin; (&minus;&pi;, &pi;].
423  **********************************************************************/
424  Math::real F(real sn, real cn, real dn) const throw();
425 
426  /**
427  * The incomplete integral of the second kind in terms of Jacobi elliptic
428  * functions.
429  *
430  * @param[in] sn = sin&phi;
431  * @param[in] cn = cos&phi;
432  * @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
433  * sin<sup>2</sup>&phi;)
434  * @return \e E(&phi;, \e k) as though &phi; &isin; (&minus;&pi;, &pi;].
435  **********************************************************************/
436  Math::real E(real sn, real cn, real dn) const throw();
437 
438  /**
439  * The incomplete integral of the third kind in terms of Jacobi elliptic
440  * functions.
441  *
442  * @param[in] sn = sin&phi;
443  * @param[in] cn = cos&phi;
444  * @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
445  * sin<sup>2</sup>&phi;)
446  * @return &Pi;(&phi;, &alpha;<sup>2</sup>, \e k) as though &phi; &isin;
447  * (&minus;&pi;, &pi;].
448  **********************************************************************/
449  Math::real Pi(real sn, real cn, real dn) const throw();
450 
451  /**
452  * Jahnke's incomplete elliptic integral in terms of Jacobi elliptic
453  * functions.
454  *
455  * @param[in] sn = sin&phi;
456  * @param[in] cn = cos&phi;
457  * @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
458  * sin<sup>2</sup>&phi;)
459  * @return \e D(&phi;, \e k) as though &phi; &isin; (&minus;&pi;, &pi;].
460  **********************************************************************/
461  Math::real D(real sn, real cn, real dn) const throw();
462 
463  /**
464  * Legendre's geodesic longitude integral in terms of Jacobi elliptic
465  * functions.
466  *
467  * @param[in] sn = sin&phi;
468  * @param[in] cn = cos&phi;
469  * @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
470  * sin<sup>2</sup>&phi;)
471  * @return \e G(&phi;, &alpha;<sup>2</sup>, \e k) as though &phi; &isin;
472  * (&minus;&pi;, &pi;].
473  **********************************************************************/
474  Math::real G(real sn, real cn, real dn) const throw();
475 
476  /**
477  * Cayley's geodesic longitude difference integral in terms of Jacobi
478  * elliptic functions.
479  *
480  * @param[in] sn = sin&phi;
481  * @param[in] cn = cos&phi;
482  * @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
483  * sin<sup>2</sup>&phi;)
484  * @return \e H(&phi;, &alpha;<sup>2</sup>, \e k) as though &phi; &isin;
485  * (&minus;&pi;, &pi;].
486  **********************************************************************/
487  Math::real H(real sn, real cn, real dn) const throw();
488  ///@}
489 
490  /** \name Periodic versions of incomplete elliptic integrals.
491  **********************************************************************/
492  ///@{
493  /**
494  * The periodic incomplete integral of the first kind.
495  *
496  * @param[in] sn = sin&phi;
497  * @param[in] cn = cos&phi;
498  * @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
499  * sin<sup>2</sup>&phi;)
500  * @return the periodic function &pi; \e F(&phi;, \e k) / (2 \e K(\e k)) -
501  * &phi;
502  **********************************************************************/
503  Math::real deltaF(real sn, real cn, real dn) const throw();
504 
505  /**
506  * The periodic incomplete integral of the second kind.
507  *
508  * @param[in] sn = sin&phi;
509  * @param[in] cn = cos&phi;
510  * @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
511  * sin<sup>2</sup>&phi;)
512  * @return the periodic function &pi; \e E(&phi;, \e k) / (2 \e E(\e k)) -
513  * &phi;
514  **********************************************************************/
515  Math::real deltaE(real sn, real cn, real dn) const throw();
516 
517  /**
518  * The periodic inverse of the incomplete integral of the second kind.
519  *
520  * @param[in] stau = sin&tau;
521  * @param[in] ctau = sin&tau;
522  * @return the periodic function <i>E</i><sup>&minus;1</sup>(&tau; (2 \e
523  * E(\e k)/&pi;), \e k) - &tau;
524  **********************************************************************/
525  Math::real deltaEinv(real stau, real ctau) const throw();
526 
527  /**
528  * The periodic incomplete integral of the third kind.
529  *
530  * @param[in] sn = sin&phi;
531  * @param[in] cn = cos&phi;
532  * @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
533  * sin<sup>2</sup>&phi;)
534  * @return the periodic function &pi; &Pi;(&phi;, \e k) / (2 &Pi;(\e k)) -
535  * &phi;
536  **********************************************************************/
537  Math::real deltaPi(real sn, real cn, real dn) const throw();
538 
539  /**
540  * The periodic Jahnke's incomplete elliptic integral.
541  *
542  * @param[in] sn = sin&phi;
543  * @param[in] cn = cos&phi;
544  * @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
545  * sin<sup>2</sup>&phi;)
546  * @return the periodic function &pi; \e D(&phi;, \e k) / (2 \e D(\e k)) -
547  * &phi;
548  **********************************************************************/
549  Math::real deltaD(real sn, real cn, real dn) const throw();
550 
551  /**
552  * Legendre's periodic geodesic longitude integral.
553  *
554  * @param[in] sn = sin&phi;
555  * @param[in] cn = cos&phi;
556  * @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
557  * sin<sup>2</sup>&phi;)
558  * @return the periodic function &pi; \e G(&phi;, \e k) / (2 \e G(\e k)) -
559  * &phi;
560  **********************************************************************/
561  Math::real deltaG(real sn, real cn, real dn) const throw();
562 
563  /**
564  * Cayley's periodic geodesic longitude difference integral.
565  *
566  * @param[in] sn = sin&phi;
567  * @param[in] cn = cos&phi;
568  * @param[in] dn = sqrt(1 &minus; <i>k</i><sup>2</sup>
569  * sin<sup>2</sup>&phi;)
570  * @return the periodic function &pi; \e H(&phi;, \e k) / (2 \e H(\e k)) -
571  * &phi;
572  **********************************************************************/
573  Math::real deltaH(real sn, real cn, real dn) const throw();
574  ///@}
575 
576  /** \name Elliptic functions.
577  **********************************************************************/
578  ///@{
579  /**
580  * The Jacobi elliptic functions.
581  *
582  * @param[in] x the argument.
583  * @param[out] sn sn(\e x, \e k).
584  * @param[out] cn cn(\e x, \e k).
585  * @param[out] dn dn(\e x, \e k).
586  **********************************************************************/
587  void sncndn(real x, real& sn, real& cn, real& dn) const throw();
588 
589  /**
590  * The &Delta; amplitude function.
591  *
592  * @param[in] sn sin&phi;
593  * @param[in] cn cos&phi;
594  * @return &Delta; = sqrt(1 &minus; <i>k</i><sup>2</sup>
595  * sin<sup>2</sup>&phi;)
596  **********************************************************************/
597  Math::real Delta(real sn, real cn) const throw()
598  { return std::sqrt(_k2 < 0 ? 1 - _k2 * sn*sn : _kp2 + _k2 * cn*cn); }
599  ///@}
600 
601  /** \name Symmetric elliptic integrals.
602  **********************************************************************/
603  ///@{
604  /**
605  * Symmetric integral of the first kind <i>R<sub>F</sub></i>.
606  *
607  * @param[in] x
608  * @param[in] y
609  * @param[in] z
610  * @return <i>R<sub>F</sub></i>(\e x, \e y, \e z)
611  *
612  * <i>R<sub>F</sub></i> is defined in http://dlmf.nist.gov/19.16.E1
613  * \f[ R_F(x, y, z) = \frac12
614  * \int_0^\infty\frac1{\sqrt{(t + x) (t + y) (t + z)}}\, dt \f]
615  * If one of the arguments is zero, it is more efficient to call the
616  * two-argument version of this function with the non-zero arguments.
617  **********************************************************************/
618  static real RF(real x, real y, real z) throw();
619 
620  /**
621  * Complete symmetric integral of the first kind, <i>R<sub>F</sub></i> with
622  * one argument zero.
623  *
624  * @param[in] x
625  * @param[in] y
626  * @return <i>R<sub>F</sub></i>(\e x, \e y, 0)
627  **********************************************************************/
628  static real RF(real x, real y) throw();
629 
630  /**
631  * Degenerate symmetric integral of the first kind <i>R<sub>C</sub></i>.
632  *
633  * @param[in] x
634  * @param[in] y
635  * @return <i>R<sub>C</sub></i>(\e x, \e y) = <i>R<sub>F</sub></i>(\e x, \e
636  * y, \e y)
637  *
638  * <i>R<sub>C</sub></i> is defined in http://dlmf.nist.gov/19.2.E17
639  * \f[ R_C(x, y) = \frac12
640  * \int_0^\infty\frac1{\sqrt{t + x}(t + y)}\,dt \f]
641  **********************************************************************/
642  static real RC(real x, real y) throw();
643 
644  /**
645  * Symmetric integral of the second kind <i>R<sub>G</sub></i>.
646  *
647  * @param[in] x
648  * @param[in] y
649  * @param[in] z
650  * @return <i>R<sub>G</sub></i>(\e x, \e y, \e z)
651  *
652  * <i>R<sub>G</sub></i> is defined in Carlson, eq 1.5
653  * \f[ R_G(x, y, z) = \frac14
654  * \int_0^\infty[(t + x) (t + y) (t + z)]^{-1/2}
655  * \biggl(
656  * \frac x{t + x} + \frac y{t + y} + \frac z{t + z}
657  * \biggr)t\,dt \f]
658  * See also http://dlmf.nist.gov/19.16.E3.
659  * If one of the arguments is zero, it is more efficient to call the
660  * two-argument version of this function with the non-zero arguments.
661  **********************************************************************/
662  static real RG(real x, real y, real z) throw();
663 
664  /**
665  * Complete symmetric integral of the second kind, <i>R<sub>G</sub></i>
666  * with one argument zero.
667  *
668  * @param[in] x
669  * @param[in] y
670  * @return <i>R<sub>G</sub></i>(\e x, \e y, 0)
671  **********************************************************************/
672  static real RG(real x, real y) throw();
673 
674  /**
675  * Symmetric integral of the third kind <i>R<sub>J</sub></i>.
676  *
677  * @param[in] x
678  * @param[in] y
679  * @param[in] z
680  * @param[in] p
681  * @return <i>R<sub>J</sub></i>(\e x, \e y, \e z, \e p)
682  *
683  * <i>R<sub>J</sub></i> is defined in http://dlmf.nist.gov/19.16.E2
684  * \f[ R_J(x, y, z, p) = \frac32
685  * \int_0^\infty[(t + x) (t + y) (t + z)]^{-1/2} (t + p)^{-1}\, dt \f]
686  **********************************************************************/
687  static real RJ(real x, real y, real z, real p) throw();
688 
689  /**
690  * Degenerate symmetric integral of the third kind <i>R<sub>D</sub></i>.
691  *
692  * @param[in] x
693  * @param[in] y
694  * @param[in] z
695  * @return <i>R<sub>D</sub></i>(\e x, \e y, \e z) = <i>R<sub>J</sub></i>(\e
696  * x, \e y, \e z, \e z)
697  *
698  * <i>R<sub>D</sub></i> is defined in http://dlmf.nist.gov/19.16.E5
699  * \f[ R_D(x, y, z) = \frac32
700  * \int_0^\infty[(t + x) (t + y)]^{-1/2} (t + z)^{-3/2}\, dt \f]
701  **********************************************************************/
702  static real RD(real x, real y, real z) throw();
703  ///@}
704 
705  };
706 
707 } // namespace GeographicLib
708 
709 #endif // GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP
#define GEOGRAPHICLIB_EXPORT
Definition: Constants.hpp:52
void Reset(real k2=0, real alpha2=0)
GeographicLib::Math::real real
Definition: GeodSolve.cpp:40
Elliptic integrals and functions.
Math::real Delta(real sn, real cn) const
Header for GeographicLib::Constants class.