252 Chapter 6.Special Functions CITED REFERENCES AND FURTHER READING: Barnett,A.R..Feng.D.H.,Steed,J.W.,and Goldfarb,L.J.B.1974,Computer Physics Commu- nications,vol.8,pp.377-395.[1] Temme,N.M.1976,Journal of Computational Physics,vol.21,pp.343-350 [2]:1975,op.cit., vol.19,pp.324-337.[3] Thompson,I.J.,and Barnett,A.R.1987,Computer Physics Communications,vol.47,pp.245- 257.[4 Barnett,A.R.1981,Computer Physics Communications,vol.21,pp.297-314. Thompson,I.J.,and Barnett,A.R.1986,Journal of Computationa/Physics,vol.64,pp.490-509. Abramowitz,M.,and Stegun,I.A.1964,Handbook of Mathematical Functions,Applied Mathe- matics Series,Volume 55(Washington:National Bureau of Standards;reprinted 1968 by Dover Publications,New York),Chapter 10. NUMERICAL 6.8 Spherical Harmonics 豆三。 令 Spherical harmonics occur in a large variety of physical problems,for ex- ample,whenever a wave equation,or Laplace's equation,is solved by separa- Press. tion of variables in spherical coordinates.The spherical harmonic Yim(,), -10.With x cos0,these are defined in terms of the ordinary Legendre polynomials (cf.$4.5 and $5.5)by P(a)=(-1)m(1-x2me domP(r) (6.8.4)
252 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 machinereadable 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). CITED REFERENCES AND FURTHER READING: Barnett, A.R., Feng, D.H., Steed, J.W., and Goldfarb, L.J.B. 1974, Computer Physics Communications, vol. 8, pp. 377–395. [1] Temme, N.M. 1976, Journal of Computational Physics, vol. 21, pp. 343–350 [2]; 1975, op. cit., vol. 19, pp. 324–337. [3] Thompson, I.J., and Barnett, A.R. 1987, Computer Physics Communications, vol. 47, pp. 245– 257. [4] Barnett, A.R. 1981, Computer Physics Communications, vol. 21, pp. 297–314. Thompson, I.J., and Barnett, A.R. 1986, Journal of Computational Physics, vol. 64, pp. 490–509. Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 by Dover Publications, New York), Chapter 10. 6.8 Spherical Harmonics Spherical harmonics occur in a large variety of physical problems, for example, whenever a wave equation, or Laplace’s equation, is solved by separation of variables in spherical coordinates. The spherical harmonic Ylm(θ, φ), −l ≤ m ≤ l, is a function of the two coordinates θ, φ on the surface of a sphere. The spherical harmonics are orthogonal for different l and m, and they are normalized so that their integrated square over the sphere is unity: 2π 0 dφ 1 −1 d(cos θ)Ylm*(θ, φ)Ylm(θ, φ) = δllδmm (6.8.1) Here asterisk denotes complex conjugation. Mathematically, the spherical harmonics are related to associated Legendre polynomials by the equation Ylm(θ, φ) = 2l + 1 4π (l − m)! (l + m)!P m l (cos θ)eimφ (6.8.2) By using the relation Yl,−m(θ, φ)=(−1)mYlm*(θ, φ) (6.8.3) we can always relate a spherical harmonic to an associated Legendre polynomial with m ≥ 0. With x ≡ cos θ, these are defined in terms of the ordinary Legendre polynomials (cf. §4.5 and §5.5) by P m l (x)=(−1)m(1 − x2) m/2 dm dxm Pl(x) (6.8.4)
6.8 Spherical Harmonics 253 The first few associated Legendre polynomials,and their corresponding nor- malized spherical harmonics,are P8(z)=1 Y00= 标 P(x)=-(1-x2)1/2 h1=-V房sin9eo PO(r)=z Y10= V层cos0 P(x)=3(1-x2) y2=号V2尝sin20e2io P(x)=-3(1-x2)1/2x Y21=-V8m 15 sin 0 cos ei 83g P(a)=是(3x2-1) Y20=√（侵cos20-) (6.8.5) 11800 There are many bad ways to evaluate associated Legendre polynomials numer- ically.For example,there are explicit expressions,such as m因=0+a-2 (0-m)(m+l+1) 2mm!(I-m)! 1(m+1) +-=m0-m+m+1+2(（） 2(m+1)(m+2) 5今 Press. (6.8.6) where the polynomial continues up through the term in (1-x)-m.(See [1]for this and related formulas.This is not a satisfactory method because evaluation 、g罗色% OF SCIENTIFIC( of the polynomial involves delicate cancellations between successive terms,which alternate in sign.For large l,the individual terms in the polynomial become very much larger than their sum,and all accuracy is lost. In practice,(6.8.6)can be used only in single precision (32-bit)for l up to 6 or 8,and in double precision(64-bit)for I up to 15 or 18,depending on the precision required for the answer.A more robust computational procedure is therefore desirable,as follows: The associated Legendre functions satisfy numerous recurrence relations,tab- 10-521 ulated in[1-21.These are recurrences on I alone,on m alone,and on both l Numerica and m simultaneously.Most of the recurrences involving m are unstable,and so 43106 dangerous for numerical work.The following recurrence on I is,however,stable (compare 5.5.1): (l-m)Pm=x(2l-1)P1-(l+m-1)P2 (6.8.7) It is useful because there is a closed-form expression for the starting value, Pm=(-1)m(2m-1)(1-x2)m/2 (6.8.8) (The notation n!!denotes the product of all odd integers less than or equal to n.) Using (6.8.7)with l =m +1,and setting Prm1=0,we find Pm+1=(2m+1)Pm (6.8.9) Equations(6.8.8)and(6.8.9)provide the two starting values required for(6.8.7) for general l. The function that implements this is
6.8 Spherical Harmonics 253 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 machinereadable 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). The first few associated Legendre polynomials, and their corresponding normalized spherical harmonics, are P0 0 (x)= 1 Y00 = 1 4π P1 1 (x) = − (1 − x2)1/2 Y11 = − 3 8π sin θeiφ P0 1 (x) = x Y10 = 3 4π cos θ P2 2 (x) = 3 (1 − x2) Y22 = 1 4 15 2π sin2 θe2iφ P1 2 (x) = −3 (1 − x2)1/2x Y21 = − 15 8π sin θ cos θeiφ P0 2 (x) = 1 2 (3x2 − 1) Y20 = 5 4π ( 3 2 cos2 θ − 1 2 ) (6.8.5) There are many bad ways to evaluate associated Legendre polynomials numerically. For example, there are explicit expressions, such as P m l (x) = (−1)m(l + m)! 2mm!(l − m)! (1 − x2) m/2 1 − (l − m)(m + l + 1) 1!(m + 1) 1 − x 2 + (l − m)(l − m − 1)(m + l + 1)(m + l + 2) 2!(m + 1)(m + 2) 1 − x 2 2 −··· (6.8.6) where the polynomial continues up through the term in (1 − x)l−m. (See [1] for this and related formulas.) This is not a satisfactory method because evaluation of the polynomial involves delicate cancellations between successive terms, which alternate in sign. For large l, the individual terms in the polynomial become very much larger than their sum, and all accuracy is lost. In practice, (6.8.6) can be used only in single precision (32-bit) for l up to 6 or 8, and in double precision (64-bit) for l up to 15 or 18, depending on the precision required for the answer. A more robust computational procedure is therefore desirable, as follows: The associated Legendre functions satisfy numerous recurrence relations, tabulated in [1-2]. These are recurrences on l alone, on m alone, and on both l and m simultaneously. Most of the recurrences involving m are unstable, and so dangerous for numerical work. The following recurrence on l is, however, stable (compare 5.5.1): (l − m)P m l = x(2l − 1)P m l−1 − (l + m − 1)P m l−2 (6.8.7) It is useful because there is a closed-form expression for the starting value, P m m = (−1)m(2m − 1)!!(1 − x2) m/2 (6.8.8) (The notation n!! denotes the product of all odd integers less than or equal to n.) Using (6.8.7) with l = m + 1, and setting P m m−1 = 0, we find P m m+1 = x(2m + 1)P m m (6.8.9) Equations (6.8.8) and (6.8.9) provide the two starting values required for (6.8.7) for general l. The function that implements this is
254 Chapter 6. Special Functions #include float plgndr(int 1,int m,float x) Computes the associated Legendre polynomial Pm(z).Here m and I are integers satisfying 0≤m≤l,while z lies in the range-1≤x≤l. void nrerror(char error_text[]); float fact,pll,pmm,pmmp1,somx2; int i,ll; if (m 1I fabs(x)>1.0) nrerror("Bad arguments in routine plgndr"); http://www.nr read able files pmm=1.0; Compute pm if(m>0)[ somx2=sqrt(1.0-x)*(1.0+x); fact=1.0; .com or call 1-800-872- (including this one) granted for 19881992 for(1=1;1<=m;i++)( pmm *=-fact*somx2; fact+=2.0; 2 1f(1=m) return pmm; else Compute P+1 pmmp1=x*(2*m+1)*pmm; 1f(1==(血+1)) 7423 (North America to any server computer,is strictly prohibited. tusers to make one paper by Cambridge University Press.Programs from NUMERICAL RECIPES IN C: THE return pmmp1; else Compute Pm,m+1. for(11=m+2;11<1;11+)[ p11=(x*(2*11-1)*pmmp1-(11+m-1)*pmm)/(11-m): only),or pmm=pmmp1; pmmp1=pl1; copy for their return pll; 22 ART OF SCIENTIFIC COMPUTING(ISBN 0-521- CITED REFERENCES AND FURTHER READING: Magnus,W.,and Oberhettinger,F.1949,Formulas and Theorems for the Functions of Mathe- v@cambr matical Physics (New York:Chelsea),pp.54ff.[1] Abramowitz,M.,and Stegun,I.A.1964,Handbook of Mathematical Functions,Applied Mathe- 1988-1992 by Numerical Recipes matics Series,Volume 55 (Washington:National Bureau of Standards;reprinted 1968 by -431085 Dover Publications,New York),Chapter 8.[2] (outside North Software. Amer visit website machine
254 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 machinereadable 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). #include float plgndr(int l, int m, float x) Computes the associated Legendre polynomial P m l (x). Here m and l are integers satisfying 0 ≤ m ≤ l, while x lies in the range −1 ≤ x ≤ 1. { void nrerror(char error_text[]); float fact,pll,pmm,pmmp1,somx2; int i,ll; if (m l || fabs(x) > 1.0) nrerror("Bad arguments in routine plgndr"); pmm=1.0; Compute P mm . if (m > 0) { somx2=sqrt((1.0-x)*(1.0+x)); fact=1.0; for (i=1;im + 1. for (ll=m+2;ll<=l;ll++) { pll=(x*(2*ll-1)*pmmp1-(ll+m-1)*pmm)/(ll-m); pmm=pmmp1; pmmp1=pll; } return pll; } } } CITED REFERENCES AND FURTHER READING: Magnus, W., and Oberhettinger, F. 1949, Formulas and Theorems for the Functions of Mathematical Physics (New York: Chelsea), pp. 54ff. [1] Abramowitz, M., and Stegun, I.A. 1964, Handbook of Mathematical Functions, Applied Mathematics Series, Volume 55 (Washington: National Bureau of Standards; reprinted 1968 by Dover Publications, New York), Chapter 8. [2]
6.9 Fresnel Integrals,Cosine and Sine Integrals 255 6.9 Fresnel Integrals,Cosine and Sine Integrals Fresnel Integrals The two Fresnel integrals are defined by ca==(后到，s回=血(g9m (6.9.1) 8= The most convenient way of evaluating these functions to arbitrary precision is to use power series for small z and a continued fraction for large z.The series are C(x)=x- /T)2x5 /T)4x9 52+(594 (6.9.2) s-(-(（)+( (T3x7 )5x11 (Nort server University 2 11.50-… America Press. THE There is a complex continued fraction that yields both S()and C()simul- taneously: Progra ca+isa=l生af,=9a-0 (6.9.3) 天复 OF SCIENTIFIC( where 61 erfc z 1 11/213/22 2+2+z+2+2十 1920 COMPUTING (ISBN (6.9.4) 22 1 1.2 3·4 222+1-222+5-222+9- In the last line we have converted the "standard"form of the continued fraction to Numerical Recipes 10-621 4310-5 its"even"form (see 85.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 inspection of the (outside terms shows that x =1.5 is a good point to switch over to the continued fraction. Note that for large x North Software. c~+am() (6.9.5) Thus the precision of the routine frenel may be limited by the precision of the library routines for sine and cosine for large x
6.9 Fresnel Integrals, Cosine and Sine Integrals 255 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 machinereadable 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). 6.9 Fresnel Integrals, Cosine and Sine Integrals Fresnel Integrals The two Fresnel integrals are defined by C(x) = x 0 cos π 2 t 2 dt, S(x) = x 0 sin π 2 t 2 dt (6.9.1) The most convenient way of evaluating these functions to arbitrary precision is to use power series for small x and a continued fraction for large x. The series are C(x) = x − π 2 2 x5 5 · 2! + π 2 4 x9 9 · 4! −··· S(x) = π 2 x3 3 · 1! − π 2 3 x7 7 · 3! + π 2 5 x11 11 · 5! −··· (6.9.2) There is a complex continued fraction that yields both S(x) and C(x) simultaneously: C(x) + iS(x) = 1 + i 2 erf z, z = √π 2 (1 − i)x (6.9.3) where ez2 erfc z = 1 √π 1 z + 1/2 z + 1 z + 3/2 z + 2 z + ··· = 2z √π 1 2z2 + 1 − 1 · 2 2z2 + 5 − 3 · 4 2z2 + 9 − ··· (6.9.4) 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.5 is a good point to switch over to the continued fraction. Note that for large x C(x) ∼ 1 2 + 1 πx sin π 2 x2 , S(x) ∼ 1 2 − 1 πx cos π 2 x2 (6.9.5) Thus the precision of the routine frenel may be limited by the precision of the library routines for sine and cosine for large x
