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