/** * \file NETGeographicLib/EllipticFunction.h * \brief Header for NETGeographicLib::EllipticFunction class * * NETGeographicLib is copyright (c) Scott Heiman (2013) * GeographicLib is Copyright (c) Charles Karney (2010-2012) * and licensed under the MIT/X11 License. * For more information, see * https://geographiclib.sourceforge.io/ **********************************************************************/ #pragma once namespace NETGeographicLib { /** * \brief .NET wrapper for GeographicLib::EllipticFunction. * * This class allows .NET applications to access GeographicLib::EllipticFunction. * * 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 * (https://dlmf.nist.gov/19.16.i) * - \e member functions to compute Legrendre's elliptic * integrals (https://dlmf.nist.gov/19.2.ii) and the * Jacobi elliptic functions (https://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 https://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 https://dlmf.nist.gov/19 and https://dlmf.nist.gov/22 * * C# Example: * \include example-EllipticFunction.cs * Managed C++ Example: * \include example-EllipticFunction.cpp * Visual Basic Example: * \include example-EllipticFunction.vb * * INTERFACE DIFFERENCES:
* The k2, kp2, alpha2, and alphap2 functions are implemented as properties. **********************************************************************/ public ref class EllipticFunction { private: // a pointer to the unmanaged GeographicLib::EllipticFunction. GeographicLib::EllipticFunction* m_pEllipticFunction; // The finalizer frees the unmanaged memory. !EllipticFunction(); 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(double k2, double 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(double k2, double alpha2, double kp2, double alphap2); /** * Destructor calls the finalizer. **********************************************************************/ ~EllipticFunction() { this->!EllipticFunction(); } /** * 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(double k2, double 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(double k2, double alpha2, double kp2, double alphap2); ///@} /** \name Inspector functions. **********************************************************************/ ///@{ /** * @return the square of the modulus k2. **********************************************************************/ property double k2 { double get(); } /** * @return the square of the complementary modulus k'2 = * 1 − k2. **********************************************************************/ property double kp2 { double get(); } /** * @return the parameter α2. **********************************************************************/ property double alpha2 { double get(); } /** * @return the complementary parameter α'2 = 1 − * α2. **********************************************************************/ property double alphap2 { double get(); } ///@} /** \name Complete elliptic integrals. **********************************************************************/ ///@{ /** * The complete integral of the first kind. * * @return \e K(\e k). * * \e K(\e k) is defined in https://dlmf.nist.gov/19.2.E4 * \f[ * K(k) = \int_0^{\pi/2} \frac1{\sqrt{1-k^2\sin^2\phi}}\,d\phi. * \f] **********************************************************************/ double K(); /** * The complete integral of the second kind. * * @return \e E(\e k) * * \e E(\e k) is defined in https://dlmf.nist.gov/19.2.E5 * \f[ * E(k) = \int_0^{\pi/2} \sqrt{1-k^2\sin^2\phi}\,d\phi. * \f] **********************************************************************/ double E(); /** * Jahnke's complete integral. * * @return \e D(\e k). * * \e D(\e k) is defined in https://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] **********************************************************************/ double D(); /** * The difference between the complete integrals of the first and second * kinds. * * @return \e K(\e k) − \e E(\e k). **********************************************************************/ double KE(); /** * The complete integral of the third kind. * * @return Π(α2, \e k) * * Π(α2, \e k) is defined in * https://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] **********************************************************************/ double Pi(); /** * 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] **********************************************************************/ double G(); /** * 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] **********************************************************************/ double H(); ///@} /** \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 https://dlmf.nist.gov/19.2.E4 * \f[ * F(\phi, k) = \int_0^\phi \frac1{\sqrt{1-k^2\sin^2\theta}}\,d\theta. * \f] **********************************************************************/ double F(double phi); /** * The incomplete integral of the second kind. * * @param[in] phi * @return \e E(φ, \e k). * * \e E(φ, \e k) is defined in https://dlmf.nist.gov/19.2.E5 * \f[ * E(\phi, k) = \int_0^\phi \sqrt{1-k^2\sin^2\theta}\,d\theta. * \f] **********************************************************************/ double E(double phi); /** * 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). **********************************************************************/ double Ed(double ang); /** * 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. **********************************************************************/ double Einv(double x); /** * The incomplete integral of the third kind. * * @param[in] phi * @return Π(φ, α2, \e k). * * Π(φ, α2, \e k) is defined in * https://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] **********************************************************************/ double Pi(double phi); /** * Jahnke's incomplete elliptic integral. * * @param[in] phi * @return \e D(φ, \e k). * * \e D(φ, \e k) is defined in https://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] **********************************************************************/ double D(double phi); /** * Legendre's geodesic longitude integral. * * @param[in] phi * @return \e G(φ, α2, \e k). * * \e G(φ, α2, \e k) is defined by * \f[ * \begin{aligned} * 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{aligned} * \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. **********************************************************************/ double G(double phi); /** * 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{aligned} * 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{aligned} * \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. **********************************************************************/ double H(double phi); ///@} /** \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 φ ∈ (−π, π]. **********************************************************************/ double F(double sn, double cn, double dn); /** * 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 φ ∈ (−π, π]. **********************************************************************/ double E(double sn, double cn, double dn); /** * 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 φ ∈ * (−π, π]. **********************************************************************/ double Pi(double sn, double cn, double dn); /** * 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 φ ∈ (−π, π]. **********************************************************************/ double D(double sn, double cn, double dn); /** * 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 φ ∈ * (−π, π]. **********************************************************************/ double G(double sn, double cn, double dn); /** * 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 φ ∈ * (−π, π]. **********************************************************************/ double H(double sn, double cn, double dn); ///@} /** \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)) - * φ **********************************************************************/ double deltaF(double sn, double cn, double dn); /** * 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)) - * φ **********************************************************************/ double deltaE(double sn, double cn, double dn); /** * 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) - τ **********************************************************************/ double deltaEinv(double stau, double ctau); /** * 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)) - * φ **********************************************************************/ double deltaPi(double sn, double cn, double dn); /** * 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)) - * φ **********************************************************************/ double deltaD(double sn, double cn, double dn); /** * 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)) - * φ **********************************************************************/ double deltaG(double sn, double cn, double dn); /** * 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)) - * φ **********************************************************************/ double deltaH(double sn, double cn, double dn); ///@} /** \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(double x, [System::Runtime::InteropServices::Out] double% sn, [System::Runtime::InteropServices::Out] double% cn, [System::Runtime::InteropServices::Out] double% dn); /** * The Δ amplitude function. * * @param[in] sn sinφ * @param[in] cn cosφ * @return Δ = sqrt(1 − k2 * sin2φ) **********************************************************************/ double Delta(double sn, double 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 https://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 double RF(double x, double y, double 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 double RF(double x, double 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 https://dlmf.nist.gov/19.2.E17 * \f[ R_C(x, y) = \frac12 * \int_0^\infty\frac1{\sqrt{t + x}(t + y)}\,dt \f] **********************************************************************/ static double RC(double x, double 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 https://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 double RG(double x, double y, double 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 double RG(double x, double 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 https://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 double RJ(double x, double y, double z, double 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 https://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 double RD(double x, double y, double z); ///@} }; } // namespace NETGeographicLib