f6-11

254 Chapter6. SpecialFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).sum=sum+c(i)*(e1/d1+1./(d2*e1)) d1=d1+2. d2=d2-2. e1=e2*e1 enddo 12 dawson=0.5641895835*sign(exp(-xp**2),x)*sum Constant is 1/√π. endifreturnEND Other methods for computingDawson’s integral are also known [2,3]. CITED REFERENCES AND FURTHER READING: Rybicki, G.B. 1989, Computers in Physics , vol. 3, no. 2, pp. 85–87. [1] Cody, W.J., Pociorek, K.A., and Thatcher, H.C. 1970, Mathematics of Computation , vol. 24, pp. 171–178. [2] McCabe, J.H. 1974, Mathematics of Computation , vol. 28, pp. 811–816. [3] 6.11 Elliptic Integrals and Jacobian Elliptic Functions Elliptic integralsoccurin manyapplications,becauseanyintegralof theform /integraldisplay R(t, s)dt (6.11.1 ) where Ris a rational function of tands, and sis the square root of a cubic or quartic polynomial in t, can be evaluated in terms of elliptic integrals. Standard references [1]describe how to carry out the reduction, which was originally done by Legendre. Legendre showed that only three basic elliptic integrals are required. The simplest of these is I1=/integraldisplayx ydt/radicalbig (a1+b1t)(a2+b2t)(a3+b3t)(a4+b4t)(6.11.2 ) wherewehavewrittenthequartic s2infactoredform. Instandardintegraltables [2], one of the limits of integration is always a zero of the quartic, while the other limit lies closer than the next zero, so that there is no singularity within the interval. To evaluate I1,wesimplybreaktheinterval [y,x]intosubintervals,eachofwhicheither beginsorendsonasingularity. Thetables,therefore,needonlydistinguishtheeightcasesinwhicheachofthefourzeros(orderedaccordingtosize)appearsastheupper or lower limit of integration. In addition, when one of the b’s in (6.11.2) tends to zero, the quartic reduces to a cubic, with the largest or smallest singularity movingto±∞; this leads to eight more cases (actually just special cases of the first eight). The sixteen cases in total are then usually tabulatedin terms of Legendre’sstandard elliptic integral of the 1st kind, which we will define below. By a change of the variable of integration t, the zeros of the quartic are mapped to standard locations 6.11EllipticIntegralsandJacobianEllipticFunctions 255Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).on the real axis. Then only two dimensionless parameters are needed to tabulate Legendre’s integral. However, the symmetry of the original integral (6.11.2) under permutation of the roots is concealed in Legendre’s notation. We will get back to Legendre’s notation below. But first, here is a better way: Carlson[3]has given a new definition of a standard elliptic integral of the first kind, RF(x, y, z )=1 2/integraldisplay∞ 0dt/radicalbig (t+x)(t+y)(t+z)(6.11.3 ) where x,y, and zare nonnegative and at most one is zero. By standardizing the range of integration, he retains permutation symmetry for the zeros. (Weierstrass’ canonical formalso has this property.) Carlson first shows that when xoryis a zero of the quartic in (6.11.2), the integral I 1can be written in terms of RFin a form that is symmetric under permutation of the remaining three zeros. In the general case when neither xnoryis a zero, two such RFfunctions can be combined into a single one by an addition theorem , leading to the fundamental formula I1=2RF(U2 12,U2 13,U2 14)( 6.11.4 ) where Uij=(XiXjYkYm+YiYjXkXm)/(x−y)( 6.11.5 ) Xi=(ai+bix)1/2,Y i=(ai+biy)1/2(6.11.6 ) andi, j, k, m is any permutation of 1,2,3,4. A short-cut in evaluating these expressions is U2 13=U2 12−(a1b4−a4b1)(a2b3−a3b2) U2 14=U2 12−(a1b3−a3b1)(a2b4−a4b2)(6.11.7 ) TheU’s correspond to the three ways of pairing the four zeros, and I1is thus manifestly symmetric under permutation of the zeros. Equation (6.11.4) therefore reproduces all sixteencases when one limit is a zero, and also includes the cases when neither limit is a zero. Thus Carlson’s function allows arbitrary ranges of integration and arbitrary positions of the branch points of the integrand relative to the interval of integration. To handle ellipticintegralsofthesecondandthirdkind,Carlsondefinesthestandardintegralofthethirdkindas R J(x, y, z, p )=3 2/integraldisplay∞ 0dt (t+p)/radicalbig (t+x)(t+y)(t+z)(6.11.8 ) which is symmetric in x,y, and z. The degenerate case when two arguments are equal is denoted RD(x, y, z )=RJ(x, y, z, z )( 6.11.9 ) and is symmetric in xandy. The function RDreplaces Legendre’s integral of the second kind. The degenerate form of RFis denoted RC(x, y )=RF(x, y, y )( 6.11.10 ) It embraces logarithmic, inverse circular, and inverse hyperbolic functions. Carlson[4-7]gives integral tables in terms of the exponents of the linear factors of the quartic in (6.11.1). For example, the integral where the exponents are (1 2,1 2,−1 2,−3 2) can be expressed as a single integral in terms of RD; it accounts for 144 separate cases in Gradshteyn and Ryzhik [2]! Refer to Carlson’s papers [3-7]for some of the practical details in reducing elliptic integrals to his standard forms, such as handling complex conjugate zeros. 256 Chapter6. SpecialFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Turn now to the numerical evaluation of elliptic integrals. The traditional methods [8] are Gauss or Landen transformations. Descending transformations decrease the modulus kof the Legendre integrals towards zero, increasing transformations increase it towards unity. In these limits the functions have simple analytic expressions. While these methodsconverge quadratically and are quite satisfactory for integrals of the first and second kinds,they generally lead to loss of significant figures in certain regimes for integrals of the thirdkind. Carlson’s algorithms [9,10], by contrast, provide a unified method for all three kinds with no significant cancellations. The key ingredient in these algorithms is the duplication theorem : RF(x, y, z )=2 RF(x+λ, y +λ, z +λ) =RF/parenleftbiggx+λ 4,y+λ 4,z+λ 4/parenrightbigg(6.11.11 ) where λ=(xy)1/2+(xz)1/2+(yz)1/2(6.11.12 ) This theorem can be proved by a simple change of variable of integration [11]. Equation (6.11.11) isiterateduntilthearguments of RFarenearly equal. Forequal argumentswehave RF(x, x, x )=x−1/2(6.11.13 ) Whentheargumentsarecloseenough, thefunctionisevaluated fromafixedTaylorexpansion about (6.11.13) through fifth-order terms. While the iterative part of the algorithm is onlylinearlyconvergent,theerrorultimatelydecreases byafactorof 4 6= 4096foreachiteration. Typically only two or three iterations are required, perhaps six or seven if the initial valuesof the arguments have huge ratios. We list the algorithm for R Fhere, and refer you to Carlson’s paper [9]for the other cases. Stage 1: For n=0,1,2,...compute µn=(xn+yn+zn)/3 Xn=1−(xn/µ n),Y n=1−(yn/µ n),Z n=1−(zn/µ n) /epsilon1n=m a x ( |Xn|,|Yn|,|Zn|) If/epsilon1n<tolgo to Stage 2; else compute λn=(xnyn)1/2+(xnzn)1/2+(ynzn)1/2 xn+1=(xn+λn)/4,y n+1=(yn+λn)/4,z n+1=(zn+λn)/4 and repeat this stage. Stage 2: Compute E2=XnYn−Z2 n,E 3=XnYnZn RF=( 1−1 10E2+1 14E3+1 24E2 2−3 44E2E3)/(µn)1/2 In some applications the argument pinRJor the argument yinRCisnegative, and the Cauchyprincipalvalueoftheintegralisrequired. Thisiseasilyhandledbyusingtheformulas RJ(x, y,z, p )= [(γ−y)RJ(x, y, z, γ )−3RF(x, y, z )+3 RC(xz/y, pγ/y )]/(y−p) (6.11.14 ) where γ≡y+(z−y)(y−x) y−p(6.11.15 ) 6.11EllipticIntegralsandJacobianEllipticFunctions 257Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).is positive if pis negative, and RC(x, y )=/parenleftbiggx x−y/parenrightbigg1/2 RC(x−y,−y)( 6.11.16 ) The Cauchy principal value of RJhas a zero at some value of p< 0, so (6.11.14) will give some loss of significant figures near the zero. FUNCTION rf(x,y,z) REAL rf,x,y,z,ERRTOL,TINY,BIG,THIRD,C1,C2,C3,C4 PARAMETER (ERRTOL=.08,TINY=1.5e-38,BIG=3.E37,THIRD=1./3., * C1=1./24.,C2=.1,C3=3./44.,C4=1./14.) Computes Carlson’s elliptic integral of the first kind, RF(x, y, z ).x,y,a n d zmust be nonnegative, and at most one can be zero. TINY must be at least 5 times the machine underflow limit, BIG at most one fifth the machine overflow limit. REAL alamb,ave,delx,dely,delz,e2,e3,sqrtx,sqrty,sqrtz,xt,yt,ztif(min(x,y,z).lt.0..or.min(x+y,x+z,y+z).lt.TINY.or. * max(x,y,z).gt.BIG)pause ’invalid arguments in rf’ xt=xyt=yzt=z 1 continue sqrtx=sqrt(xt)sqrty=sqrt(yt) sqrtz=sqrt(zt) alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtzxt=.25*(xt+alamb)yt=.25*(yt+alamb) zt=.25*(zt+alamb) ave=THIRD*(xt+yt+zt)delx=(ave-xt)/ave dely=(ave-yt)/ave delz=(ave-zt)/ave if(max(abs(delx),abs(dely),abs(delz)).gt.ERRTOL)goto 1e2=delx*dely-delz**2 e3=delx*dely*delz rf=(1.+(C1*e2-C2-C3*e3)*e2+C4*e3)/sqrt(ave)return END A value of 0.08 for the error tolerance parameter is adequate for single precision (7 significant digits). Since the error scales as /epsilon16 n,we see that 0.0025 willyield double precision (16 significant digits) and require at most two or three more iterations. Since the coefficientsofthesixth-ordertruncation erroraredifferentfortheother ellipticfunctions, thesevaluesfortheerrortoleranceshouldbechangedto0.04and0.0012inthealgorithmfor R C,and0.05and 0.0015for RJandRD. Aswellasbeinganalgorithminitsownrightforcertaincombinations of elementary functions, the algorithm for RCis used repeatedly in the computation of RJ. The Fortran implementations test the input arguments against two machine-dependent constants, TINYandBIG, to ensure that there will be no underflow or overflow during the computation. We have chosen conservative values, corresponding to a machine minimum of3×10 −39and a machine maximum of 1.7×1038. You can always extend the range of admissible argument values by using the homogeneity relations (6.11.22), below. FUNCTION rd(x,y,z) REAL rd,x,y,z,ERRTOL,TINY,BIG,C1,C2,C3,C4,C5,C6 PARAMETER (ERRTOL=.05,TINY=1.e-25,BIG=4.5E21,C1=3./14.,C2=1./6., * C3=9./22.,C4=3./26.,C5=.25*C3,C6=1.5*C4) Computes Carlson’s elliptic integral of the second kind, RD(x, y, z ).xandymust be nonnegative, and at most one can be zero. zmust 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. REAL alamb,ave,delx,dely,delz,ea,eb,ec,ed,ee,fac,sqrtx,sqrty, 258 Chapter6. SpecialFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).* sqrtz,sum,xt,yt,zt if(min(x,y).lt.0..or.min(x+y,z).lt.TINY.or. * max(x,y,z).gt.BIG)pause ’invalid arguments in rd’ xt=xyt=y zt=z sum=0.fac=1. 1 continue sqrtx=sqrt(xt) sqrty=sqrt(yt)sqrtz=sqrt(zt)alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtz sum=sum+fac/(sqrtz*(zt+alamb)) fac=.25*facxt=.25*(xt+alamb) yt=.25*(yt+alamb) zt=.25*(zt+alamb)ave=.2*(xt+yt+3.*zt)delx=(ave-xt)/ave dely=(ave-yt)/ave delz=(ave-zt)/ave if(max(abs(delx),abs(dely),abs(delz)).gt.ERRTOL)goto 1 ea=delx*dely eb=delz*delzec=ea-ebed=ea-6.*eb ee=ed+ec+ec rd=3.*sum+fac*(1.+ed*(-C1+C5*ed-C6*delz*ee) * +delz*(C2*ee+delz*(-C3*ec+delz*C4*ea)))/(ave*sqrt(ave)) return END FUNCTION rj(x,y,z,p) REAL rj,p,x,y,z,ERRTOL,TINY,BIG,C1,C2,C3,C4,C5,C6,C7,C8PARAMETER (ERRTOL=.05,TINY=2.5e-13,BIG=9.E11,C1=3./14.,C2=1./3., * C3=3./22.,C4=3./26.,C5=.75*C3,C6=1.5*C4,C7=.5*C2,C8=C3+C3) C USES rc,rf Computes Carlson’s elliptic integral of the third kind, RJ(x, y, z, p ).x,y,a n d zmust be nonnegative, and at most one can be zero. pmust 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. REAL a,alamb,alpha,ave,b,beta,delp,delx,dely,delz,ea,eb,ec, * ed,ee,fac,pt,rcx,rho,sqrtx,sqrty,sqrtz,sum,tau,xt, * yt,zt,rc,rf if(min(x,y,z).lt.0..or.min(x+y,x+z,y+z,abs(p)).lt.TINY.or. * max(x,y,z,abs(p)).gt.BIG)pause ’invalid arguments in rj’ sum=0. fac=1.if(p.gt.0.)then xt=x yt=y zt=zpt=p else xt=min(x,y,z)zt=max(x,y,z)yt=x+y+z-xt-zt a=1./(yt-p) b=a*(zt-yt)*(yt-xt)pt=yt+b rho=xt*zt/yt 6.11EllipticIntegralsandJacobianEllipticFunctions 259Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).tau=p*pt/yt rcx=rc(rho,tau) endif 1 continue sqrtx=sqrt(xt) sqrty=sqrt(yt) sqrtz=sqrt(zt)alamb=sqrtx*(sqrty+sqrtz)+sqrty*sqrtzalpha=(pt*(sqrtx+sqrty+sqrtz)+sqrtx*sqrty*sqrtz)**2 beta=pt*(pt+alamb)**2 sum=sum+fac*rc(alpha,beta)fac=.25*facxt=.25*(xt+alamb) yt=.25*(yt+alamb) zt=.25*(zt+alamb)pt=.25*(pt+alamb) ave=.2*(xt+yt+zt+pt+pt) delx=(ave-xt)/avedely=(ave-yt)/avedelz=(ave-zt)/ave delp=(ave-pt)/ave if(max(abs(delx),abs(dely),abs(delz),abs(delp)).gt.ERRTOL)goto 1ea=delx*(dely+delz)+dely*delz eb=delx*dely*delz ec=delp**2ed=ea-3.*ecee=eb+2.*delp*(ea-ec) rj=3.*sum+fac*(1.+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.le.0.) rj=a*(b*rj+3.*(rcx-rf(xt,yt,zt))) return END FUNCTION rc(x,y) REAL rc,x,y,ERRTOL,TINY,SQRTNY,BIG,TNBG,COMP1,COMP2,THIRD, * C1,C2,C3,C4 PARAMETER (ERRTOL=.04,TINY=1.69e-38,SQRTNY=1.3e-19,BIG=3.E37, * TNBG=TINY*BIG,COMP1=2.236/SQRTNY,COMP2=TNBG*TNBG/25.,* THIRD=1./3.,C1=.3,C2=1./7.,C3=.375,C4=9./22.) Computes Carlson’s degenerate elliptic integral, R C(x, y).xmust 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. REAL alamb,ave,s,w,xt,yt if(x.lt.0..or.y.eq.0..or.(x+abs(y)).lt.TINY.or.(x+abs(y)).gt.BIG * .or.(y.lt.-COMP1.and.x.gt.0..and.x.lt.COMP2))* pause ’invalid arguments in rc’ if(y.gt.0.)then xt=xyt=yw=1. else xt=x-yyt=-y w=sqrt(x)/sqrt(xt) endif 1 continue alamb=2.*sqrt(xt)*sqrt(yt)+yt xt=.25*(xt+alamb) yt=.25*(yt+alamb)ave=THIRD*(xt+yt+yt) s=(yt-ave)/ave 260 Chapter6. SpecialFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).if(abs(s).gt.ERRTOL)goto 1 rc=w*(1.+s*s*(C1+s*(C2+s*(C3+s*C4))))/sqrt(ave) return END 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. TheLegendre elliptic integral of the 1st kind is defined as F(φ, k)≡/integraldisplayφ 0dθ/radicalbig 1−k2sin2θ(6.11.17 ) Thecomplete elliptic integral of the 1st kind is given by K(k)≡F(π/2,k)( 6.11.18 ) In terms of RF, F(φ, k)=s i n φRF(cos2φ,1−k2sin2φ,1) K(k)=RF(0,1−k2,1)(6.11.19 ) TheLegendre elliptic integral of the 2nd kind and thecomplete elliptic integral of the 2nd kind are given by E(φ, k)≡/integraldisplayφ 0/radicalbig 1−k2sin2θd θ =s i n φRF(cos2φ,1−k2sin2φ,1) −1 3k2sin3φRD(cos2φ,1−k2sin2φ,1) E(k)≡E(π/2,k)=RF(0,1−k2,1)−1 3k2RD(0,1−k2,1)(6.11.20 ) Finally, the Legendre elliptic integral of the 3rd kind is Π(φ, n, k )≡/integraldisplayφ 0dθ (1 + nsin2θ)/radicalbig 1−k2sin2θ =s i n φRF(cos2φ,1−k2sin2φ,1) −1 3nsin3φRJ(cos2φ,1−k2sin2φ,1,1+nsin2φ)(6.11.21 ) (Notethatthissignconventionfor nisoppositethatofAbramowitzandStegun [12], and that their sinαis our k.) FUNCTION ellf(phi,ak) REAL ellf,ak,phi C USES rf Legendre elliptic integral of the 1st kind F(φ, k), evaluated using Carlson’s function RF. The argument ranges are 0≤φ≤π/2,0≤ksinφ≤1. REAL s,rf s=sin(phi) ellf=s*rf(cos(phi)**2,(1.-s*ak)*(1.+s*ak),1.)returnEND 6.11EllipticIntegralsandJacobianEllipticFunctions 261Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).FUNCTION elle(phi,ak) REAL elle,ak,phi C USES rd,rf Legendre elliptic integral of the 2nd kind E(φ, k), evaluated using Carlson’s functions RD andRF. The argument ranges are 0≤φ≤π/2,0≤ksinφ≤1. REAL cc,q,s,rd,rf s=sin(phi)cc=cos(phi)**2q=(1.-s*ak)*(1.+s*ak) elle=s*(rf(cc,q,1.)-((s*ak)**2)*rd(cc,q,1.)/3.) returnEND FUNCTION ellpi(phi,en,ak) REAL ellpi,ak,en,phi C USES rf,rj Legendre elliptic integral of the 3rd kind Π(φ, n, k ), evaluated using Carlson’s functions RJ andRF. (Note that the sign convention on nis opposite that of Abramowitz and Stegun.) The ranges of φandkare0≤φ≤π/2,0≤ksinφ≤1. REAL cc,enss,q,s,rf,rjs=sin(phi) enss=en*s*s cc=cos(phi)**2q=(1.-s*ak)*(1.+s*ak) ellpi=s*(rf(cc,q,1.)-enss*rj(cc,q,1.,1.+enss)/3.) returnEND Carlson’s functions are homogeneous of degree −1 2and−3 2,s o RF(λx, λy, λz )=λ−1/2RF(x, y, z ) RJ(λx, λy, λz, λp )=λ−3/2RJ(x, y, z, p )(6.11.22 ) Thus to express a Carlson function in Legendre’s notation, permute the first three arguments into ascending order, use homogeneity to scale the third argument to be 1, and then use equations (6.11.19)–(6.11.21). Jacobian EllipticFunctions The Jacobian elliptic function sn is defined as follows: instead of considering the elliptic integral u(y,k)≡u=F(φ, k)( 6.11.23 ) consider the inversefunction y=s i n φ=sn(u, k)( 6.11.24 ) Equivalently, u=/integraldisplaysn 0dy/radicalbig (1−y2)(1−k2y2)(6.11.25 ) When k=0, sn is just sin. The functionscn and dn are defined bythe relations sn2+cn2=1,k2sn2+dn2=1 ( 6.11.26 ) The routine given below actually takes mc≡k2 c=1−k2as an input parameter. It also computes all three functions sn, cn, and dn since computing all three is no harderthancomputingany oneof them. For a descriptionof the method,see [8]. 262 Chapter6. SpecialFunctionsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).SUBROUTINE sncndn(uu,emmc,sn,cn,dn) REAL cn,dn,emmc,sn,uu,CA PARAMETER (CA=.0003) The accuracy is the square of CA. Returns the Jacobian elliptic functions sn (u, k c),c n(u, k c), and dn (u, k c).H e r e uu=u, while emmc =k2 c. INTEGER i,ii,l REAL a,b,c,d,emc,u,em(13),en(13)LOGICAL boemc=emmc u=uu if(emc.ne.0.)then bo=(emc.lt.0.)if(bo)then d=1.-emc emc=-emc/dd=sqrt(d) u=d*u endifa=1.dn=1. do 11i=1,13 l=iem(i)=a emc=sqrt(emc) en(i)=emcc=0.5*(a+emc)if(abs(a-emc).le.CA*a)goto 1 emc=a*emc a=c enddo 11 1 u=c*u sn=sin(u)cn=cos(u)if(sn.eq.0.)goto 2 a=cn/sn c=a*cdo 12ii=l,1,-1 b=em(ii) a=c*a c=dn*cdn=(en(ii)+a)/(b+a) a=c/b enddo 12 a=1./sqrt(c**2+1.)if(sn.lt.0.)then sn=-a else sn=a endif cn=c*sn 2 if(bo)then a=dn dn=cn cn=asn=sn/d endif else cn=1./cosh(u)dn=cn sn=tanh(u) endifreturnEND 6.12HypergeometricFunctions 263Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).CITED REFERENCES AND FURTHER READING: Erd´elyi, A., Magnus, W., Oberhettinger, F., and Tricomi, F.G. 1953, Higher Transcendental Functions , Vol. II, (New York: McGraw-Hill). [1] Gradshteyn, I.S., and Ryzhik, I.W. 1980, Table of Integrals, Series, and Products (New York: Academic Press). [2] Carlson, B.C. 1977, SIAM Journal on Mathematical Analysis , vol. 8, pp. 231–242. [3] Carlson,B.C.1987, MathematicsofComputation ,vol.49,pp.595–606[4];1988, op.cit.,vol.51, pp.267–280[5];1989, op.cit.,vol.53,pp.327–333[6];1991, op.cit.,vol.56,pp.267–280. [7] Bulirsch,R.1965, NumerischeMathematik ,vol.7,pp.78–90;1965, op.cit.,vol.7,pp.353–354; 1969,op. cit., vol. 13, pp. 305–315. [8] Carlson, B.C. 1979, Numerische Mathematik , vol. 33, pp. 1–16. [9] Carlson, B.C., and Notis, E.M. 1981, ACM Transactions on Mathematical Software , vol. 7, pp. 398–403. [10] Carlson, B.C. 1978, SIAM Journal on Mathematical Analysis , vol. 9, p. 524–528. [11] Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions , Applied Mathe- matics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 byDover Publications, New York), Chapter 17. [12] Mathews, J., and Walker, R.L. 1970, Mathematical Methods of Physics , 2nd ed. (Reading, MA: W.A. Benjamin/Addison-Wesley), pp. 78–79. 6.12 Hypergeometric Functions As was discussed in §5.14,a fast, generalroutineforthe the complexhyperge- ometricfunction 2F1(a, b, c ;z),is difficultorimpossible. Thefunctionis definedas the analytic continuation of the hypergeometric series, 2F1(a, b, c ;z)=1+ab cz 1!+a(a+1 )b(b+1 ) c(c+1 )z2 2!+··· +a(a+1 ) ...(a+j−1)b(b+1 ) ...(b+j−1) c(c+1 ) ...(c+j−1)zj j!+··· (6.12.1 ) This series converges only within the unit circle |z|<1(see[1]), but one’s interest in the function is not confined to this region. Section 5.14 discussed the method of evaluating this function by direct path integrationin the complex plane. We here merely list the routines that result. Implementation of the function hypgeois straightforward, and is described by comments in the program. The machinery associated with Chapter 16’s routine forintegratingdifferentialequations, odeint, is onlyminimallyintrusive,and need not even be completely understood: use of odeintrequires a common block with one zeroed variable, one subroutine call, and a prescribed format for the derivative routine hypdrv. The subroutine hypgeowill fail, of course, for values of ztoo close to the singularity at 1. (If you need to approach this singularity, or the one at ∞, use the “lineartransformationformulas”in §15.3of[1].) Awayfrom z=1,andformoderate values of a, b, c, it is often remarkable how few steps are required to integrate the equations. A half-dozen is typical.