6.11 Elliptic Integrals and Jacobian Elliptic Functions 265 #define TINY 1.0e-25 #define BIG 4.5e21 #define C1 (3.0/14.0) #define C2 (1.0/6.0) 姓daf1na3 (9.0/22.0) #def1nec4(3.0/26.0) #define C5 (0.25*C3) #define C6 (1.5*C4) float rd(float x,float y,float z) Computes Carlson's elliptic integral of the second kind,Rp(z,y,2).x and y must be non- negative,and at most one can be zero.2 must be positive.TINY must be at least twice the http://www.n Permission is read able files negative 2/3 power of the machine overflow limit.BIG must be at most 0.1 x ERRTOL times the negative 2/3 power of the machine underflow limit. 83g float alamb,ave,delx,dely,delz,ea,eb,ec,ed,ee,fac,sqrtx,sqrty (including this one) granted fori sqrtz,sum,xt,yt,zt; if (FMIN(x,y)<0.0 II FMIN(x+y,z)<TINY II FMAX(FMAX(x,y),z)>BIG) 11-800-872 internet nrerror("invalid arguments in rd"); xt=x; yt-y; from NUMERICAL RECIPES IN 2t=2 sum=0.0; fac=1.0; do 7423 (North America sqrtx=sgrt(xt); sqrty=sqrt(yt); sqrtz=sqrt(zt); alamb=sgrtx*(sqrty+sqrtz)+sqrty.sqrtz; to any server computer,is strictly prohibited. t users to make one paper copy for their THE sum +fac/(sqrtz*(zt+alamb)); fac=0.25*fac; xt=0.25*(xt+alamb) yt=0.25*(yt+alamb); zt=0.25*(zt+alamb); own ave=0.2*(xt+yt+3.0*zt) delx=(ave-xt)/ave; dely=(ave-yt)/ave; only),or send email to directcustsen delz=(ave-zt)/ave; while (FMAX(FMAX(fabs(delx),fabs(dely)),fabs(delz))>ERRTOL); ART OF SCIENTIFIC COMPUTING(ISBN 0-521 ea=delx*dely; 1988-1992 by Cambridge University Press.Programs Copyright (C)1988-1992 by Numerical Recipes eb=delz*delz; ecsea-eb; ed=ea-6.0*eb; eesedtectec; return 3.0*sum+fac*(1.0+ed*(-C1+C5*ed-C6*delz*ee) @cambridge.org personal use.Further reproduction,or +delz*(C2*ee+delz*(-C3*ec+delz*C4*ea)))/(ave*sqrt(ave)): 1-431085 Software. #include <math.h> #include "nrutil.h" (outside North America) #define ERRTOL 0.05 #define TINY 2.5e-13 ying of machine #define BIG 9.0e11 #define C1 (3.0/14.0) #def1neC2(1.0/3.0) #define C3 (3.0/22.0) #define C4 (3.0/26.0) #define C5 (0.75*C3) #define C6 (1.5米C4 #define C7 (0.5*C2) #define C8 (C3+C3)6.11 Elliptic Integrals and Jacobian Elliptic Functions 265 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). #define TINY 1.0e-25 #define BIG 4.5e21 #define C1 (3.0/14.0) #define C2 (1.0/6.0) #define C3 (9.0/22.0) #define C4 (3.0/26.0) #define C5 (0.25*C3) #define C6 (1.5*C4) float rd(float x, float y, float z) Computes Carlson’s elliptic integral of the second kind, RD(x, y, z). x and y must be non￾negative, and at most one can be zero. z must be positive. TINY must be at least twice the negative 2/3 power of the machine overflow limit. BIG must be at most 0.1 × ERRTOL times the negative 2/3 power of the machine underflow limit. { float alamb,ave,delx,dely,delz,ea,eb,ec,ed,ee,fac,sqrtx,sqrty, sqrtz,sum,xt,yt,zt; if (FMIN(x,y) < 0.0 || FMIN(x+y,z) < TINY || FMAX(FMAX(x,y),z) > BIG) nrerror("invalid arguments in rd"); xt=x; yt=y; zt=z; sum=0.0; fac=1.0; do { sqrtx=sqrt(xt); sqrty=sqrt(yt); sqrtz=sqrt(zt); alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz; sum += fac/(sqrtz*(zt+alamb)); fac=0.25*fac; xt=0.25*(xt+alamb); yt=0.25*(yt+alamb); zt=0.25*(zt+alamb); ave=0.2*(xt+yt+3.0*zt); delx=(ave-xt)/ave; dely=(ave-yt)/ave; delz=(ave-zt)/ave; } while (FMAX(FMAX(fabs(delx),fabs(dely)),fabs(delz)) > ERRTOL); ea=delx*dely; eb=delz*delz; ec=ea-eb; ed=ea-6.0*eb; ee=ed+ec+ec; return 3.0*sum+fac*(1.0+ed*(-C1+C5*ed-C6*delz*ee) +delz*(C2*ee+delz*(-C3*ec+delz*C4*ea)))/(ave*sqrt(ave)); } #include <math.h> #include "nrutil.h" #define ERRTOL 0.05 #define TINY 2.5e-13 #define BIG 9.0e11 #define C1 (3.0/14.0) #define C2 (1.0/3.0) #define C3 (3.0/22.0) #define C4 (3.0/26.0) #define C5 (0.75*C3) #define C6 (1.5*C4) #define C7 (0.5*C2) #define C8 (C3+C3)