f6-9

248 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).pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m) pmm=pmmp1 pmmp1=pll enddo 12 plgndr=pll endif endifreturnEND CITED REFERENCES AND FURTHER READING: Magnus, W., and Oberhettinger, F. 1949, Formulas and Theorems for the Functions of Mathe- matical Physics (New York: Chelsea), pp. 54ff. [1] 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 8. [2] 6.9 FresnelIntegrals,CosineandSineIntegrals Fresnel Integrals The two Fresnel integrals are defined by C(x)=/integraldisplayx 0cos/parenleftBigπ 2t2/parenrightBig dt, S (x)=/integraldisplayx 0sin/parenleftBigπ 2t2/parenrightBig dt (6.9.1 ) The mostconvenientway of evaluatingthese functionsto arbitraryprecisionis to use powerseries forsmall xanda continuedfractionforlarge x. The series are C(x)=x−/parenleftBigπ 2/parenrightBig2x5 5·2!+/parenleftBigπ 2/parenrightBig4x9 9·4!−··· S(x)=/parenleftBigπ 2/parenrightBigx3 3·1!−/parenleftBigπ 2/parenrightBig3x7 7·3!+/parenleftBigπ 2/parenrightBig5x11 11·5!−···(6.9.2 ) There is a complex continued fraction that yields both S(x)andC(x)simul- taneously: C(x)+iS(x)=1+i 2erfz, z =√π 2(1−i)x (6.9.3 ) where ez2erfcz=1√π/parenleftbigg1 z+1/2 z+1 z+3/2 z+2 z+···/parenrightbigg =2z√π/parenleftbigg1 2z2+1−1·2 2z2+5−3·4 2z2+9−···/parenrightbigg (6.9.4 ) 6.9FresnelIntegrals,CosineandSineIntegrals 249Sample 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).In the last line we have converted the “standard” form of the continued fraction to its “even” form (see §5.2), which converges twice as fast. We must be careful not to evaluate the alternating series (6.9.2) at too large a value of x; inspection of the terms shows that x=1.5is a goodpointto switch overto the continuedfraction. Note that for large x C(x)∼1 2+1 πxsin/parenleftBigπ 2x2/parenrightBig ,S (x)∼1 2−1 πxcos/parenleftBigπ 2x2/parenrightBig (6.9.5 ) Thus the precision of the routine frenelmay be limited by the precision of the library routines for sine and cosine for large x. SUBROUTINE frenel(x,s,c) INTEGER MAXIT REAL c,s,x,EPS,FPMIN,PI,PIBY2,XMIN PARAMETER (EPS=6.e-8,MAXIT=100,FPMIN=1.e-30,XMIN=1.5, * PI=3.1415927,PIBY2=1.5707963) Computes the Fresnel integrals S(x)andC(x)for all real x. Parameters: EPS is the relative error; MAXIT is the maximum number of iterations allowed; FPMIN is a number near the smallest representable floating-point number; XMIN is the dividing line between using the series and continued fraction; PI =π;PIBY2 =π/2. INTEGER k,n REAL a,absc,ax,fact,pix2,sign,sum,sumc,sums,term,testCOMPLEX b,cc,d,h,del,csLOGICAL odd absc(h)=abs(real(h))+abs(aimag(h)) Statement function. ax=abs(x)if(ax.lt.sqrt(FPMIN))then Special case: avoid failure of convergence test because of underflow. s=0. c=ax else if(ax.le.XMIN)then Evaluate both series simultaneously. sum=0. sums=0. sumc=axsign=1. fact=PIBY2*ax*ax odd=.true.term=axn=3 do 11k=1,MAXIT term=term*fact/ksum=sum+sign*term/ntest=abs(sum)*EPS if(odd)then sign=-signsums=sum sum=sumc else sumc=sumsum=sums endif if(term.lt.test)goto 1odd=.not.odd n=n+2 enddo 11 pause ’series failed in frenel’ 1 s=sums c=sumc else Evaluate continued fraction by modified Lentz’s method ( §5.2). pix2=PI*ax*ax b=cmplx(1.,-pix2) 250 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).cc=1./FPMIN d=1./b h=d n=-1do 12k=2,MAXIT n=n+2 a=-n*(n+1)b=b+4.d=1./(a*d+b) Denominators cannot be zero. cc=b+a/cc del=cc*dh=h*delif(absc(del-1.).lt.EPS)goto 2 enddo 12 pause ’cf failed in frenel’ 2 h=h*cmplx(ax,-ax) cs=cmplx(.5,.5)*(1.-cmplx(cos(.5*pix2),sin(.5*pix2))*h) c=real(cs)s=aimag(cs) endif if(x.lt.0.)then Use antisymmetry. c=-cs=-s endif returnEND Cosine and Sine Integrals The cosine and sine integrals are defined by Ci(x)=γ+l n x+/integraldisplayx 0cost−1 tdt Si(x)=/integraldisplayx 0sint tdt(6.9.6 ) Here γ≈0.5772 ...is Euler’s constant. We only need a way to calculate the functions for x> 0, because Si(−x)=−Si(x), Ci(−x)=C i ( x)−iπ (6.9.7 ) Onceagainwecanevaluatethesefunctionsbyajudiciouscombinationofpower series and complex continued fraction. The series are Si(x)=x−x3 3·3!+x5 5·5!−··· Ci(x)=γ+l n x+/parenleftbigg −x2 2·2!+x4 4·4!−···/parenrightbigg (6.9.8 ) The continued fraction for the exponential integral E1(ix)is E1(ix)=−Ci(x)+i[Si(x)−π/2] =e−ix/parenleftbigg1 ix+1 1+1 ix+2 1+2 ix+···/parenrightbigg =e−ix/parenleftbigg1 1+ix−12 3+ix−22 5+ix−···/parenrightbigg(6.9.9 ) 6.9FresnelIntegrals,CosineandSineIntegrals 251Sample 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).The “even” form of the continued fraction is given in the last line and converges twice as fast for about the same amount of computation. A good crossover pointfrom the alternating series to the continued fraction is x=2in this case. As for the Fresnel integrals, for large xthe precision may be limited by the precision of the sine and cosine routines. SUBROUTINE cisi(x,ci,si) INTEGER MAXIT REAL ci,si,x,EPS,EULER,PIBY2,FPMIN,TMIN PARAMETER (EPS=6.e-8,EULER=.57721566,MAXIT=100,PIBY2=1.5707963, * FPMIN=1.e-30,TMIN=2.) Computes the cosine and sine integrals Ci(x)and Si(x).Ci(0) is returned as a large negative number and no error message is generated. For x< 0the routine returns Ci(−x)and you must supply the −iπyourself. Parameters: EPS is the relative error, or absolute error near a zero of Ci(x);EULER =γ; MAXIT is the maximum number of iterations allowed; PIBY2 =π/2;FPMIN is a number near the smallest representable floating-point number; TMIN is the dividing line between using the series and continued fraction. INTEGER i,k REAL a,err,fact,sign,sum,sumc,sums,t,term,absc COMPLEX h,b,c,d,delLOGICAL odd absc(h)=abs(real(h))+abs(aimag(h)) Statement function. t=abs(x)if(t.eq.0.)then Special case. si=0. ci=-1./FPMIN return endif if(t.gt.TMIN)then Evaluate continued fraction by modified Lentz’s method ( §5.2). b=cmplx(1.,t) c=1./FPMINd=1./b h=d do 11i=2,MAXIT a=-(i-1)**2 b=b+2. d=1./(a*d+b) Denominators cannot be zero. c=b+a/cdel=c*d h=h*del if(absc(del-1.).lt.EPS)goto 1 enddo 11 pause ’cf failed in cisi’ 1 continue h=cmplx(cos(t),-sin(t))*hci=-real(h) si=PIBY2+aimag(h) else Evaluate both series simultaneously. if(t.lt.sqrt(FPMIN))then Special case: avoid failure of convergence test because of underflow. sumc=0. sums=t else sum=0. sums=0. sumc=0.sign=1.fact=1. odd=.true. do 12k=1,MAXIT fact=fact*t/k term=fact/k 252 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+sign*term err=term/abs(sum) if(odd)then sign=-signsums=sum sum=sumc else sumc=sumsum=sums endif if(err.lt.EPS)goto 2odd=.not.odd enddo 12 pause ’maxits exceeded in cisi’ endif 2 si=sums ci=sumc+log(t)+EULER endifif(x.lt.0.)si=-sireturn END CITED REFERENCES AND FURTHER READING: Stegun, I.A., and Zucker, R. 1976, Journal of Research of the National Bureau of Standards , vol. 80B, pp. 291–311; 1981, op. cit., vol. 86, pp. 661–686. 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), Chapters 5 and 7. 6.10 Dawson’s Integral Dawson’s Integral F(x)is defined by F(x)=e−x2/integraldisplayx 0et2dt (6.10.1 ) The function can also be related to the complex error function by F(z)=i√π 2e−z2[1−erfc (−iz)]. (6.10.2 ) A remarkable approximation for F(z), due to Rybicki [1],i s F(z) = lim h→01√π/summationdisplay nodde−(z−nh)2 n(6.10.3 ) What makes equation (6.10.3) unusual is that its accuracy increases exponentially ashgets small, so that quite moderatevalues of h(and correspondinglyquite rapid convergence of the series) give very accurate approximations.