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