6.9 Fresnel Integrals,Cosine and Sine Integrals 257 cc=Cadd(b,Cdiv(Complex(a,0.0),cc)); del-Cmul(cc,d); h=Cmul(h,del); if (fabs(del.r-1.0)+fabs(del.i)<EPS)break; if (k MAXIT)nrerror("cf failed in frenel"); h=Cmul(Complex(ax,-ax),h) cs=Cmul(Complex(0.5,0.5), Csub(ONE,Cmul(Complex(cos(0.5*pix2),sin(0.5*pix2)),h))) *c=cs.ri *s=Cs.1; d http://www.nr read able files 1f(x<0.0)[ Use antisymmetry. 米c=-(*c); .com or call 1-800-872- granted for (including this one) Cosine and Sine Integrals to any /Cambridge -7423 (North users to make one paper from NUMERICAL RECIPES IN C: 19881992 The cosine and sine integrals are defined by cost-1 Ci(x)=y+lnx+ dt America server computer, e University Press. THE 0 (6.9.6) ART sin t Si(z)= dt Jo t 9 Programs Here 0.5772...is Euler's constant. We only need a way to calculate the functions for x>0.because OF SCIENTIFIC Si(-z)=-Si(x), Ci(-x)=Ci(x)-iπ (6.9.7) 6 Once again we can evaluate these functions by a judicious combination of power series and complex continued fraction.The series are COMPUTING (ISBN 188810920 Si(x)=- 3.3! 5.51 23 (6.9.8) Ci(z)=y+Inz+ 22+ 4.4 Numerical Recipes 10-621 -43106 The continued fraction for the exponential integral E1(ix)is (outside B1(ix)=-Ci(x)+[Si(x)-π/2] Software. =e-i证 1112 2 ix+1+ix+1+ix+ (6.9.9) 1 12 22 1+ix-3+ix-5+ix- 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 point from the alternating series to the continued fraction is =2 in this case.As for the Fresnel integrals,for large the precision may be limited by the precision of the sine and cosine routines.6.9 Fresnel Integrals, Cosine and Sine Integrals 257 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). cc=Cadd(b,Cdiv(Complex(a,0.0),cc)); del=Cmul(cc,d); h=Cmul(h,del); if (fabs(del.r-1.0)+fabs(del.i) < EPS) break; } if (k > MAXIT) nrerror("cf failed in frenel"); h=Cmul(Complex(ax,-ax),h); cs=Cmul(Complex(0.5,0.5), Csub(ONE,Cmul(Complex(cos(0.5*pix2),sin(0.5*pix2)),h))); *c=cs.r; *s=cs.i; } if (x < 0.0) { Use antisymmetry. *c = -(*c); *s = -(*s); } } Cosine and Sine Integrals The cosine and sine integrals are defined by Ci(x) = γ + ln x +
x 0 cost − 1 t dt Si(x) =
x 0 sin t t dt (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) = Ci(x) − iπ (6.9.7) Once again we can evaluate these functions by a judicious combination of power series and complex continued fraction. The series are Si(x) = x − x3 3 · 3! + x5 5 · 5! −··· Ci(x) = γ + ln x +  − x2 2 · 2! + x4 4 · 4! −··· (6.9.8) The continued fraction for the exponential integral E1(ix) is E1(ix) = − Ci(x) + i[Si(x) − π/2] = e−ix  1 ix + 1 1 + 1 ix + 2 1 + 2 ix + ··· = e−ix  1 1 + ix − 12 3 + ix − 22 5 + ix − ··· (6.9.9) 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 point from the alternating series to the continued fraction is x = 2 in this case. As for the Fresnel integrals, for large x the precision may be limited by the precision of the sine and cosine routines