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