GeographicLib
1.21
|
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