00001
00002
00003
00004
00005
00006
00007
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
00032
00033
00034
00035
00036
00037
00038 Math::real EllipticFunction::RF(real x, real y, real z) throw() {
00039
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
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
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
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
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
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
00128
00129 , _init(false)
00130 {}
00131
00132 bool EllipticFunction::Init() const throw() {
00133
00134 _kc = RF(real(0), _m1, real(1));
00135
00136 _ec = 2 * RG0(_m1, real(1));
00137
00138 _kec = _m / 3 * RD(real(0), _m1, real(1));
00139 return _init = true;
00140 }
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150 void EllipticFunction::sncndn(real x, real& sn, real& cn, real& dn)
00151 const throw() {
00152
00153
00154
00155
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
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
00201 ei = abs(sn) * (RF(cn2, dn2, real(1)) -
00202 (_m / 3) * sn2 * RD(cn2, dn2, real(1)));
00203
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 }