266 Chapter 6. Special Functions float rj(float x,float y,float z,float p) Computes Carlson's elliptic integral of the third kind,RJ(,y,2,p).x.y,and 2 must be nonnegative,and at most one can be zero.p must be nonzero.If p<0,the Cauchy principal value is returned.TINY must be at least twice the cube root of the machine underflow limit, BIG at most one fifth the cube root of the machine overflow limit. float rc(float x,float y); float rf(float x,float y,float z); float a,alamb,alpha,ans,ave,b,beta,delp,delx,dely,delz,ea,eb,ec, ed,ee,fac,pt,rcx,rho,sqrtx,sqrty,sqrtz,sum,tau,xt,yt,zt; http://www.n Permission is read able files if (FMIN(FMIN(x,y),z)<0.0 II FMIN(FMIN(FMIN(x+y,x+z),y+z),fabs(p))TINY II FMAX(FMAX (FMAX(x,y),z),fabs(p))>BIG) nrerror("invalid arguments in rj"); 83 gm=0,0 fac=1.0; 1f(p>0.0) 11-800-872 (including this one) xt=xi internet yt=y; zt=z; pt=pi else xt=FMIN(FMIN(x,y),z) zt=FMAX(FMAX(x,y),z) yt=x+y+z-xt-zt; -7423 (North America to any server computer,is strictly prohibited. from NUMERICAL RECIPES IN C: THE a=1.0/(yt-p); b=a*(zt-yt)*(yt-xt); pt=yt+b; rho=xt*zt/yt; tau=p*pt/yt; rcx=rc(rho,tau); do sqrtx=sqrt(xt); sqrty=sqrt(yt); sqrtz=sgrt(zt); alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz; alpha-SQR(pt*(sqrtx+sqrty+sqrtz)+sqrtx*sqrty*sqrtz); beta=pt*SQR(pt+alamb); sum +fac*rc(alpha,beta); fac=0.25*fac; xt=0.25*(xt+alamb); only),orsend email to directcustserv@cambridge.org 1988-1992 by Cambridge University Press.Programs Copyright (C)1988-1992 by Numerical Recipes ART OF SCIENTIFIC COMPUTING(ISBN 0-521-43108-5) yt=0.25*(yt+alamb): zt=0.25*(zt+alamb) pt=0.25*(pt+alamb); personal use.Further reproduction,or ave=0.2*(xt+yt+zt+pt+pt); delx=(ave-xt)/ave; dely=(ave-yt)/ave; delz=(ave-zt)/ave; delp=(ave-pt)/ave: (outside North America) un Software. while (FMAX(FMAX(FMAX(fabs(delx),fabs(dely)), fabs(delz)),fabs(delp))>ERRTOL); ea=delx*(dely+delz)+dely*delz: eb=delx*dely*delz; ying of machine ec=delp*delp; ed=ea-3.0*ec; ee=eb+2.0*delp*(ea-ec); ans=3.0*sum+fac*(1.0+ed*(-C1+C5*ed-C6*ee)+eb*(C7+delp*(-C8+delp*C4)) +delp*ea*(C2-delp*C3)-C2*delp*ec)/(ave*sgrt(ave)); if (p <=0.0)ans=a*(b*ans+3.0*(rcx-rf(xt,yt,zt))); return ans;266 Chapter 6. Special Functions 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). float rj(float x, float y, float z, float p) Computes Carlson’s elliptic integral of the third kind, RJ (x, y, z, p). x, y, and z must be nonnegative, and at most one can be zero. p must be nonzero. If p < 0, the Cauchy principal value is returned. TINY must be at least twice the cube root of the machine underflow limit, BIG at most one fifth the cube root of the machine overflow limit. { float rc(float x, float y); float rf(float x, float y, float z); float a,alamb,alpha,ans,ave,b,beta,delp,delx,dely,delz,ea,eb,ec, ed,ee,fac,pt,rcx,rho,sqrtx,sqrty,sqrtz,sum,tau,xt,yt,zt; if (FMIN(FMIN(x,y),z) < 0.0 || FMIN(FMIN(FMIN(x+y,x+z),y+z),fabs(p)) < TINY || FMAX(FMAX(FMAX(x,y),z),fabs(p)) > BIG) nrerror("invalid arguments in rj"); sum=0.0; fac=1.0; if (p > 0.0) { xt=x; yt=y; zt=z; pt=p; } else { xt=FMIN(FMIN(x,y),z); zt=FMAX(FMAX(x,y),z); yt=x+y+z-xt-zt; a=1.0/(yt-p); b=a*(zt-yt)*(yt-xt); pt=yt+b; rho=xt*zt/yt; tau=p*pt/yt; rcx=rc(rho,tau); } do { sqrtx=sqrt(xt); sqrty=sqrt(yt); sqrtz=sqrt(zt); alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz; alpha=SQR(pt*(sqrtx+sqrty+sqrtz)+sqrtx*sqrty*sqrtz); beta=pt*SQR(pt+alamb); sum += fac*rc(alpha,beta); fac=0.25*fac; xt=0.25*(xt+alamb); yt=0.25*(yt+alamb); zt=0.25*(zt+alamb); pt=0.25*(pt+alamb); ave=0.2*(xt+yt+zt+pt+pt); delx=(ave-xt)/ave; dely=(ave-yt)/ave; delz=(ave-zt)/ave; delp=(ave-pt)/ave; } while (FMAX(FMAX(FMAX(fabs(delx),fabs(dely)), fabs(delz)),fabs(delp)) > ERRTOL); ea=delx*(dely+delz)+dely*delz; eb=delx*dely*delz; ec=delp*delp; ed=ea-3.0*ec; ee=eb+2.0*delp*(ea-ec); ans=3.0*sum+fac*(1.0+ed*(-C1+C5*ed-C6*ee)+eb*(C7+delp*(-C8+delp*C4)) +delp*ea*(C2-delp*C3)-C2*delp*ec)/(ave*sqrt(ave)); if (p <= 0.0) ans=a*(b*ans+3.0*(rcx-rf(xt,yt,zt))); return ans; }