/** * \file EllipticFunction.hpp * \brief Header for GeographicLib::EllipticFunction class * * Copyright (c) Charles Karney (2008-2012) and licensed * under the MIT/X11 License. For more information, see * http://geographiclib.sourceforge.net/ **********************************************************************/ #if !defined(GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP) #define GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP 1 #include namespace GeographicLib { /** * \brief Elliptic integrals and functions * * This provides the elliptic functions and integrals needed for Ellipsoid, * GeodesicExact, and TransverseMercatorExact. Two categories of function * are provided: * - \e static functions to compute symmetric elliptic integrals * (http://dlmf.nist.gov/19.16.i) * - \e member functions to compute Legrendre's elliptic * integrals (http://dlmf.nist.gov/19.2.ii) and the * Jacobi elliptic functions (http://dlmf.nist.gov/22.2). * . * In the latter case, an object is constructed giving the modulus \e k (and * optionally the parameter α2). The modulus is always * passed as its square k2 which allows \e k to be pure * imaginary (k2 < 0). (Confusingly, Abramowitz and * Stegun call \e m = k2 the "parameter" and \e n = * α2 the "characteristic".) * * In geodesic applications, it is convenient to separate the incomplete * integrals into secular and periodic components, e.g., * \f[ * E(\phi, k) = (2 E(\phi) / \pi) [ \phi + \delta E(\phi, k) ] * \f] * where δ\e E(φ, \e k) is an odd periodic function with period * π. * * The computation of the elliptic integrals uses the algorithms given in * - B. C. Carlson, * Computation of real or * complex elliptic integrals, Numerical Algorithms 10, 13--26 (1995) * . * with the additional optimizations given in http://dlmf.nist.gov/19.36.i. * The computation of the Jacobi elliptic functions uses the algorithm given * in * - R. Bulirsch, * Numerical Calculation of * Elliptic Integrals and Elliptic Functions, Numericshe Mathematik 7, * 78--90 (1965). * . * The notation follows http://dlmf.nist.gov/19 and http://dlmf.nist.gov/22 * * Example of use: * \include example-EllipticFunction.cpp **********************************************************************/ class GEOGRAPHICLIB_EXPORT EllipticFunction { private: typedef Math::real real; enum { num_ = 13 }; // Max depth required for sncndn. Probably 5 is enough. real _k2, _kp2, _alpha2, _alphap2, _eps; real _Kc, _Ec, _Dc, _Pic, _Gc, _Hc; public: /** \name Constructor **********************************************************************/ ///@{ /** * Constructor specifying the modulus and parameter. * * @param[in] k2 the square of the modulus k2. * k2 must lie in (-∞, 1). (No checking is * done.) * @param[in] alpha2 the parameter α2. * α2 must lie in (-∞, 1). (No checking is done.) * * If only elliptic integrals of the first and second kinds are needed, * then set α2 = 0 (the default value); in this case, we * have Π(φ, 0, \e k) = \e F(φ, \e k), \e G(φ, 0, \e k) = \e * E(φ, \e k), and \e H(φ, 0, \e k) = \e F(φ, \e k) - \e * D(φ, \e k). **********************************************************************/ EllipticFunction(real k2 = 0, real alpha2 = 0) { Reset(k2, alpha2); } /** * Constructor specifying the modulus and parameter and their complements. * * @param[in] k2 the square of the modulus k2. * k2 must lie in (-∞, 1). (No checking is * done.) * @param[in] alpha2 the parameter α2. * α2 must lie in (-∞, 1). (No checking is done.) * @param[in] kp2 the complementary modulus squared k'2 = * 1 − k2. * @param[in] alphap2 the complementary parameter α'2 = 1 * − α2. * * The arguments must satisfy \e k2 + \e kp2 = 1 and \e alpha2 + \e alphap2 * = 1. (No checking is done that these conditions are met.) This * constructor is provided to enable accuracy to be maintained, e.g., when * \e k is very close to unity. **********************************************************************/ EllipticFunction(real k2, real alpha2, real kp2, real alphap2) { Reset(k2, alpha2, kp2, alphap2); } /** * Reset the modulus and parameter. * * @param[in] k2 the new value of square of the modulus * k2 which must lie in (-∞, 1). (No checking is * done.) * @param[in] alpha2 the new value of parameter α2. * α2 must lie in (-∞, 1). (No checking is done.) **********************************************************************/ void Reset(real k2 = 0, real alpha2 = 0) { Reset(k2, alpha2, 1 - k2, 1 - alpha2); } /** * Reset the modulus and parameter supplying also their complements. * * @param[in] k2 the square of the modulus k2. * k2 must lie in (-∞, 1). (No checking is * done.) * @param[in] alpha2 the parameter α2. * α2 must lie in (-∞, 1). (No checking is done.) * @param[in] kp2 the complementary modulus squared k'2 = * 1 − k2. * @param[in] alphap2 the complementary parameter α'2 = 1 * − α2. * * The arguments must satisfy \e k2 + \e kp2 = 1 and \e alpha2 + \e alphap2 * = 1. (No checking is done that these conditions are met.) This * constructor is provided to enable accuracy to be maintained, e.g., when * is very small. **********************************************************************/ void Reset(real k2, real alpha2, real kp2, real alphap2); ///@} /** \name Inspector functions. **********************************************************************/ ///@{ /** * @return the square of the modulus k2. **********************************************************************/ Math::real k2() const { return _k2; } /** * @return the square of the complementary modulus k'2 = * 1 − k2. **********************************************************************/ Math::real kp2() const { return _kp2; } /** * @return the parameter α2. **********************************************************************/ Math::real alpha2() const { return _alpha2; } /** * @return the complementary parameter α'2 = 1 − * α2. **********************************************************************/ Math::real alphap2() const { return _alphap2; } ///@} /** \name Complete elliptic integrals. **********************************************************************/ ///@{ /** * The complete integral of the first kind. * * @return \e K(\e k). * * \e K(\e k) is defined in http://dlmf.nist.gov/19.2.E4 * \f[ * K(k) = \int_0^{\pi/2} \frac1{\sqrt{1-k^2\sin^2\phi}}\,d\phi. * \f] **********************************************************************/ Math::real K() const { return _Kc; } /** * The complete integral of the second kind. * * @return \e E(\e k) * * \e E(\e k) is defined in http://dlmf.nist.gov/19.2.E5 * \f[ * E(k) = \int_0^{\pi/2} \sqrt{1-k^2\sin^2\phi}\,d\phi. * \f] **********************************************************************/ Math::real E() const { return _Ec; } /** * Jahnke's complete integral. * * @return \e D(\e k). * * \e D(\e k) is defined in http://dlmf.nist.gov/19.2.E6 * \f[ * D(k) = \int_0^{\pi/2} \frac{\sin^2\phi}{\sqrt{1-k^2\sin^2\phi}}\,d\phi. * \f] **********************************************************************/ Math::real D() const { return _Dc; } /** * The difference between the complete integrals of the first and second * kinds. * * @return \e K(\e k) − \e E(\e k). **********************************************************************/ Math::real KE() const { return _k2 * _Dc; } /** * The complete integral of the third kind. * * @return Π(α2, \e k) * * Π(α2, \e k) is defined in * http://dlmf.nist.gov/19.2.E7 * \f[ * \Pi(\alpha^2, k) = \int_0^{\pi/2} * \frac1{\sqrt{1-k^2\sin^2\phi}(1 - \alpha^2\sin^2\phi)}\,d\phi. * \f] **********************************************************************/ Math::real Pi() const { return _Pic; } /** * Legendre's complete geodesic longitude integral. * * @return \e G(α2, \e k) * * \e G(α2, \e k) is given by * \f[ * G(\alpha^2, k) = \int_0^{\pi/2} * \frac{\sqrt{1-k^2\sin^2\phi}}{1 - \alpha^2\sin^2\phi}\,d\phi. * \f] **********************************************************************/ Math::real G() const { return _Gc; } /** * Cayley's complete geodesic longitude difference integral. * * @return \e H(α2, \e k) * * \e H(α2, \e k) is given by * \f[ * H(\alpha^2, k) = \int_0^{\pi/2} * \frac{\cos^2\phi}{(1-\alpha^2\sin^2\phi)\sqrt{1-k^2\sin^2\phi}} * \,d\phi. * \f] **********************************************************************/ Math::real H() const { return _Hc; } ///@} /** \name Incomplete elliptic integrals. **********************************************************************/ ///@{ /** * The incomplete integral of the first kind. * * @param[in] phi * @return \e F(φ, \e k). * * \e F(φ, \e k) is defined in http://dlmf.nist.gov/19.2.E4 * \f[ * F(\phi, k) = \int_0^\phi \frac1{\sqrt{1-k^2\sin^2\theta}}\,d\theta. * \f] **********************************************************************/ Math::real F(real phi) const; /** * The incomplete integral of the second kind. * * @param[in] phi * @return \e E(φ, \e k). * * \e E(φ, \e k) is defined in http://dlmf.nist.gov/19.2.E5 * \f[ * E(\phi, k) = \int_0^\phi \sqrt{1-k^2\sin^2\theta}\,d\theta. * \f] **********************************************************************/ Math::real E(real phi) const; /** * The incomplete integral of the second kind with the argument given in * degrees. * * @param[in] ang in degrees. * @return \e E(π ang/180, \e k). **********************************************************************/ Math::real Ed(real ang) const; /** * The inverse of the incomplete integral of the second kind. * * @param[in] x * @return φ = E−1(\e x, \e k); i.e., the * solution of such that \e E(φ, \e k) = \e x. **********************************************************************/ Math::real Einv(real x) const; /** * The incomplete integral of the third kind. * * @param[in] phi * @return Π(φ, α2, \e k). * * Π(φ, α2, \e k) is defined in * http://dlmf.nist.gov/19.2.E7 * \f[ * \Pi(\phi, \alpha^2, k) = \int_0^\phi * \frac1{\sqrt{1-k^2\sin^2\theta}(1 - \alpha^2\sin^2\theta)}\,d\theta. * \f] **********************************************************************/ Math::real Pi(real phi) const; /** * Jahnke's incomplete elliptic integral. * * @param[in] phi * @return \e D(φ, \e k). * * \e D(φ, \e k) is defined in http://dlmf.nist.gov/19.2.E4 * \f[ * D(\phi, k) = \int_0^\phi * \frac{\sin^2\theta}{\sqrt{1-k^2\sin^2\theta}}\,d\theta. * \f] **********************************************************************/ Math::real D(real phi) const; /** * Legendre's geodesic longitude integral. * * @param[in] phi * @return \e G(φ, α2, \e k). * * \e G(φ, α2, \e k) is defined by * \f[ * \begin{align} * G(\phi, \alpha^2, k) &= * \frac{k^2}{\alpha^2} F(\phi, k) + * \biggl(1 - \frac{k^2}{\alpha^2}\biggr) \Pi(\phi, \alpha^2, k) \\ * &= \int_0^\phi * \frac{\sqrt{1-k^2\sin^2\theta}}{1 - \alpha^2\sin^2\theta}\,d\theta. * \end{align} * \f] * * Legendre expresses the longitude of a point on the geodesic in terms of * this combination of elliptic integrals in Exercices de Calcul * Intégral, Vol. 1 (1811), p. 181, * https://books.google.com/books?id=riIOAAAAQAAJ&pg=PA181. * * See \ref geodellip for the expression for the longitude in terms of this * function. **********************************************************************/ Math::real G(real phi) const; /** * Cayley's geodesic longitude difference integral. * * @param[in] phi * @return \e H(φ, α2, \e k). * * \e H(φ, α2, \e k) is defined by * \f[ * \begin{align} * H(\phi, \alpha^2, k) &= * \frac1{\alpha^2} F(\phi, k) + * \biggl(1 - \frac1{\alpha^2}\biggr) \Pi(\phi, \alpha^2, k) \\ * &= \int_0^\phi * \frac{\cos^2\theta}{(1-\alpha^2\sin^2\theta)\sqrt{1-k^2\sin^2\theta}} * \,d\theta. * \end{align} * \f] * * Cayley expresses the longitude difference of a point on the geodesic in * terms of this combination of elliptic integrals in Phil. Mag. 40 * (1870), p. 333, https://books.google.com/books?id=Zk0wAAAAIAAJ&pg=PA333. * * See \ref geodellip for the expression for the longitude in terms of this * function. **********************************************************************/ Math::real H(real phi) const; ///@} /** \name Incomplete integrals in terms of Jacobi elliptic functions. **********************************************************************/ /** * The incomplete integral of the first kind in terms of Jacobi elliptic * functions. * * @param[in] sn = sinφ * @param[in] cn = cosφ * @param[in] dn = sqrt(1 − k2 * sin2φ) * @return \e F(φ, \e k) as though φ ∈ (−π, π]. **********************************************************************/ Math::real F(real sn, real cn, real dn) const; /** * The incomplete integral of the second kind in terms of Jacobi elliptic * functions. * * @param[in] sn = sinφ * @param[in] cn = cosφ * @param[in] dn = sqrt(1 − k2 * sin2φ) * @return \e E(φ, \e k) as though φ ∈ (−π, π]. **********************************************************************/ Math::real E(real sn, real cn, real dn) const; /** * The incomplete integral of the third kind in terms of Jacobi elliptic * functions. * * @param[in] sn = sinφ * @param[in] cn = cosφ * @param[in] dn = sqrt(1 − k2 * sin2φ) * @return Π(φ, α2, \e k) as though φ ∈ * (−π, π]. **********************************************************************/ Math::real Pi(real sn, real cn, real dn) const; /** * Jahnke's incomplete elliptic integral in terms of Jacobi elliptic * functions. * * @param[in] sn = sinφ * @param[in] cn = cosφ * @param[in] dn = sqrt(1 − k2 * sin2φ) * @return \e D(φ, \e k) as though φ ∈ (−π, π]. **********************************************************************/ Math::real D(real sn, real cn, real dn) const; /** * Legendre's geodesic longitude integral in terms of Jacobi elliptic * functions. * * @param[in] sn = sinφ * @param[in] cn = cosφ * @param[in] dn = sqrt(1 − k2 * sin2φ) * @return \e G(φ, α2, \e k) as though φ ∈ * (−π, π]. **********************************************************************/ Math::real G(real sn, real cn, real dn) const; /** * Cayley's geodesic longitude difference integral in terms of Jacobi * elliptic functions. * * @param[in] sn = sinφ * @param[in] cn = cosφ * @param[in] dn = sqrt(1 − k2 * sin2φ) * @return \e H(φ, α2, \e k) as though φ ∈ * (−π, π]. **********************************************************************/ Math::real H(real sn, real cn, real dn) const; ///@} /** \name Periodic versions of incomplete elliptic integrals. **********************************************************************/ ///@{ /** * The periodic incomplete integral of the first kind. * * @param[in] sn = sinφ * @param[in] cn = cosφ * @param[in] dn = sqrt(1 − k2 * sin2φ) * @return the periodic function π \e F(φ, \e k) / (2 \e K(\e k)) - * φ **********************************************************************/ Math::real deltaF(real sn, real cn, real dn) const; /** * The periodic incomplete integral of the second kind. * * @param[in] sn = sinφ * @param[in] cn = cosφ * @param[in] dn = sqrt(1 − k2 * sin2φ) * @return the periodic function π \e E(φ, \e k) / (2 \e E(\e k)) - * φ **********************************************************************/ Math::real deltaE(real sn, real cn, real dn) const; /** * The periodic inverse of the incomplete integral of the second kind. * * @param[in] stau = sinτ * @param[in] ctau = sinτ * @return the periodic function E−1(τ (2 \e * E(\e k)/π), \e k) - τ **********************************************************************/ Math::real deltaEinv(real stau, real ctau) const; /** * The periodic incomplete integral of the third kind. * * @param[in] sn = sinφ * @param[in] cn = cosφ * @param[in] dn = sqrt(1 − k2 * sin2φ) * @return the periodic function π Π(φ, \e k) / (2 Π(\e k)) - * φ **********************************************************************/ Math::real deltaPi(real sn, real cn, real dn) const; /** * The periodic Jahnke's incomplete elliptic integral. * * @param[in] sn = sinφ * @param[in] cn = cosφ * @param[in] dn = sqrt(1 − k2 * sin2φ) * @return the periodic function π \e D(φ, \e k) / (2 \e D(\e k)) - * φ **********************************************************************/ Math::real deltaD(real sn, real cn, real dn) const; /** * Legendre's periodic geodesic longitude integral. * * @param[in] sn = sinφ * @param[in] cn = cosφ * @param[in] dn = sqrt(1 − k2 * sin2φ) * @return the periodic function π \e G(φ, \e k) / (2 \e G(\e k)) - * φ **********************************************************************/ Math::real deltaG(real sn, real cn, real dn) const; /** * Cayley's periodic geodesic longitude difference integral. * * @param[in] sn = sinφ * @param[in] cn = cosφ * @param[in] dn = sqrt(1 − k2 * sin2φ) * @return the periodic function π \e H(φ, \e k) / (2 \e H(\e k)) - * φ **********************************************************************/ Math::real deltaH(real sn, real cn, real dn) const; ///@} /** \name Elliptic functions. **********************************************************************/ ///@{ /** * The Jacobi elliptic functions. * * @param[in] x the argument. * @param[out] sn sn(\e x, \e k). * @param[out] cn cn(\e x, \e k). * @param[out] dn dn(\e x, \e k). **********************************************************************/ void sncndn(real x, real& sn, real& cn, real& dn) const; /** * The Δ amplitude function. * * @param[in] sn sinφ * @param[in] cn cosφ * @return Δ = sqrt(1 − k2 * sin2φ) **********************************************************************/ Math::real Delta(real sn, real cn) const { using std::sqrt; return sqrt(_k2 < 0 ? 1 - _k2 * sn*sn : _kp2 + _k2 * cn*cn); } ///@} /** \name Symmetric elliptic integrals. **********************************************************************/ ///@{ /** * Symmetric integral of the first kind RF. * * @param[in] x * @param[in] y * @param[in] z * @return RF(\e x, \e y, \e z) * * RF is defined in http://dlmf.nist.gov/19.16.E1 * \f[ R_F(x, y, z) = \frac12 * \int_0^\infty\frac1{\sqrt{(t + x) (t + y) (t + z)}}\, dt \f] * If one of the arguments is zero, it is more efficient to call the * two-argument version of this function with the non-zero arguments. **********************************************************************/ static real RF(real x, real y, real z); /** * Complete symmetric integral of the first kind, * RF with one argument zero. * * @param[in] x * @param[in] y * @return RF(\e x, \e y, 0) **********************************************************************/ static real RF(real x, real y); /** * Degenerate symmetric integral of the first kind * RC. * * @param[in] x * @param[in] y * @return RC(\e x, \e y) = * RF(\e x, \e y, \e y) * * RC is defined in http://dlmf.nist.gov/19.2.E17 * \f[ R_C(x, y) = \frac12 * \int_0^\infty\frac1{\sqrt{t + x}(t + y)}\,dt \f] **********************************************************************/ static real RC(real x, real y); /** * Symmetric integral of the second kind RG. * * @param[in] x * @param[in] y * @param[in] z * @return RG(\e x, \e y, \e z) * * RG is defined in Carlson, eq 1.5 * \f[ R_G(x, y, z) = \frac14 * \int_0^\infty[(t + x) (t + y) (t + z)]^{-1/2} * \biggl( * \frac x{t + x} + \frac y{t + y} + \frac z{t + z} * \biggr)t\,dt \f] * See also http://dlmf.nist.gov/19.16.E3. * If one of the arguments is zero, it is more efficient to call the * two-argument version of this function with the non-zero arguments. **********************************************************************/ static real RG(real x, real y, real z); /** * Complete symmetric integral of the second kind, * RG with one argument zero. * * @param[in] x * @param[in] y * @return RG(\e x, \e y, 0) **********************************************************************/ static real RG(real x, real y); /** * Symmetric integral of the third kind RJ. * * @param[in] x * @param[in] y * @param[in] z * @param[in] p * @return RJ(\e x, \e y, \e z, \e p) * * RJ is defined in http://dlmf.nist.gov/19.16.E2 * \f[ R_J(x, y, z, p) = \frac32 * \int_0^\infty[(t + x) (t + y) (t + z)]^{-1/2} (t + p)^{-1}\, dt \f] **********************************************************************/ static real RJ(real x, real y, real z, real p); /** * Degenerate symmetric integral of the third kind * RD. * * @param[in] x * @param[in] y * @param[in] z * @return RD(\e x, \e y, \e z) = * RJ(\e x, \e y, \e z, \e z) * * RD is defined in http://dlmf.nist.gov/19.16.E5 * \f[ R_D(x, y, z) = \frac32 * \int_0^\infty[(t + x) (t + y)]^{-1/2} (t + z)^{-3/2}\, dt \f] **********************************************************************/ static real RD(real x, real y, real z); ///@} }; } // namespace GeographicLib #endif // GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP