  GNU-Darwin Web
i/elliptic Home Manual Packages Global Index Keywords Quick Reference ``` /* * elliptic.i -- \$Id\$ * elliptic functions * * Copyright (c) 1999. See accompanying LEGAL file for details. */ /* Abramowitz and Stegun, sections 16.4, 17.2.17-19, 17.6 */ local elliptic ; /* DOCUMENT elliptic, ell_am, ell_f, ell_e, dn_, ellip_k, ellip_e The elliptic integral of the first kind is: u = integral[0 to phi]( dt / sqrt(1-m*sin(t)^2) ) The functions ell_f and ell_am compute this integral and its inverse: u = ell_f(phi, m) phi = ell_am(u, m) The Jacobian elliptic functions can be computed from the "amplitude" ell_am by means of: sn(u|m) = sin(ell_am(u,m)) cn(u|m) = cos(ell_am(u,m)) dn(u|m) = dn_(ell_am(u,m)) = sqrt(1-m*sn(u|m)^2) The other nine functions are sc=sn/cn, cs=cn/sn, nd=1/dn, cd=cn/dn, dc=dn/cn, ns=1/sn, sd=sn/dn, nc=1/cn, and ds=dn/sn. (The notation u|m does not means yorick's | operator; it is the mathematical notation, not valid yorick code!) The parameter M is given in three different notations: as M, the "parameter", as k, the "modulus", or as alpha, the "modular angle", which are related by: M = k^2 = sin(alpha)^2. The yorick elliptic functions in terms of M may need to be written ell_am(u,k^2) or ell_am(u,sin(alpha)^2) in order to agree with the definitions in other references. Sections 17.2.17-19 of Abramowitz and Stegun explains these notations, and chapters 16 and 17 present a compact overview of the subject of elliptic functions in general. The parameter M must be a scalar; U may be an array. The exceptions are the complete elliptic integrals ellip_k and ellip_e which accept an array of M values. The ell_am function uses the external variable ell_m if M is omitted, otherwise stores M in ell_m. Hence, you may set ell_m, then simply call ell_am(u) if you have a series of calls with the same value of M; this also allows the dn_ function to work without a second specification of M. The elliptic integral of the second kind is: u = integral[0 to phi]( dt * sqrt(1-m*sin(t)^2) ) The function ell_e computes this integral: u = ell_e(phi, m) The special values ell_f(pi/2,m) and ell_e(pi/2,m) are the complete elliptic integrals of the first and second kinds; separate functions ellip_k and ellip_e are provided to compute them. Note that the function ellip_k is infinite for M=1 and for large negative M. The "natural" range for M is 0<=M<=1; all other real values can be "reduced" to this range by various transformations; the logarithmic singularity of ellip_k is actually very mild, and other functions such as ell_am are perfectly well-defined there. Here are the sum formulas for elliptic functions: sn(u+v) = ( sn(u)*cn(v)*dn(v) + sn(v)*cn(u)*dn(u) ) / ( 1 - m*sn(u)^2*sn(v)^2 ) cn(u+v) = ( cn(u)*cn(v) - sn(u)*dn(u)*sn(v)*dn(v) ) / ( 1 - m*sn(u)^2*sn(v)^2 ) dn(u+v) = ( dn(u)*dn(v) - m*sn(u)*cn(u)*sn(v)*cn(v) ) / ( 1 - m*sn(u)^2*sn(v)^2 ) And the formulas for pure imaginary values: sn(1i*u,m) = 1i * sc(u,1-m) cn(1i*u,m) = nc(u,1-m) dn(1i*u,m) = dc(u,1-m) SEE ALSO: ell_am, ell_f, ell_e, dn_, ellip_k, ellip_e */ func ell_am (u,m) /* DOCUMENT ell_am(u) or ell_am(u,m) returns the "amplitude" (an angle in radians) for the Jacobi elliptic functions at U, with parameter M. That is, phi = ell_am(u,m) means that u = integral[0 to phi]( dt / sqrt(1-m*sin(t)^2) ) Thus ell_am is the inverse of the incomplete elliptic function of the first kind ell_f. See help,elliptic for more. SEE ALSO: elliptic */ { /* set up the arithmetic-geometric mean scale */ extern ell_m, _agm_m, _agm_n, _agm_coa, _agm_a, _agm_sn; if (is_void(m)) m= ell_m; if (m != ell_m) { ell_m= m= double(m); if (m<0.) { _agm_sn= 1./(1.-m); m*= -_agm_sn; _agm_sn= sqrt(_agm_sn); } else if (m>1.) { m= 1./m; _agm_sn= sqrt(m); } _agm_m= m; _agm_coa()= 0.; _agm_n= 0; _agm_a= 1.; if (m!=1.) { b= sqrt(1.-m); for (;;) { /* maximum of 8 passes for 64-bit double */ c= 0.5*(_agm_a-b); if (!c) break; am= _agm_a-c; /* arithmetic mean */ _agm_coa(++_agm_n)= c/am; if (am==_agm_a) break; b= sqrt(_agm_a*b); /* geometric mean */ _agm_a= am; } _agm_a*= 2.^_agm_n; } } phi= _agm_a*u; if (m<0. || m>1.) phi/= _agm_sn; if (m!=1.) for (n=_agm_n ; n>0 ; n--) phi= 0.5*(phi+asin(_agm_coa(n)*sin(phi))); else phi= atan(tanh(phi), sech(phi)); if (m<0.) { cn= cos(phi); phi= sin(phi); nd= 1./sqrt(1.-m*phi*phi); phi= atan(_agm_sn*phi*nd, cn*nd); } else if (m>1.) { phi= sin(phi); phi= atan(_agm_sn*phi, sqrt(1.-_agm_m*phi*phi)); } return phi; } _agm_coa= array(0.,16); _agm_n= 0; _agm_m= 0.0; _agm_a= 1.0; _agm_sn= 1.0; ell_m= 0.0; func dn_ (phi, m) /* DOCUMENT dn_(ell_am(u,m)) return the Jacobian elliptic function dn(u|m). The external variable ell_m must be set properly before calling dn_. SEE ALSO: elliptic, ell_am */ { if (is_void(m)) m= ell_m; phi= sin(phi); return sqrt(1.-m*phi*phi); } func ell_f (phi,m) /* DOCUMENT ell_f(phi,m) returns the incomplete elliptic integral of the first kind F(phi|M). That is, u = ell_f(phi,m) means that u = integral[0 to phi]( dt / sqrt(1-m*sin(t)^2) ) See help,elliptic for more. SEE ALSO: elliptic, ell_e */ { orig_m= m= double(m); if (m>1.) { scale= sqrt(m); phi= asin(scale*sin(phi)); m= 1./m; } else if (m<0.) { scale= sqrt(1.-m); phi= 0.5*pi - phi; m/= m-1.; } if (m==1.) { phi= atanh(sin(phi)); a= 0.; } else { /* compute using arithmetic-geometric mean */ sgn= sign(phi); phi= abs(phi); a= 1.; b= sqrt(1.-m); twon= 1.0; pi2= 0.5*pi; for (;;) { /* maximum of 8 passes for 64-bit double */ c= 0.5*(a-b); if (!c) break; phase= (phi+pi2)/pi; cycle= floor(phase); phi*= 1. + 1.e-15*(cycle==phase); phi+= atan((b/a)*tan(phi)) + pi*cycle; twon*= 2.0; am= a-c; if (am==a) break; b= sqrt(a*b); a= am; } phi/= twon*a*sgn; } if (orig_m>1.) { phi/= scale; } else if (orig_m<0.) { phi= (0.5*pi/a - phi)/scale; } return phi; } func ell_e (phi,m) /* DOCUMENT ell_e(phi,m) returns the incomplete elliptic integral of the second kind E(phi|M). That is, u = ell_e(phi,m) means that u = integral[0 to phi]( dt * sqrt(1-m*sin(t)^2) ) See help,elliptic for more. SEE ALSO: elliptic, ell_f */ { orig_m= m= double(m); if (m>1.) { scale= sqrt(m); phi= asin(scale*sin(phi)); m= 1./m; } else if (m<0.) { scale= sqrt(1.-m); phi= 0.5*pi - phi; m/= m-1.; } if (m==1.) { per= floor((phi+0.5*pi)/pi); phi= sin(phi)*sign(0.5-abs(per%2.)) + 2.*per; a= 0.; } else { /* compute using arithmetic-geometric mean */ sgn= sign(phi); phi= abs(phi); a= 1.; b= sqrt(1.-m); eok= 1.-0.5*m; cs= 0.; twon= 1.0; pi2= 0.5*pi; for (;;) { /* maximum of 8 passes for 64-bit double */ c= 0.5*(a-b); if (!c) break; phi+= atan((b/a)*tan(phi)) + pi*floor((phi+pi2)/pi); cs+= c*sin(phi); eok-= twon*c*c; twon*= 2.0; am= a-c; if (am==a) break; b= sqrt(a*b); a= am; } f= sgn*phi/(twon*a); phi= eok*f + sgn*cs; } if (orig_m>1.) { phi= (phi-(1.-m)*f)/scale; } else if (orig_m<0.) { phi= (eok*0.5*pi/a - phi)*scale; } return phi; } func ellip_k (m) /* DOCUMENT ellip_k(m) returns the complete elliptic integral of the first kind K(M): K(M) = integral[0 to pi/2]( dt / sqrt(1-M*sin(t)^2) ) See help,elliptic for more. SEE ALSO: elliptic, ellip_e, ell_f */ { if (anyof(m>=1.)) error, "ellip_k(m) not computed for m>=1"; m= double(m); mask= (m>=0.); list= where(mask); if (numberof(list)) scale= array(0.5*pi,numberof(list)); list= where(!mask); if (numberof(list)) { sm= 1./(1.-m(list)); m(list)*= -sm; sm= 0.5*pi*sqrt(sm); } scale= merge(scale,sm,mask); a= array(1.,dimsof(m)); b= sqrt(1.-m); for (;;) { c= 0.5*(a-b); am= a-c; if (allof(am==a)) break; b= sqrt(a*b); a= am; } return scale/a; } func ellip_e (m) /* DOCUMENT ellip_e(m) returns the complete elliptic integral of the second kind E(M): E(M) = integral[0 to pi/2]( dt * sqrt(1-M*sin(t)^2) ) See help,elliptic for more. SEE ALSO: elliptic, ellip_k, ell_e */ { if (anyof(m>1.)) error, "ellip_e(m) not computed for m>1"; m= double(m); mask= (m>=0.); list= where(mask); if (numberof(list)) scale= array(0.5*pi,numberof(list)); list= where(!mask); if (numberof(list)) { sm= 1.-m(list); m(list)/= -sm; sm= 0.5*pi*sqrt(sm); } scale= merge(scale,sm,mask); mask= (m!=1.); list= where(mask); if (numberof(list)) { m= m(list); a= array(1.,numberof(m)); b= sqrt(1.-m); e= 1.-0.5*m; twon= 1.0; for (n=0 ;;) { /* maximum of 8 passes for 64-bit double */ c= 0.5*(a-b); am= a-c; if (allof(am==a)) break; e-= twon*c*c; twon*= 2.0; b= sqrt(a*b); a= am; } e/= a; } list= where(!mask); if (numberof(list)) em= array(2./pi,numberof(list)); return scale*merge(e,em,mask); } func ellip2_k (m) /* DOCUMENT ellip2_k(m) returns the complete elliptic integral of the first kind K(M): K(M) = integral[0 to pi/2]( dt / sqrt(1-M*sin(t)^2) ) accurate to 2e-8 for 0<=M<1 SEE ALSO: elliptic, ellip_e, ell_f */ { m= 1.-m; return poly(m,1.38629436112,0.09666344259,0.03590092383, 0.03742563713,0.01451196212) + log(1./m)*poly(m,0.5,0.12498593597,0.06880248576, 0.03328355346,0.00441787012); } func ellip2_e (m) /* DOCUMENT ellip2_e(m) returns the complete elliptic integral of the second kind E(M): E(M) = integral[0 to pi/2]( dt * sqrt(1-M*sin(t)^2) ) accurate to 2e-8 for 0<=M<=1 SEE ALSO: elliptic, ellip_k, ell_e */ { m= 1.-m; return poly(m,1.,0.44325141463,0.06260601220, 0.04757383546,0.01736506451) + log(1./m)*poly(m,0.,0.24998368310,0.09200180037, 0.04069697526,0.00526449639); } ```