Computation of Scattering Kernels in Radiative Transfer Hans Engler Department of Mathematics and Statistics Georgetown University Washington, D.C. 20057, U.S.A. 5 1 January 13, 2015 0 2 n Abstract a phase function p directly leads to expansions for the P . If m J p(x)=(cid:80)∞ α L (x) for all x∈[−1,1], then in particular n=0 n n 0 Thisnoteproposesrapidlyconvergentcomputationalformu- 1 lae for evaluating scattering kernels from radiative transfer (cid:88)∞ theory. The approach used here does not rely on Legen- P0(x,y)= αnLn(x)Ln(y) (3) ] dreexpansions,butratherusesexponentiallyconvergentnu- n=0 A merical integration rules. A closed form for the Henyey- forall−1≤x,y≤1. ThehigherordertermsP , m>0can m N Greenstein scattering kernel in terms of complete elliptic beexpandedintoassociatedLegendrefunctions,usingagain . integrals is also derived. theLegendrecoefficientsofp. Theclassicalreference[2]con- h tains a mathematical presentation of the theory of radiative t Keywords. Plane-parallelscattering,scatteringkernel,phase a transfer. A modern account with more physical details may m function, trapezoidal rule, complete elliptic integral. be found in [6]. [ While the evaluation of eq. (3) and its higher order ver- 1 Introduction and Background sions is in theory straightforward, there are several practical 1 difficulties. Firstly, the formula requires knowledge of the v 7 In theories of radiative transfer and of neutron transport, Legendre expansion of p. While the methods to find p for 0 theinteractionbetweenradiationandascatteringmediumis aparticularscatteringmediumindeedproduceLegendreex- 4 described by a phase function p : [−1,1] → R+ with the pansions(see[11]),itmaybedesirabletohavemorecompact 2 property (cid:82)1 p(t)dt = 1. Let S2 denote the unit sphere representations of a scattering function, in which case the −1 0 in R3. In a single scattering event, a photon or neutron Legendre expansion is not readily available and not easy to 1. that arrives from a given direction Θ ∈ S2 is scattered into compute (there is no “fast Legendre transform”). Secondly, 0 the direction Θ(cid:48) ∈ S2 according to the probability density for cases of strongly forward-peaked scattering, several hun- 5 Θ(cid:48) (cid:55)→ 1 p(Θ·Θ(cid:48)). This probability density then appears dredtermsoftheLegendreexpansionmaybeneededtoeval- 1 as an in4tπegral kernel in an integro-differential equation that uate each Pm even to modest accuracy. This requires care : when evaluating the Legendre polynomials in eq. (3), and it v describes multiple scattering. In the special case of plane- is computationally expensive in any case. i parallelscattering,onerewritesthisequationinsphericalco- X Inthisnote,thedirectnumericalintegrationschemeknown ordinatesandexpandsitasaFourierseriesintheazimuthal as trapezoidal rule is proposed to evaluate P from eq. (2), r angle φ. This requires the evaluation of m a as an alternative to the Legendre expansion in eq. (3). The (cid:90) 2π methodisknown([10])toconvergeexponentiallyiftheinte- c p(cosθcosµ+sinθsinµ·cosφ(cid:48))cosmφ(cid:48)dφ(cid:48) (1) m grandisanalytic. Thishighlydesirablepropertyisexploited 0 systematically in this note. The relation between the con- for 0≤θ, µ≤π, m=0, 1, 2,..., where cm = 21π for m>0 vergence rate and the location of the singularities (points or and c0 = π1. Using the notation x = cosθ, y = cosµ, one regionsofnon-analyticity)ofthephasefunctionisexplained. therefore arrives at the problem of evaluating the scattering The second contribution of this note is the derivation of a kernels Pm(x,y), defined for −1 ≤ x, y ≤ 1 and for m = closedformofthescatteringkernelP0,inthespecialcaseof 0, 1,... as the Henyey-Greenstein phase function ([5]). It expresses P 0 in terms of a complete elliptic integral and can be evaluated c (cid:90) 2πp(xy+(cid:112)1−x2(cid:112)1−y2coss)cosmsds (2) very rapidly without any expansions or numerical integra- m 0 tion, even for cases of extremely forward-peaked scattering. Numerical examples are given to demonstrate the approach. LetL denotethen-thLegendrepolynomial. Fromspherical n harmonics one obtains that the Legendre expansion of the 1 1.1 Henyey-Greenstein Scattering Func- tion The Henyey-Greenstein phase function ([5]) 1 was proposed todescribeinterstellarscatteringandisgivenbytheformula p (x,g)= 1 1−g2 =(cid:88)∞ 2(cid:96)+1g(cid:96)L (x) (4) HG 2(1+g2−2gx)3/2 2 (cid:96) (cid:96)=0 where −1 < g < 1 is known as the asymmetry factor. This phasefunctionhassincebeenusedinareasasdiverseasscat- tering in cloudy and hazy atmospheres ([4]), light scattering in seawater ([3]) and in tissue ([7]), and even in computer graphics. Theneq.(2)togetherwitheq.(3)leadtotheprob- lem of evaluating ∞ H(x,y;g)=(cid:88)2(cid:96)+1g(cid:96)L (x)L (y) (5) 2 (cid:96) (cid:96) (cid:96)=0 for −1≤x,y≤1. The series may be evaluated by using the firstN terms. Since|L (x)L (y)|∼ 2 withindeterminate (cid:96) (cid:96) 2(cid:96)+1 sign, the direct evaluation of the series leads to problems Figure 2: Log Errors for evaluating Henyey-Greenstein when g is close to 1, because then the series converges very scattering kernel in fig. 1. Blue: Eq. (5) with N = 40 slowly. This is illustrated in fig. 1. terms. Black, red, green: Eq. (19) with N =40, 80, 160 terms. Thick: Computed errors. Thin: Error estimates from eq. (20). 2 An Exact Formula We start with the product formula (see 18.17.6 in [8]) for Legendre polynomials 1 (cid:90) π L (cosθ)L (cosµ) = L (cosθcosµ n n π n 0 +sinθsinµcoss)ds. Let ∞ H (x,y;g)=(cid:88)g(cid:96)L (x)L (y). (6) 0 (cid:96) (cid:96) (cid:96)=0 Using the generating function for Legendre polynomials (see 18.12.11 in [8]), we therefore obtain H (cosθ,cosµ;g) (7) 0 Figure 1: Henyey-Greenstein scattering kernel = 1 (cid:90) π(cid:88)∞ gnL (cosθcosµ+sinθsinµcoss) ds (8) H(x,y0;g) for g = 0.95, y0 = 0.4, −1 ≤ x ≤ 1. π 0 k=0 n Black: exact evaluation using eq. (15). Blue: Eq. (5) 1 (cid:90) π ds = (.9) with N =40 terms. Red: Eq. (19) with N =40 terms. π (cid:112)1−2g(cosθcosµ+sinθsinµcoss)+g2 0 By eq. (33) with 1It appears that Chandrasekhar was not aware of this work α = 1+g2−2gcosθcosµ whenhewrotehisclassicaltreatise[2]in1950. β = 2gsinθsinµ α+β = 1+g2−2gcos(θ+µ) this implies the strip S = {z ∈ C||(cid:61)z| < α} and satisfies |f(z)| ≤ M α forsomeconstantM ≥0there. ChooseapositiveintegerN (cid:18) (cid:19) 2 2β H (cosθ,cosµ;g)= √ K (10) and set 0 π α+β 0 α+β 2π (cid:88)N (cid:18)2πk(cid:19) I = f , (19) where K is the complete elliptic integral of the first kind N N N 0 k=1 defined in eq. (28). Setting then u± =1+g2−2gcos(θ±µ) (11) (cid:12)(cid:12)(cid:90) 2πf(t)dt−IN(cid:12)(cid:12)≤ eα4NπM−1 =e−αN1−4πeM−αN . (20) this may also be written as 0 (cid:18) (cid:19) The sum is just an approximation of the integral with the 2 u −u H0(cosθ,cosµ;g)= π√u K0 +u − . composite trapezoidal rule. The rate of convergence thus is + + much faster than for ordinary smooth (twice differentiable) functions, for which the error estimate has the form This formula is closely related to formula (5.10.2.1)2 in [9] which has the form (cid:12)(cid:90) 2π (cid:12) 2π3|f(ξ)| H0(cosθ,cosµ;g)= π(√u 4+√u )K(cid:18)√√uu+−+√√uu−(cid:19) . (cid:12) 0 f(t)dt−IN(cid:12)= 3N2 (21) + − + − for some ξ∈(0,2π). (12) Tousethisresultinthecomputationofeq.(2),notefirst From eq. (10) and eq. (32) we obtain finally that for any α∈R, the function that takes u+i·v=z to H(cosθ,cosµ;g) (13) A+Bcosz=A+Bcosucoshv+Bsinusinhv·i (cid:18) (cid:19) (cid:18) (cid:19) d 1 2 u −u = gdg + 2 π√u K0 +u − (14) maps the strip S to an ellipse about [A−B,A+B] in the + + α (1−g2) (cid:18)u −u (cid:19) complex plane which has focal points A±B and major and = πu √u E0 +u − (15) minoraxeswithlengthsBcoshαandBsinhα. Therefore,if − + + thephasefunctionpisanalyticinaneighborhoodsurround- where E0 is the complete elliptic integral of the second kind ingtheset[−1,1]⊂C,thenhm isanalyticinasuitablestrip defined in eq. (30). In terms of x, y, this becomes Sα with α > 0. The domain of analyticity of p is automat- ically symmetric with respect to the real axis. It should be (cid:32) √ (cid:112) (cid:33) (1−g2) 4g 1−x2 1−y2 emphasized that it is of course not necessary to determine H(x,y;g)= πw √w E0 w (16) this domain in order to use eq. (19). − + + For an illustration, refer to fig. 3. The plot shows the do- with w = 1+g2 −2g(cid:16)xy∓√1−x2(cid:112)1−y2(cid:17). The for- mainofanassumed phase functionthatisoriginallydefined ± on the interval [−1,1] (thin horizontal black line) and that mula may be used to evaluate H(·,·;g) reliably even if g is can be extended into the complex plane (everywhere except extremelycloseto1. Ifg=1−εandε ≈10−16 is“machine 0 at singularities shown as colored lines and circles). For a epsilon”inIEEEarithmetic,thenH(x,y;g)canbeevaluated particular choice of x, y, the integral in eq. (2) extends over to relative accuracy ε /ε. 0 [A−B,A+B]⊂[−1,1] (black circles). A suitable strip S α ismappedtotheellipsesurroundingthissetwherethephase 3 Fast Evaluation function is analytic. Consequently, the trapezoidal approxi- mationconvergesexponentiallywitharategivenbyeq.(20). Wenowturntothegeneralcase. Foragivenphasefunction The rate of convergence depends on α which in turn comes p and given x, y ∈ [−1,1], m ∈ {0, 1,2,...}, we need to from the location of [A−B,A+B] relative to the set of evaluate the integral given by eq. (2). Let singularities of the phase function. h (z) = p(A+Bcosz)cosmz (17) Example 1 Henyey-Greenstein Phase Function. Con- m (cid:112) (cid:112) sider a phase function that is analytic in C minus the ray A = xy, B = 1−x2 1−y2. (18) [1+g2,∞). Then the integrand in eq. (2) is analytic in the 2g Notethat|A±B|≤1,withequalityinoneofthecaseswhere strip defined by x(a[b1=o0u])±ttyth.haetTrihfeeaalfufaunxncictsti,oiontnhhefnmitsihspe2eπtrri-oapdpeiercizooadindicdaloanrnuaRlley.tfiIoctriinaspakpnsrotowrxiinp- cosh(cid:61)z<(cid:12)(cid:12)(cid:12)(cid:12)(1+g2B)/2g−A(cid:12)(cid:12)(cid:12)(cid:12)=(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(√11+−gx2)2/(cid:112)2g1−−xyy2(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . mating the integral (cid:82)2πf(t)dt converges exponentially fast. 0 More precisely, assume that f is 2π-periodic and analytic in This case is illustrated by the magenta line in fig. 3. Similar comments apply to the case where g < 0 (red line). The 2Thiswaspointedouttomebyananonymouscontributoron Henyey-Greensteinphasefunctiondefinedineq.(4)isofthis math.stackexchange.com. Figure 3: Integration domain for a phase function. Figure 4: Multimodal phase function given by eq. (26). form. Theerrorbehaviorofthetrapezoidalrule(19)isillus- wheretheproportionalityconstantsarechosensuchthatthe trated in fig. 2 which shows the log-errors (thick lines) when integrals over [−1,1] equal 1. Both functions have single this rule is used to approximate the integral in eq. (2) for maxima (peaks) at x=x and the width of the peak is pro- 0 x∈[−1,1] and for fixed y=y0, g, for three different choices portional to δ. The Legendre expansions of these functions ofN. Theplotshowsthaterrorsaremaximalforx≈y0 and aregenerallynotavailableinclosedform. Thefunctionf1 is that the errors have the slowest decrease near this value as analyticinthecomplexplaneminusbranchcutsfromx ±iδ 0 N increases. Also plotted (thin lines, same colors) are error tox ±i∞(greenlinesinfig.3). Thefunctionf isanalyticin 0 2 estimatesfromeq.(20),usingM =1forsimplicity. Itisseen thecomplexplaneminuspolesatz=x ±iδ(π/2+nπ), n∈ 0 that the actual errors track the estimate very closely. Z+ (bluecirclesinfig.3). Thereforeifp=f orp=f ,then 1 2 the integrand in eq. (2) is analytic in any strip S where α Example 2 Consider a phase function that is analytic ex- α=arsinhS and S is as in example 2. ceptpossiblyonraysx0±iδ tox0±i∞,wherex0 ∈R, δ>0. Thereadermaynotethattheintegrandineq.(2)contains Then the integrand in eq. (2) is analytic in the strip defined factors cosmz. These terms grow like em|(cid:61)z| away from the by real axis and their second derivatives contain the factor m2. −S 0, x ∈R,γ >0andconsider 0 3.1 Multimodal Phase Functions functions of x∈R (cid:18) (cid:16)x−x (cid:17)2(cid:19)−γ Inpractice,phasefunctionsareobtainedfromscatteringcal- f (x;x ,m,γ) ∝ 1+ 0 (23) culationsusingMietheory,seee.g. [11]. Suchphasefunctions 1 0 δ may have multiple local extrema. An artificial example (not (cid:16)x−x (cid:17) f (x;x ,m) ∝ sech 0 (24) obtained from Mie theory) is given in fig. 4. It uses the 2 0 δ function p(x) = 0.8p (x;.9)+0.1p (x;−.6) (25) HG HG +0.04f (x;.2,.01,3)+0.06f (x;.6,.02).(26) 1 2 wheref andf areasineq.(23,24). Theintegrandineq.(2) 1 2 turns out to be analytic in the strip S with α≈10−2. The α scatteringkernelsP andP forthisphasefunctionwerecom- 0 7 putedat200×200pointswitheq.(19),usingN =128terms ineachcase. AlogarithmicheatmapofP isshowninfig.5 0 andaheatmapofP isshowninfig.6. Thecalculationtook 7 about 10 seconds per scattering kernel on a laptop equipped withadual-coreprocessorrunningat1.40GHz. Therelative accuracy of each result is about 10−3. Figure6: HeatmapofscatteringkernelP forthephase 7 function given by eq. (26). beamswhichmaybecomputedfromscatteringkernels. This is where a fast and accurate computational scheme such as the one presented here will hopefully be of use. 5 Appendix: Complete elliptic integrals Legendre’s complete elliptic integrals K and E of the first 0 0 and second kind are defined as Figure 5: Logarithmic heat map of scattering kernel P 0 for the phase function given by eq. (26). (cid:90) π/2 ds K (m) = (27) 0 (cid:112) 0 1−msin2s 1(cid:90) π ds = (28) 2 (cid:112)1− m ± mcoss 0 2 2 4 Conclusion (cid:90) π/2(cid:112) E (m) = 1−msin2sds (29) 0 0 Adirectnumericalintegrationmethodusingthetrapezoidal 1(cid:90) π(cid:114) m m rule has been presented for the evaluation of scattering ker- = 1− ± cossds. (30) 2 2 2 nels that arise in plane-parallel radiation transfer equations. 0 Its convergence is exponential, and the relation between the where m ∈ C is known as the modulus. Note that usually convergence rate and the domain of analyticity of the phase these integrals are expressed in terms of the parameter k functionisexplained. Thenotealsopresentsaclosedformof where m=k2, leading to the more common notation the scattering kernel for the Henyey-Greenstein phase function,intermsofcompleteellipticintegralsofthesecondkind. K(k)=K (k2), E(k)=E (k2) (31) 0 0 The closed form can be used to assess the accuracy of the proposed numerical integration scheme. For example, the treatment in [8] is in terms of K, E while Mostcomputationalapproachestoplane-parallelradiative Mathematica(cid:13)R uses K0, E0. These integrals converge for transfer use discretizations based on truncated versions of k ∈ C with (cid:60)m < 1 and the functions can be continued eq.(3)(Nystro¨m’smethod). However,someproblemsofthis analyticallytoCminusabranchcutalong[1,∞). Itisknown form also require the evaluation of intensities from scattered (see eq. (19.4.1) in [8]) that [5] Louis G Henyey and Jesse L Greenstein. Diffuse radia- tioninthegalaxy. TheAstrophysicalJournal,93:70–83, d E (m) 2m K (m)= 0 −K (m). (32) 1941. dm 0 1−m 0 [6] Kuo-Nan Liou. An Introduction to Atmospheric Radia- Therefore, for α, β ∈C, we can set m= 2β and obtain in tion, volume 84. Academic press, 2002. α+β the case when (cid:60) 2β <1 [7] Markolf H Niemz. Laser-Tissue Interactions: Funda- α+β mentals and Applications. Springer, 2007. (cid:90) π ds 2 (cid:18) 2β (cid:19) √α−βcoss = √α+βK0 α+β (33) [8] FrankWJOlver.NISTHandbookofMathematicalFunc- 0 tions. Cambridge University Press, 2010. (cid:90) π(cid:112) (cid:112) (cid:18) 2β (cid:19) α−βcossds = 2 α+βE0 α+β (34) [9] Anatoli˘ı Platonovich Prudnikov, Iuri˘ı Aleksandrovich 0 Brychkov, and Oleg Igorevich Marichev. Integrals and where the principal branch of the square root is used. Series: Special Functions, volume 2. CRC Press, 1986. Given real 0 < m < 1, then E (m) and K (m) may be 0 0 [10] Lloyd N Trefethen and J A C Weideman. The expo- evaluated rapidly using the iterations nentially convergent trapezoidal rule. SIAM Review, √ √ 56:385–458, 2014. a = 1, g = 1−m, c = m (35) 0 0 0 a +g √ [11] Warren J Wiscombe. Improved mie scattering algo- a = n n g = a g (36) n+1 2 n+1 n n rithms. Applied Optics, 19(9):1505–1509, 1980. c2 c = n (37) n+1 4a n+1 for n=0, 1, 2,.... Then lim a = lim g =M n n n→∞ n→∞ where M = AGM(1,g ) is known as Gauss’s arithmetic- 0 geometricmean;see[8]. Theconvergenceisquadratic. More- over, lim c =0, and the convergence is also quadratic. n→∞ n One then obtains the values of K and E from 0 0 (cid:32) ∞ (cid:33) K (m)= 2 , E (m)= 2 1−(cid:88)2n−1c2 . (38) 0 πM 0 πM n n=0 Due to the rapid convergence, only a few terms need to be evaluated. An alternative fast computation method for E (M) is given in [1]. 0 Acknowledgement This research was supported by the CooperativeInstituteforResearchintheAtmosphere(CIRA) at Colorado State University. Part of this work was car- ried out at the Joint Center for Satellite Data Assimilation (JCSDA) at NCWCP, College Park, MD. References [1] SemjonAdlaj. Aneloquentformulafortheperimeterof an ellipse. Notices of the AMS, 59(8), 2012. [2] Subrahmanyan Chandrasekhar. Radiative Transfer. Courier Dover Publications, 1960. [3] Vladimir I Haltrin. One-parameter two-term henyey- greenstein phase function for light scattering in seawater. Applied Optics, 41(6):1022–1028, 2002. [4] James Edward Hansen. Exact and approximate solu- tions for multiple scattering by cloudy and hazy plane- taryatmospheres. Journal of the Atmospheric Sciences, 26(3):478–487, 1969.

Computation of Scattering Kernels in Radiative Transfer

Similar Computation of Scattering Kernels in Radiative Transfer