GeographicLib  1.21
EllipticFunction.cpp
Go to the documentation of this file.
00001 /**
00002  * \file EllipticFunction.cpp
00003  * \brief Implementation for GeographicLib::EllipticFunction class
00004  *
00005  * Copyright (c) Charles Karney (2008-2011) <charles@karney.com> and licensed
00006  * under the MIT/X11 License.  For more information, see
00007  * http://geographiclib.sourceforge.net/
00008  **********************************************************************/
00009 
00010 #include <GeographicLib/EllipticFunction.hpp>
00011 
00012 #define GEOGRAPHICLIB_ELLIPTICFUNCTION_CPP \
00013   "$Id: 00b30b3d051fce1da7eb0c7e74c1c03854de6ea3 $"
00014 
00015 RCSID_DECL(GEOGRAPHICLIB_ELLIPTICFUNCTION_CPP)
00016 RCSID_DECL(GEOGRAPHICLIB_ELLIPTICFUNCTION_HPP)
00017 
00018 namespace GeographicLib {
00019 
00020   using namespace std;
00021 
00022   const Math::real EllipticFunction::tol_ =
00023     numeric_limits<real>::epsilon() * real(0.01);
00024   const Math::real EllipticFunction::tolRF_ = pow(3 * tol_, 1/real(6));
00025   const Math::real EllipticFunction::tolRD_ =
00026     pow(real(0.25) * tol_, 1/real(6));
00027   const Math::real EllipticFunction::tolRG0_ = real(2.7) * sqrt(tol_);
00028   const Math::real EllipticFunction::tolJAC_ = sqrt(tol_);
00029   const Math::real EllipticFunction::tolJAC1_ = sqrt(6 * tol_);
00030 
00031   /*
00032    * Implementation of methods given in
00033    *
00034    *   B. C. Carlson
00035    *   Computation of elliptic integrals
00036    *   Numerical Algorithms 10, 13-26 (1995)
00037    */
00038 
00039   Math::real EllipticFunction::RF(real x, real y, real z) throw() {
00040     // Carlson, eqs 2.2 - 2.7
00041     real
00042       a0 = (x + y + z)/3,
00043       an = a0,
00044       q = max(max(abs(a0-x), abs(a0-y)), abs(a0-z)) / tolRF_,
00045       x0 = x,
00046       y0 = y,
00047       z0 = z,
00048       mul = 1;
00049     while (q >= mul * abs(an)) {
00050       // Max 6 trips
00051       real ln = sqrt(x0)*sqrt(y0) + sqrt(y0)*sqrt(z0) + sqrt(z0)*sqrt(x0);
00052       an = (an + ln)/4;
00053       x0 = (x0 + ln)/4;
00054       y0 = (y0 + ln)/4;
00055       z0 = (z0 + ln)/4;
00056       mul *= 4;
00057     }
00058     real
00059       xx = (a0 - x) / (mul * an),
00060       yy = (a0 - y) / (mul * an),
00061       zz = - xx - yy,
00062       e2 = xx * yy - zz * zz,
00063       e3 = xx * yy * zz;
00064     return (1 - e2 / 10 + e3 / 14 + e2 * e2 / 24 - 3 * e2 * e3 / 44) / sqrt(an);
00065   }
00066 
00067   Math::real EllipticFunction::RD(real x, real y, real z) throw() {
00068     // Carlson, eqs 2.28 - 2.34
00069     real
00070       a0 = (x + y + 3 * z)/5,
00071       an = a0,
00072       q = max(max(abs(a0-x), abs(a0-y)), abs(a0-z)) / tolRD_,
00073       x0 = x,
00074       y0 = y,
00075       z0 = z,
00076       mul = 1,
00077       s = 0;
00078     while (q >= mul * abs(an)) {
00079       // Max 7 trips
00080       real ln = sqrt(x0)*sqrt(y0) +
00081         sqrt(y0)*sqrt(z0) +
00082         sqrt(z0)*sqrt(x0);
00083       s += 1/(mul * sqrt(z0) * (z0 + ln ));
00084       an = (an + ln)/4;
00085       x0 = (x0 + ln)/4;
00086       y0 = (y0 + ln)/4;
00087       z0 = (z0 + ln)/4;
00088       mul *= 4;
00089     }
00090     real
00091       xx = (a0 - x) / (mul * an),
00092       yy = (a0 - y) / (mul * an),
00093       zz = -(xx + yy) / 3,
00094       e2 = xx * yy - 6 * zz * zz,
00095       e3 = (3 * xx * yy - 8 * zz * zz)*zz,
00096       e4 = 3 * (xx * yy - zz * zz) * zz * zz,
00097       e5 = xx * yy * zz * zz * zz;
00098     return (1 - 3 * e2 / 14 + e3 / 6 + 9 * e2 * e2 / 88 - 3 * e4 / 22
00099             - 9 * e2 * e3 / 52 + 3 * e5 / 26) / (mul * an * sqrt(an))
00100       + 3 * s;
00101   }
00102 
00103   Math::real EllipticFunction::RG0(real x, real y) throw() {
00104     // Carlson, eqs 2.36 - 2.39
00105     real
00106       x0 = sqrt(x),
00107       y0 = sqrt(y),
00108       xn = x0,
00109       yn = y0,
00110       s = 0,
00111       mul = real(0.25);
00112     while (abs(xn-yn) >= tolRG0_ * abs(xn)) {
00113       // Max 4 trips
00114       real t = (xn + yn) /2;
00115       yn = sqrt(xn * yn);
00116       xn = t;
00117       mul *= 2;
00118       t = xn - yn;
00119       s += mul * t * t;
00120     }
00121     x0 = (x0 + y0)/2;
00122     return (x0 * x0 - s) * Math::pi<real>() / (2 * (xn + yn));
00123   }
00124 
00125   EllipticFunction::EllipticFunction(real m) throw()
00126     : _m(m)
00127     , _m1(1 - m)
00128       // Don't initialize _kc, _ec, _kec since this constructor might be called
00129       // before the static real constants tolRF_, etc., are initialized.
00130     , _init(false)
00131   {}
00132 
00133   bool EllipticFunction::Init() const throw() {
00134     // Complete elliptic integral K(m), Carlson eq. 4.1
00135     _kc = RF(real(0), _m1, real(1));
00136     // Complete elliptic integral E(m), Carlson eq. 4.2
00137     _ec = 2 * RG0(_m1, real(1));
00138     // K - E, Carlson eq.4.3
00139     _kec = _m / 3 * RD(real(0), _m1, real(1));
00140     return _init = true;
00141   }
00142 
00143   /*
00144    * Implementation of methods given in
00145    *
00146    *   R. Bulirsch
00147    *   Numerical Calculation of Elliptic Integrals and Elliptic Functions
00148    *   Numericshe Mathematik 7, 78-90 (1965)
00149    */
00150 
00151   void EllipticFunction::sncndn(real x, real& sn, real& cn, real& dn)
00152     const throw() {
00153     // Bulirsch's sncndn routine, p 89.
00154     //
00155     // Assume _m1 is in [0, 1].  See Bulirsch article for code to treat
00156     // negative _m1.
00157     if (_m1 != 0) {
00158       real mc = _m1;
00159       real c;
00160       real m[num_], n[num_];
00161       unsigned l = 0;
00162       for (real a = 1; l < num_; ++l) {
00163         // Max 5 trips
00164         m[l] = a;
00165         n[l] = mc = sqrt(mc);
00166         c = (a + mc) / 2;
00167         if (!(abs(a - mc) > tolJAC_ * a)) {
00168           ++l;
00169           break;
00170         }
00171         mc = a * mc;
00172         a = c;
00173       }
00174       x = c * x;
00175       sn = sin(x);
00176       cn = cos(x);
00177       dn = 1;
00178       if (sn != 0) {
00179         real a = cn / sn;
00180         c = a * c;
00181         while (l--) {
00182           real b = m[l];
00183           a = c * a;
00184           c = dn * c;
00185           dn = (n[l] + a) / (b + a);
00186           a = c / b;
00187         }
00188         a = 1 / sqrt(c * c + 1);
00189         sn = sn < 0 ? -a : a;
00190         cn = c * sn;
00191       }
00192     } else {
00193       sn = tanh(x);
00194       dn = cn = 1 / cosh(x);
00195     }
00196   }
00197 
00198   Math::real EllipticFunction::E(real sn, real cn, real dn) const throw() {
00199     real
00200       cn2 = cn * cn, dn2 = dn * dn, sn2 = sn * sn,
00201       // Carlson, eq. 4.6
00202       ei = abs(sn) * (RF(cn2, dn2, real(1)) -
00203                       (_m / 3) * sn2 * RD(cn2, dn2, real(1)));
00204     // Enforce usual trig-like symmetries
00205     if (cn < 0) {
00206       ei = 2 * E() - ei;
00207     }
00208     if (sn < 0)
00209       ei = -ei;
00210     return ei;
00211   }
00212 
00213   Math::real EllipticFunction::E(real phi) const throw() {
00214     real sn = sin(phi);
00215     return E(sn, cos(phi), sqrt(1 - _m * sn * sn));
00216   }
00217 
00218 } // namespace GeographicLib