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