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