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