1 /*
2 * $Id: elliptic.i,v 1.2 2007-11-10 20:03:49 dhmunro Exp $
3 * elliptic functions
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 /* Abramowitz and Stegun, sections 16.4, 17.2.17-19, 17.6 */
12
13 local elliptic;
14 /* DOCUMENT elliptic, ell_am, ell_f, ell_e, dn_, ellip_k, ellip_e
15
16 The elliptic integral of the first kind is:
17
18 u = integral[0 to phi]( dt / sqrt(1-m*sin(t)^2) )
19
20 The functions ell_f and ell_am compute this integral and its
21 inverse:
22
23 u = ell_f(phi, m)
24 phi = ell_am(u, m)
25
26 The Jacobian elliptic functions can be computed from the
27 "amplitude" ell_am by means of:
28
29 sn(u|m) = sin(ell_am(u,m))
30 cn(u|m) = cos(ell_am(u,m))
31 dn(u|m) = dn_(ell_am(u,m)) = sqrt(1-m*sn(u|m)^2)
32
33 The other nine functions are sc=sn/cn, cs=cn/sn, nd=1/dn,
34 cd=cn/dn, dc=dn/cn, ns=1/sn, sd=sn/dn, nc=1/cn, and ds=dn/sn.
35 (The notation u|m does not means yorick's | operator; it is
36 the mathematical notation, not valid yorick code!)
37
38 The parameter M is given in three different notations:
39 as M, the "parameter",
40 as k, the "modulus", or
41 as alpha, the "modular angle",
42 which are related by: M = k^2 = sin(alpha)^2. The yorick elliptic
43 functions in terms of M may need to be written
44 ell_am(u,k^2) or ell_am(u,sin(alpha)^2)
45 in order to agree with the definitions in other references.
46 Sections 17.2.17-19 of Abramowitz and Stegun explains these notations,
47 and chapters 16 and 17 present a compact overview of the subject of
48 elliptic functions in general.
49
50 The parameter M must be a scalar; U may be an array. The
51 exceptions are the complete elliptic integrals ellip_k and
52 ellip_e which accept an array of M values.
53
54 The ell_am function uses the external variable ell_m if M is
55 omitted, otherwise stores M in ell_m. Hence, you may set ell_m,
56 then simply call ell_am(u) if you have a series of calls with
57 the same value of M; this also allows the dn_ function to work
58 without a second specification of M.
59
60 The elliptic integral of the second kind is:
61
62 u = integral[0 to phi]( dt * sqrt(1-m*sin(t)^2) )
63
64 The function ell_e computes this integral:
65
66 u = ell_e(phi, m)
67
68 The special values ell_f(pi/2,m) and ell_e(pi/2,m) are the complete
69 elliptic integrals of the first and second kinds; separate functions
70 ellip_k and ellip_e are provided to compute them.
71
72 Note that the function ellip_k is infinite for M=1 and for large
73 negative M. The "natural" range for M is 0<=M<=1; all other real
74 values can be "reduced" to this range by various transformations;
75 the logarithmic singularity of ellip_k is actually very mild, and
76 other functions such as ell_am are perfectly well-defined there.
77
78 Here are the sum formulas for elliptic functions:
79
80 sn(u+v) = ( sn(u)*cn(v)*dn(v) + sn(v)*cn(u)*dn(u) ) /
81 ( 1 - m*sn(u)^2*sn(v)^2 )
82 cn(u+v) = ( cn(u)*cn(v) - sn(u)*dn(u)*sn(v)*dn(v) ) /
83 ( 1 - m*sn(u)^2*sn(v)^2 )
84 dn(u+v) = ( dn(u)*dn(v) - m*sn(u)*cn(u)*sn(v)*cn(v) ) /
85 ( 1 - m*sn(u)^2*sn(v)^2 )
86
87 And the formulas for pure imaginary values:
88
89 sn(1i*u,m) = 1i * sc(u,1-m)
90 cn(1i*u,m) = nc(u,1-m)
91 dn(1i*u,m) = dc(u,1-m)
92
93 SEE ALSO: ell_am, ell_f, ell_e, dn_, ellip_k, ellip_e
94 */
95
ell_am(u,m)96 func ell_am(u,m)
97 /* DOCUMENT ell_am(u)
98 or ell_am(u,m)
99
100 returns the "amplitude" (an angle in radians) for the Jacobi
101 elliptic functions at U, with parameter M. That is,
102 phi = ell_am(u,m)
103 means that
104 u = integral[0 to phi]( dt / sqrt(1-m*sin(t)^2) )
105
106 Thus ell_am is the inverse of the incomplete elliptic function
107 of the first kind ell_f. See help,elliptic for more.
108
109 SEE ALSO: elliptic
110 */
111 {
112 if (structof(u)==complex || structof(m)==complex)
113 error, "elliptic function not valid for complex";
114 /* set up the arithmetic-geometric mean scale */
115 extern ell_m, _agm_m, _agm_n, _agm_coa, _agm_a, _agm_sn;
116 if (is_void(m)) m= ell_m;
117 if (m != ell_m) {
118 ell_m= m= double(m);
119 if (m<0.) {
120 _agm_sn= 1./(1.-m);
121 m*= -_agm_sn;
122 _agm_sn= sqrt(_agm_sn);
123 } else if (m>1.) {
124 m= 1./m;
125 _agm_sn= sqrt(m);
126 }
127 _agm_m= m;
128 _agm_coa()= 0.;
129 _agm_n= 0;
130 _agm_a= 1.;
131 if (m!=1.) {
132 b= sqrt(1.-m);
133 for (;;) { /* maximum of 8 passes for 64-bit double */
134 c= 0.5*(_agm_a-b);
135 if (!c) break;
136 am= _agm_a-c; /* arithmetic mean */
137 _agm_coa(++_agm_n)= c/am;
138 if (am==_agm_a) break;
139 b= sqrt(_agm_a*b); /* geometric mean */
140 _agm_a= am;
141 }
142 _agm_a*= 2.^_agm_n;
143 }
144 }
145
146 phi= _agm_a*u;
147 if (m<0. || m>1.) phi/= _agm_sn;
148 if (m!=1.)
149 for (n=_agm_n ; n>0 ; n--) phi= 0.5*(phi+asin(_agm_coa(n)*sin(phi)));
150 else
151 phi= atan(tanh(phi), sech(phi));
152 if (m<0.) {
153 cn= cos(phi);
154 phi= sin(phi);
155 nd= 1./sqrt(1.-m*phi*phi);
156 phi= atan(_agm_sn*phi*nd, cn*nd);
157 } else if (m>1.) {
158 phi= sin(phi);
159 phi= atan(_agm_sn*phi, sqrt(1.-_agm_m*phi*phi));
160 }
161 return phi;
162 }
163
164 _agm_coa= array(0.,16);
165 _agm_n= 0;
166 _agm_m= 0.0;
167 _agm_a= 1.0;
168 _agm_sn= 1.0;
169 ell_m= 0.0;
170
dn_(phi,m)171 func dn_(phi, m)
172 /* DOCUMENT dn_(ell_am(u,m))
173
174 return the Jacobian elliptic function dn(u|m). The external
175 variable ell_m must be set properly before calling dn_.
176
177 SEE ALSO: elliptic, ell_am
178 */
179 {
180 if (is_void(m)) m= ell_m;
181 phi= sin(phi);
182 return sqrt(1.-m*phi*phi);
183 }
184
ell_f(phi,m)185 func ell_f(phi,m)
186 /* DOCUMENT ell_f(phi,m)
187
188 returns the incomplete elliptic integral of the first kind F(phi|M).
189 That is,
190 u = ell_f(phi,m)
191 means that
192 u = integral[0 to phi]( dt / sqrt(1-m*sin(t)^2) )
193
194 See help,elliptic for more.
195
196 SEE ALSO: elliptic, ell_e
197 */
198 {
199 if (structof(phi)==complex || structof(m)==complex)
200 error, "elliptic integral not valid for complex";
201 orig_m= m= double(m);
202 if (m>1.) {
203 scale= sqrt(m);
204 phi= asin(scale*sin(phi));
205 m= 1./m;
206 } else if (m<0.) {
207 scale= sqrt(1.-m);
208 phi= 0.5*pi - phi;
209 m/= m-1.;
210 }
211
212 if (m==1.) {
213 phi= atanh(sin(phi));
214 a= 0.;
215 } else { /* compute using arithmetic-geometric mean */
216 sgn= sign(phi);
217 phi= abs(phi);
218 a= 1.;
219 b= sqrt(1.-m);
220 twon= 1.0;
221 pi2= 0.5*pi;
222 for (;;) { /* maximum of 8 passes for 64-bit double */
223 c= 0.5*(a-b);
224 if (!c) break;
225 phase= (phi+pi2)/pi;
226 cycle= floor(phase);
227 phi*= 1. + 1.e-15*(cycle==phase);
228 phi+= atan((b/a)*tan(phi)) + pi*cycle;
229 twon*= 2.0;
230 am= a-c;
231 if (am==a) break;
232 b= sqrt(a*b);
233 a= am;
234 }
235 phi/= twon*a*sgn;
236 }
237
238 if (orig_m>1.) {
239 phi/= scale;
240 } else if (orig_m<0.) {
241 phi= (0.5*pi/a - phi)/scale;
242 }
243
244 return phi;
245 }
246
ell_e(phi,m)247 func ell_e(phi,m)
248 /* DOCUMENT ell_e(phi,m)
249
250 returns the incomplete elliptic integral of the second kind E(phi|M).
251 That is,
252 u = ell_e(phi,m)
253 means that
254 u = integral[0 to phi]( dt * sqrt(1-m*sin(t)^2) )
255
256 See help,elliptic for more.
257
258 SEE ALSO: elliptic, ell_f
259 */
260 {
261 if (structof(phi)==complex || structof(m)==complex)
262 error, "elliptic integral not valid for complex";
263 orig_m= m= double(m);
264 if (m>1.) {
265 scale= sqrt(m);
266 phi= asin(scale*sin(phi));
267 m= 1./m;
268 } else if (m<0.) {
269 scale= sqrt(1.-m);
270 phi= 0.5*pi - phi;
271 m/= m-1.;
272 }
273
274 if (m==1.) {
275 per= floor((phi+0.5*pi)/pi);
276 phi= sin(phi)*sign(0.5-abs(per%2.)) + 2.*per;
277 a= 0.;
278 } else { /* compute using arithmetic-geometric mean */
279 sgn= sign(phi);
280 phi= abs(phi);
281 a= 1.;
282 b= sqrt(1.-m);
283 eok= 1.-0.5*m;
284 cs= 0.;
285 twon= 1.0;
286 pi2= 0.5*pi;
287 for (;;) { /* maximum of 8 passes for 64-bit double */
288 c= 0.5*(a-b);
289 if (!c) break;
290 phi+= atan((b/a)*tan(phi)) + pi*floor((phi+pi2)/pi);
291 cs+= c*sin(phi);
292 eok-= twon*c*c;
293 twon*= 2.0;
294 am= a-c;
295 if (am==a) break;
296 b= sqrt(a*b);
297 a= am;
298 }
299 f= sgn*phi/(twon*a);
300 phi= eok*f + sgn*cs;
301 }
302
303 if (orig_m>1.) {
304 phi= (phi-(1.-m)*f)/scale;
305 } else if (orig_m<0.) {
306 phi= (eok*0.5*pi/a - phi)*scale;
307 }
308
309 return phi;
310 }
311
ellip_k(m)312 func ellip_k(m)
313 /* DOCUMENT ellip_k(m)
314
315 returns the complete elliptic integral of the first kind K(M):
316 K(M) = integral[0 to pi/2]( dt / sqrt(1-M*sin(t)^2) )
317
318 See help,elliptic for more.
319
320 SEE ALSO: elliptic, ellip_e, ell_f
321 */
322 {
323 if (anyof(m>=1.)) error, "ellip_k(m) not computed for m>=1";
324
325 m= double(m);
326 mask= (m>=0.);
327 list= where(mask);
328 if (numberof(list)) scale= array(0.5*pi,numberof(list));
329 list= where(!mask);
330 if (numberof(list)) {
331 sm= 1./(1.-m(list));
332 m(list)*= -sm;
333 sm= 0.5*pi*sqrt(sm);
334 }
335 scale= merge(scale,sm,mask);
336
337 a= array(1.,dimsof(m));
338 b= sqrt(1.-m);
339 for (;;) {
340 c= 0.5*(a-b);
341 am= a-c;
342 if (allof(am==a)) break;
343 b= sqrt(a*b);
344 a= am;
345 }
346
347 return scale/a;
348 }
349
ellip_e(m)350 func ellip_e(m)
351 /* DOCUMENT ellip_e(m)
352
353 returns the complete elliptic integral of the second kind E(M):
354 E(M) = integral[0 to pi/2]( dt * sqrt(1-M*sin(t)^2) )
355
356 See help,elliptic for more.
357
358 SEE ALSO: elliptic, ellip_k, ell_e
359 */
360 {
361 if (anyof(m>1.)) error, "ellip_e(m) not computed for m>1";
362
363 m= double(m);
364 mask= (m>=0.);
365 list= where(mask);
366 if (numberof(list)) scale= array(0.5*pi,numberof(list));
367 list= where(!mask);
368 if (numberof(list)) {
369 sm= 1.-m(list);
370 m(list)/= -sm;
371 sm= 0.5*pi*sqrt(sm);
372 }
373 scale= merge(scale,sm,mask);
374
375 mask= (m!=1.);
376 list= where(mask);
377 if (numberof(list)) {
378 m= m(list);
379 a= array(1.,numberof(m));
380 b= sqrt(1.-m);
381 e= 1.-0.5*m;
382 twon= 1.0;
383 for (n=0 ;;) { /* maximum of 8 passes for 64-bit double */
384 c= 0.5*(a-b);
385 am= a-c;
386 if (allof(am==a)) break;
387 e-= twon*c*c;
388 twon*= 2.0;
389 b= sqrt(a*b);
390 a= am;
391 }
392 e/= a;
393 }
394 list= where(!mask);
395 if (numberof(list)) em= array(2./pi,numberof(list));
396
397 return scale*merge(e,em,mask);
398 }
399
ellip2_k(m)400 func ellip2_k(m)
401 /* DOCUMENT ellip2_k(m)
402
403 returns the complete elliptic integral of the first kind K(M):
404 K(M) = integral[0 to pi/2]( dt / sqrt(1-M*sin(t)^2) )
405 accurate to 2e-8 for 0<=M<1
406
407 SEE ALSO: elliptic, ellip_e, ell_f
408 */
409 {
410 m= 1.-m;
411 return poly(m,1.38629436112,0.09666344259,0.03590092383,
412 0.03742563713,0.01451196212) +
413 log(1./m)*poly(m,0.5,0.12498593597,0.06880248576,
414 0.03328355346,0.00441787012);
415 }
416
ellip2_e(m)417 func ellip2_e(m)
418 /* DOCUMENT ellip2_e(m)
419
420 returns the complete elliptic integral of the second kind E(M):
421 E(M) = integral[0 to pi/2]( dt * sqrt(1-M*sin(t)^2) )
422 accurate to 2e-8 for 0<=M<=1
423
424 SEE ALSO: elliptic, ellip_k, ell_e
425 */
426 {
427 m= 1.-m;
428 return poly(m,1.,0.44325141463,0.06260601220,
429 0.04757383546,0.01736506451) +
430 log(1./m)*poly(m,0.,0.24998368310,0.09200180037,
431 0.04069697526,0.00526449639);
432 }
433