1 /** 2 * Elliptic integrals. 3 * The functions are named similarly to the names used in Mathematica. 4 * 5 * Eric W. Weisstein. "Elliptic Integral of the First Kind." From MathWorld--A Wolfram Web Resource. $(LINK http://mathworld.wolfram.com/EllipticIntegraloftheFirstKind.html) 6 * 7 * $(LINK http://www.netlib.org/cephes/ldoubdoc.html) 8 * 9 * Copyright: 10 * Based on the CEPHES math library, which is 11 * Copyright (C) 1994 Stephen L. Moshier (moshier@world.std.com). 12 * Some parts copyright (c) 2009-2016 dunnhumby Germany GmbH. 13 * All rights reserved. 14 * 15 * License: 16 * Tango Dual License: 3-Clause BSD License / Academic Free License v3.0. 17 * See LICENSE_TANGO.txt for details. 18 * 19 * Authors: Stephen L. Moshier (original C code). Conversion to D by Don Clugston 20 * 21 * References: $(LINK http://en.wikipedia.org/wiki/Elliptic_integral) 22 * 23 * Macros: 24 * TABLE_SV = <table border=1 cellpadding=4 cellspacing=0> 25 * <caption>Special Values</caption> 26 * $0</table> 27 * SVH = $(TR $(TH $1) $(TH $2)) 28 * SV = $(TR $(TD $1) $(TD $2)) 29 * GAMMA = Γ 30 * INTEGRATE = $(BIG ∫<sub>$(SMALL $1)</sub><sup>$2</sup>) 31 * POWER = $1<sup>$2</sup> 32 * NAN = $(RED NAN) 33 */ 34 /** 35 * Macros: 36 * TABLE_SV = <table border=1 cellpadding=4 cellspacing=0> 37 * <caption>Special Values</caption> 38 * $0</table> 39 * SVH = $(TR $(TH $1) $(TH $2)) 40 * SV = $(TR $(TD $1) $(TD $2)) 41 * 42 * NAN = $(RED NAN) 43 * SUP = <span style="vertical-align:super;font-size:smaller">$0</span> 44 * GAMMA = Γ 45 * INTEGRAL = ∫ 46 * INTEGRATE = $(BIG ∫<sub>$(SMALL $1)</sub><sup>$2</sup>) 47 * POWER = $1<sup>$2</sup> 48 * BIGSUM = $(BIG Σ <sup>$2</sup><sub>$(SMALL $1)</sub>) 49 * CHOOSE = $(BIG () <sup>$(SMALL $1)</sup><sub>$(SMALL $2)</sub> $(BIG )) 50 */ 51 52 module ocean.math.Elliptic; 53 54 import ocean.math.Math; 55 import ocean.math.IEEE; 56 import ocean.core.Verify; 57 58 version (unittest) import ocean.core.Test; 59 60 /* These functions are based on code from: 61 Cephes Math Library, Release 2.3: October, 1995 62 Copyright 1984, 1987, 1995 by Stephen L. Moshier 63 */ 64 65 66 /** 67 * Incomplete elliptic integral of the first kind 68 * 69 * Approximates the integral 70 * F(phi | m) = $(INTEGRATE 0, phi) dt/ (sqrt( 1- m $(POWER sin, 2) t)) 71 * 72 * of amplitude phi and modulus m, using the arithmetic - 73 * geometric mean algorithm. 74 */ 75 76 real ellipticF(real phi, real m ) 77 { 78 real a, b, c, e, temp, t, K; 79 int d, mod, sign, npio2; 80 81 if( m == 0.0L ) 82 return phi; 83 a = 1.0L - m; 84 if( a == 0.0L ) { 85 if ( fabs(phi) >= PI_2 ) return real.infinity; 86 return log( tan( 0.5L*(PI_2 + phi) ) ); 87 } 88 npio2 = cast(int)floor( phi/PI_2 ); 89 if ( npio2 & 1 ) 90 npio2 += 1; 91 if ( npio2 ) { 92 K = ellipticKComplete( a ); 93 phi = phi - npio2 * PI_2; 94 } else 95 K = 0.0L; 96 if( phi < 0.0L ){ 97 phi = -phi; 98 sign = -1; 99 } else sign = 0; 100 b = sqrt(a); 101 t = tan( phi ); 102 if( fabs(t) > 10.0L ) { 103 /* Transform the amplitude */ 104 e = 1.0L/(b*t); 105 /* ... but avoid multiple recursions. */ 106 if( fabs(e) < 10.0L ){ 107 e = atan(e); 108 if( npio2 == 0 ) 109 K = ellipticKComplete( a ); 110 temp = K - ellipticF( e, m ); 111 goto done; 112 } 113 } 114 a = 1.0L; 115 c = sqrt(m); 116 d = 1; 117 mod = 0; 118 119 while( fabs(c/a) > real.epsilon ) { 120 temp = b/a; 121 phi = phi + atan(t*temp) + mod * PI; 122 mod = cast(int)((phi + PI_2)/PI); 123 t = t * ( 1.0L + temp )/( 1.0L - temp * t * t ); 124 c = 0.5L * ( a - b ); 125 temp = sqrt( a * b ); 126 a = 0.5L * ( a + b ); 127 b = temp; 128 d += d; 129 } 130 temp = (atan(t) + mod * PI)/(d * a); 131 132 done: 133 if ( sign < 0 ) 134 temp = -temp; 135 temp += npio2 * K; 136 return temp; 137 } 138 139 140 /** 141 * Incomplete elliptic integral of the second kind 142 * 143 * Approximates the integral 144 * 145 * E(phi | m) = $(INTEGRATE 0, phi) sqrt( 1- m $(POWER sin, 2) t) dt 146 * 147 * of amplitude phi and modulus m, using the arithmetic - 148 * geometric mean algorithm. 149 */ 150 151 real ellipticE(real phi, real m) 152 { 153 real a, b, c, e, temp, t, E; 154 int d, mod, npio2, sign; 155 156 if ( m == 0.0L ) return phi; 157 real lphi = phi; 158 npio2 = cast(int)floor( lphi/PI_2 ); 159 if( npio2 & 1 ) 160 npio2 += 1; 161 lphi = lphi - npio2 * PI_2; 162 if( lphi < 0.0L ){ 163 lphi = -lphi; 164 sign = -1; 165 } else { 166 sign = 1; 167 } 168 a = 1.0L - m; 169 E = ellipticEComplete( a ); 170 if( a == 0.0L ) { 171 temp = sin( lphi ); 172 goto done; 173 } 174 t = tan( lphi ); 175 b = sqrt(a); 176 if ( fabs(t) > 10.0L ) { 177 /* Transform the amplitude */ 178 e = 1.0L/(b*t); 179 /* ... but avoid multiple recursions. */ 180 if( fabs(e) < 10.0L ){ 181 e = atan(e); 182 temp = E + m * sin( lphi ) * sin( e ) - ellipticE( e, m ); 183 goto done; 184 } 185 } 186 c = sqrt(m); 187 a = 1.0L; 188 d = 1; 189 e = 0.0L; 190 mod = 0; 191 192 while( fabs(c/a) > real.epsilon ) { 193 temp = b/a; 194 lphi = lphi + atan(t*temp) + mod * PI; 195 mod = cast(int)((lphi + PI_2)/PI); 196 t = t * ( 1.0L + temp )/( 1.0L - temp * t * t ); 197 c = 0.5L*( a - b ); 198 temp = sqrt( a * b ); 199 a = 0.5L*( a + b ); 200 b = temp; 201 d += d; 202 e += c * sin(lphi); 203 } 204 205 temp = E / ellipticKComplete( 1.0L - m ); 206 temp *= (atan(t) + mod * PI)/(d * a); 207 temp += e; 208 209 done: 210 if( sign < 0 ) 211 temp = -temp; 212 temp += npio2 * E; 213 return temp; 214 } 215 216 217 /** 218 * Complete elliptic integral of the first kind 219 * 220 * Approximates the integral 221 * 222 * K(m) = $(INTEGRATE 0, π/2) dt/ (sqrt( 1- m $(POWER sin, 2) t)) 223 * 224 * where m = 1 - x, using the approximation 225 * 226 * P(x) - log x Q(x). 227 * 228 * The argument x is used rather than m so that the logarithmic 229 * singularity at x = 1 will be shifted to the origin; this 230 * preserves maximum accuracy. 231 * 232 * x must be in the range 233 * 0 <= x <= 1 234 * 235 * This is equivalent to ellipticF(PI_2, 1-x). 236 * 237 * K(0) = π/2. 238 */ 239 240 real ellipticKComplete(real x) 241 { 242 // verify(x>=0.0L && x<=1.0L); 243 244 static immutable real [] P = [ 245 0x1.62e42fefa39ef35ap+0, // 1.3862943611198906189 246 0x1.8b90bfbe8ed811fcp-4, // 0.096573590279993142323 247 0x1.fa05af797624c586p-6, // 0.030885144578720423267 248 0x1.e979cdfac7249746p-7, // 0.01493761594388688915 249 0x1.1f4cc8890cff803cp-7, // 0.0087676982094322259125 250 0x1.7befb3bb1fa978acp-8, // 0.0057973684116620276454 251 0x1.2c2566aa1d5fe6b8p-8, // 0.0045798659940508010431 252 0x1.7333514e7fe57c98p-8, // 0.0056640695097481470287 253 0x1.09292d1c8621348cp-7, // 0.0080920667906392630755 254 0x1.b89ab5fe793a6062p-8, // 0.0067230886765842542487 255 0x1.28e9c44dc5e26e66p-9, // 0.002265267575136470585 256 0x1.c2c43245d445addap-13, // 0.00021494216542320112406 257 0x1.4ee247035a03e13p-20 // 1.2475397291548388385e-06 258 ]; 259 260 static immutable real [] Q = [ 261 0x1p-1, // 0.5 262 0x1.fffffffffff635eap-4, // 0.12499999999999782631 263 0x1.1fffffff8a2bea1p-4, // 0.070312499993302227507 264 0x1.8ffffe6f40ec2078p-5, // 0.04882812208418620146 265 0x1.323f4dbf7f4d0c2ap-5, // 0.037383701182969303058 266 0x1.efe8a028541b50bp-6, // 0.030267864612427881354 267 0x1.9d58c49718d6617cp-6, // 0.025228683455123323041 268 0x1.4d1a8d2292ff6e2ep-6, // 0.020331037356569904872 269 0x1.b637687027d664aap-7, // 0.013373304362459048444 270 0x1.687a640ae5c71332p-8, // 0.0055004591221382442135 271 0x1.0f9c30a94a1dcb4ep-10, // 0.001036110372590318803 272 0x1.d321746708e92d48p-15 // 5.568631677757315399e-05 273 ]; 274 275 static immutable real LOG4 = 0x1.62e42fefa39ef358p+0; // log(4) 276 277 if( x > real.epsilon ) 278 return poly(x,P) - log(x) * poly(x,Q); 279 if ( x == 0.0L ) 280 return real.infinity; 281 return LOG4 - 0.5L * log(x); 282 } 283 284 /** 285 * Complete elliptic integral of the second kind 286 * 287 * Approximates the integral 288 * 289 * E(m) = $(INTEGRATE 0, π/2) sqrt( 1- m $(POWER sin, 2) t) dt 290 * 291 * where m = 1 - x, using the approximation 292 * 293 * P(x) - x log x Q(x). 294 * 295 * Though there are no singularities, the argument m1 is used 296 * rather than m for compatibility with ellipticKComplete(). 297 * 298 * E(1) = 1; E(0) = π/2. 299 * m must be in the range 0 <= m <= 1. 300 */ 301 302 real ellipticEComplete(real x) 303 { 304 verify(x>=0 && x<=1.0); 305 static immutable real [] P = [ 306 0x1.c5c85fdf473f78f2p-2, // 0.44314718055994670505 307 0x1.d1591f9e9a66477p-5, // 0.056805192715569305834 308 0x1.65af6a7a61f587cp-6, // 0.021831373198011179718 309 0x1.7a4d48ed00d5745ap-7, // 0.011544857605264509506 310 0x1.d4f5fe4f93b60688p-8, // 0.0071557756305783152481 311 0x1.4cb71c73bac8656ap-8, // 0.0050768322432573952962 312 0x1.4a9167859a1d0312p-8, // 0.0050440671671840438539 313 0x1.dd296daa7b1f5b7ap-8, // 0.0072809117068399675418 314 0x1.04f2c29224ba99b6p-7, // 0.0079635095646944542686 315 0x1.0f5820e2d80194d8p-8, // 0.0041403847015715420009 316 0x1.95ee634752ca69b6p-11, // 0.00077425232385887751162 317 0x1.0c58aa9ab404f4fp-15 // 3.1989378120323412946e-05 318 ]; 319 320 static immutable real [] Q = [ 321 0x1.ffffffffffffb1cep-3, // 0.24999999999999986434 322 0x1.7ffffffff29eaa0cp-4, // 0.093749999999239422678 323 0x1.dfffffbd51eb098p-5, // 0.058593749514839092674 324 0x1.5dffd791cb834c92p-5, // 0.04272453406734691973 325 0x1.1397b63c2f09a8ep-5, // 0.033641677787700181541 326 0x1.c567cde5931e75bcp-6, // 0.02767367465121309044 327 0x1.75e0cae852be9ddcp-6, // 0.022819708015315777007 328 0x1.12bb968236d4e434p-6, // 0.016768357258894633433 329 0x1.1f6572c1c402d07cp-7, // 0.0087706384979640787504 330 0x1.452c6909f88b8306p-9, // 0.0024808767529843311337 331 0x1.1f7504e72d664054p-12, // 0.00027414045912208516032 332 0x1.ad17054dc46913e2p-18 // 6.3939381343012054851e-06 333 ]; 334 if (x==0) 335 return 1.0L; 336 return 1.0L + x * poly(x,P) - log(x) * (x * poly(x,Q) ); 337 } 338 339 unittest { 340 test( ellipticF(1, 0)==1); 341 test(ellipticEComplete(0)==1); 342 test(ellipticEComplete(1)==PI_2); 343 test(feqrel(ellipticKComplete(1),PI_2)>= real.mant_dig-1); 344 test(ellipticKComplete(0)==real.infinity); 345 // assert(ellipticKComplete(1)==0); //-real.infinity); 346 347 real x=0.5653L; 348 test(ellipticKComplete(1-x) == ellipticF(PI_2, x) ); 349 test(ellipticEComplete(1-x) == ellipticE(PI_2, x) ); 350 } 351 352 /** 353 * Incomplete elliptic integral of the third kind 354 * 355 * Approximates the integral 356 * 357 * PI(n; phi | m) = $(INTEGRATE t=0, phi) dt/((1 - n $(POWER sin,2)t) * sqrt( 1- m $(POWER sin, 2) t)) 358 * 359 * of amplitude phi, modulus m, and characteristic n using Gauss-Legendre 360 * quadrature. 361 * 362 * Note that ellipticPi(PI_2, m, 1) is infinite for any m. 363 */ 364 real ellipticPi(real phi, real m, real n) 365 { 366 // BUGS: This implementation suffers from poor precision. 367 static immutable double [] t = [ 368 0.9931285991850949, 0.9639719272779138, 369 0.9122344282513259, 0.8391169718222188, 370 0.7463319064601508, 0.6360536807265150, 371 0.5108670019508271, 0.3737060887154195, 372 0.2277858511416451, 0.7652652113349734e-1 373 ]; 374 static immutable double [] w =[ 375 0.1761400713915212e-1, 0.4060142980038694e-1, 376 0.6267204833410907e-1, 0.8327674157670475e-1, 377 0.1019301198172404, 0.1181945319615184, 378 0.1316886384491766, 0.1420961093183820, 379 0.1491729864726037, 0.1527533871307258 380 ]; 381 bool b1 = (m==1) && abs(phi-90)<=1e-8; 382 bool b2 = (n==1) && abs(phi-90)<=1e-8; 383 if (b1 || b2) return real.infinity; 384 real c1 = 0.87266462599716e-2 * phi; 385 real c2 = c1; 386 double x = 0; 387 for (int i=0; i< t.length; ++i) { 388 real c0 = c2 * t[i]; 389 real t1 = c1 + c0; 390 real t2 = c1 - c0; 391 real s1 = sin(t1); // sin(c1 * (1 + t[i])) 392 real s2 = sin(t2); // sin(c1 * (1 - t[i])) 393 real f1 = 1.0 / ((1.0 - n * s1 * s1) * sqrt(1.0 - m * s1 * s1)); 394 real f2 = 1.0 / ((1.0 - n * s2 * s2) * sqrt(1.0 - m * s2 * s2)); 395 x+= w[i]*(f1+f2); 396 } 397 return c1 * x; 398 } 399 400 /** 401 * Complete elliptic integral of the third kind 402 */ 403 real ellipticPiComplete(real m, real n) 404 { 405 verify(m>=-1.0 && m<=1.0); 406 return ellipticPi(PI_2, m, n); 407 }