f6-8

246 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).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-tion of variables in spherical coordinates. The spherical harmonic Y lm(θ, φ ), −l≤m≤l,is a functionof the two coordinates θ, φon the surface of a sphere. The spherical harmonics are orthogonal for different landm, and they are normalized so that their integrated square over the sphere is unity: /integraldisplay2π 0dφ/integraldisplay1 −1d(cosθ)Yl/primem/prime*(θ, φ )Ylm(θ, φ )=δl/primelδm/primem (6.8.1 ) Here asterisk denotes complex conjugation. Mathematically, the spherical harmonics are related to associated Legendre polynomials by the equation Ylm(θ, φ )=/radicalBigg 2l+1 4π(l−m)! (l+m)!Pm 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 Pm l(x)=(−1)m(1−x2)m/ 2dm dxmPl(x)( 6.8.4 ) The first few associated Legendre polynomials, and their corresponding nor- malized spherical harmonics, are P0 0(x)= 1 Y00=/radicalBig 1 4π P1 1(x)=−(1−x2)1/2Y11=−/radicalBig 3 8πsinθeiφ P0 1(x)= xY 10=/radicalBig 3 4πcosθ P2 2(x)= 3( 1 −x2) Y22=1 4/radicalBig 15 2πsin2θe2iφ P1 2(x)=−3( 1−x2)1/2xY 21=−/radicalBig 15 8πsinθcosθeiφ P0 2(x)=1 2(3x2−1) Y20=/radicalBig 5 4π(3 2cos2θ−1 2) (6.8.5 ) Thereare manybad ways to evaluateassociated Legendrepolynomialsnumer- ically. For example, there are explicit expressions, such as Pm l(x)=(−1)m(l+m)! 2mm!(l−m)!(1−x2)m/ 2/bracketleftbigg 1−(l−m)(m+l+1 ) 1!(m+1 )/parenleftbigg1−x 2/parenrightbigg +(l−m)(l−m−1)(m+l+1 ) ( m+l+2 ) 2!(m+1 ) ( m+2 )/parenleftbigg1−x 2/parenrightbigg2 −···/bracketrightBigg (6.8.6 ) 6.8SphericalHarmonics 247Sample 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).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 evaluationof 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 lup to 6 or 8, and in double precision (64-bit) for lup 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- ulated in [1-2]. These are recurrences on lalone, on malone, and on both l andmsimultaneously. Most of the recurrences involving mare unstable, and so dangerous for numerical work. The following recurrence on lis, however, stable (compare 5.5.1): (l−m)Pm l=x(2l−1)Pm l−1−(l+m−1)Pm l−2 (6.8.7 ) It is useful because there is a closed-formexpression for the starting value, Pm m=(−1)m(2m−1)!!(1 −x2)m/ 2(6.8.8 ) (The notation n!!denotes the product of all oddintegers less than or equal to n.) Using (6.8.7) with l=m+1, and setting Pm m−1=0, we find Pm m+1=x(2m+1 )Pm 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 FUNCTION plgndr(l,m,x) INTEGER l,mREAL plgndr,x Computes the associated Legendre polynomial P m l(x).H e r e mand lare integers satisfying 0≤m≤l, while xlies in the range −1≤x≤1. INTEGER i,llREAL fact,pll,pmm,pmmp1,somx2 if(m.lt.0.or.m.gt.l.or.abs(x).gt.1.)pause ’bad arguments in plgndr’ pmm=1. Compute P m m. if(m.gt.0) then somx2=sqrt((1.-x)*(1.+x)) fact=1.do 11i=1,m pmm=-pmm*fact*somx2 fact=fact+2. enddo 11 endif if(l.eq.m) then plgndr=pmm else pmmp1=x*(2*m+1)*pmm Compute Pm m+1. if(l.eq.m+1) then plgndr=pmmp1 else Compute Pm l,l>m +1 . do12ll=m+2,l 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 )