6.11 Elliptic Integrals and Jacobian Elliptic Functions 267 #include <math.h> #include "nrutil.h" #define ERRTOL 0.04 #define TINY 1.69e-38 #define SQRTNY 1.3e-19 #define BIG 3.e37 #define TNBG (TINYBIG) #define COMP1 (2.236/SQRTNY) #define COMP2 (TNBG*TNBG/25.0) #define THIRD (1.0/3.0) #define C1 0.3 #define C2(1.0/7.0) http://www.nr. #define C3 0.375 #define C4(9.0/22.0) 83 float rc(float x,float y) 鱼 granted for 19881992 Computes Carlson's degenerate elliptic integral,Rc(x,y).x must be nonnegative and y must be nonzero.If y 0,the Cauchy principal value is returned.TINY must be at least 5 times 11800 the machine underflow limit,BIG at most one fifth the machine maximum overflow limit. float alamb,ave,s,w,xt,yt; Cambridge from NUMERICAL RECIPESI if (x 0.0 II y =0.0 II (x+fabs(y))<TINY II (x+fabs(y))>BIG II (y<-COP122x>0.0&x<C0MP2)) nrerror("invalid arguments in rc"); if(y>0.0) xt=x; yt=y; server computer, (North America to make one paper UnN电.t THE w=1.0; ART 】e1se{ xt=x-yi 是 yt =-y; Programs w=sqrt(x)/sqrt(xt); send do alamb=2.0*sgrt(xt)*sqrt(yt)+yt; email st st Copyright (C) xt=0.25*(xt+alamb); yt=0.25*(yt+alamb); ave=THIRD*(xt+yt+yt); s=(yt-ave)/ave; OF SCIENTIFIC COMPUTING(ISBN while (fabs(s)>ERRTOL); return w*(1.0+s*s*(C1+s*(C2+s*(C3+s*C4))))/sgrt(ave); v@cam 10-621 At times you may want to express your answer in Legendre's notation.Alter- natively,you may be given results in that notation and need to compute their values 1988-1992 by Numerical Recipes 43108 with the programs given above.It is a simple matter to transform back and forth. The Legendre elliptic integral of the Ist kind is defined as (outside do North Software. F(中，k)≡ (6.11.17) /1-k2sin20 The complete elliptic integral of the Ist kind is given by K(k)≡F(π/2，) (6.11.18) In terms of RF. F(,k)=sinoRg(cos2p,1-k2sin2中，1) (6.11.19) K(k)=RF(0,1-k2,1)6.11 Elliptic Integrals and Jacobian Elliptic Functions 267 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). #include <math.h> #include "nrutil.h" #define ERRTOL 0.04 #define TINY 1.69e-38 #define SQRTNY 1.3e-19 #define BIG 3.e37 #define TNBG (TINY*BIG) #define COMP1 (2.236/SQRTNY) #define COMP2 (TNBG*TNBG/25.0) #define THIRD (1.0/3.0) #define C1 0.3 #define C2 (1.0/7.0) #define C3 0.375 #define C4 (9.0/22.0) float rc(float x, float y) Computes Carlson’s degenerate elliptic integral, RC(x, y). x must be nonnegative and y must be nonzero. If y < 0, the Cauchy principal value is returned. TINY must be at least 5 times the machine underflow limit, BIG at most one fifth the machine maximum overflow limit. { float alamb,ave,s,w,xt,yt; if (x < 0.0 || y == 0.0 || (x+fabs(y)) < TINY || (x+fabs(y)) > BIG || (y<-COMP1 && x > 0.0 && x < COMP2)) nrerror("invalid arguments in rc"); if (y > 0.0) { xt=x; yt=y; w=1.0; } else { xt=x-y; yt = -y; w=sqrt(x)/sqrt(xt); } do { alamb=2.0*sqrt(xt)*sqrt(yt)+yt; xt=0.25*(xt+alamb); yt=0.25*(yt+alamb); ave=THIRD*(xt+yt+yt); s=(yt-ave)/ave; } while (fabs(s) > ERRTOL); return w*(1.0+s*s*(C1+s*(C2+s*(C3+s*C4))))/sqrt(ave); } At times you may want to express your answer in Legendre’s notation. Alter￾natively, you may be given results in that notation and need to compute their values with the programs given above. It is a simple matter to transform back and forth. The Legendre elliptic integral of the 1st kind is defined as F(φ, k) ≡
φ 0 dθ 1 − k2 sin2 θ (6.11.17) The complete elliptic integral of the 1st kind is given by K(k) ≡ F(π/2, k) (6.11.18) In terms of RF , F(φ, k) = sin φRF (cos2 φ, 1 − k2 sin2 φ, 1) K(k) = RF (0, 1 − k2, 1) (6.11.19)