Computation of Scattering Kernels in Radiative Transfer

0.45 MB·English

by  Hans Engler

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 <sinh(cid:61)z<S (22) Whentheintegrandisanalyticandeq.(20)canbeused,the
where S is the unique positive solution of the equation erroristhereforeproportionaltoe−αNeαmM =e−α(N−m)M.
Thustheadditionalfactorcosmzhasthesameeffectasusing
(A−x0)2 + δ2 =B2. N−mpointsinsteadofN pointsfortheevaluationofeq.(19),
S2+1 S2 resulting in a modest loss of accuracy. On the other hand,
iftheintegrand ismerelytwicedifferentiable, theerrorfrom
eq.(21)becomesproportionaltom2/N2. Thustheadditional
To see this, find α such that the ellipse parametrized by factor cosmz now has the same effect as using only N/m
points instead of N points, leading to a much larger loss of
s(cid:55)→A+B(cosscoshα+i·sinssinhα) accuracy. This illustrates the powerful effect of having an
analytic integrand.
passes through the points x ±i·δ, and set S =sinhα.
0
Toobtainexamples,fixδ>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 in-
tegrals
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 func-
tion,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 seawa-
ter. 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.