Angular and spatial moments for half-space albedo transport problems

Annals of Nuclear Energy 86 (2015) 72–79

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

Angular and spatial moments for half-space albedo transport problems Norman J. McCormick Department of Mechanical Engineering, University of Washington, Seattle, WA 98195-2600, USA

a r t i c l e

i n f o

Article history: Received 1 September 2014 Received in revised form 20 December 2014 Accepted 22 December 2014

In memory of Ivan Kušcˇer 1918–2000 and dedicated to M. M. R. Williams.

a b s t r a c t Two classic problems of radiative transfer and neutron transport are solved for a spatially-uniform semiinfinite medium with isotropic scattering. General analytical equations are derived for (1) angular moments of the outward current and (2) spatial moments of the total flux within the half-space. Such moments, for example, can provide analytically explicit equations for determining the surface albedo of the medium as well as the mean depth and mean square distance of travel within the medium. The analysis is done with the Case-style eigenmode method as expressed in terms of the Chandrasekhar Hfunction and its moments. Ó 2014 Elsevier Ltd. All rights reserved.

Keywords: Radiative transfer Neutron transport Angular moments Spatial moments Case’s method Chandrasekhar H-function

1. Introduction We consider the elementary imaging problem of light emerging from a one-dimensional half-space due to uniform illumination over its surface. This problem was elegantly treated in an analytical manner by Chandrasekhar (1960) who obtained the emerging angular distribution in terms of his classic H-function and a few of its angular moments. In the 1960s the eigenmode expansion technique (Case and Zweifel, 1967) was developed that enabled more easily the analytical determination of the light field within, as well as emerging from, the half-space. The eigenmode expansion technique is used here to solve for the angular moments of the outward angular current at the surface and for the spatial moments of the total flux within the half-space for two general incident illumination boundary conditions from which even more general boundary conditions can be constructed. The derivation here illustrates the elegance of employing the eigenmode expansion technique with half-range orthogonality relations for isotropic scattering that contain the weight function lHðlÞ, rather than the weight function used in their original derivation (Case and Zweifel, 1967; Kušcˇer et al., 1964). Moment-driven approximations have a long history in applied transport theory. The explicit expressions for angular and spatial moments in a semi-infinite medium derived here, with the complication of leakage from the medium-vacuum interface of E-mail address: [email protected] http://dx.doi.org/10.1016/j.anucene.2014.12.021 0306-4549/Ó 2014 Elsevier Ltd. All rights reserved.

the half-space, are analogous to those for a localized plane source in an infinite medium with isotropic scattering developed by Ganapol (1985). The infinite medium problem with lethargydependence and anisotropic scattering was solved by Cacuci and Goldstein (1977). Larsen (1996) has developed a moments approach to solving infinite medium 1D, 2D, and 3D problems that includes energy and time dependence; see also an extensive discussion of the moments method for such problems by Sanchez (2003). Recent examples of the use of analytic moments is illustrated in the verification and validation of energy moments for diffusion theory calculations (Crawford and Ring, 2012a; Crawford and Ring, 2012b). When applying approximate analytic searchlight solutions based on the method of images, some limitations of diffusion theory can be avoided by configuring positive and negative source configurations using known properties of the corresponding 1D problem (d’Eon, 2014a). The nonlinear optimization problem in such an approach potentially can be avoided by solving for a source configuration that satisfies either angular or spatial moments derived here. Additionally, zero-variance based variance-reduction schemes for Monte Carlo methods require efficient and accurate approximate solutions to guide random walks. It is likely that moment-driven approximations for the interior solution in a half space could provide improved accuracy while still permitting analytic sampling for 1D and searchlight problems (Krˇivánek and d’Eon, 2014b). Another application of spatial moments is the computation of the mean distance of travel within an idealized, spatially-uniform atmosphere or ocean. The mean depth before a photon is absorbed

73

N.J. McCormick / Annals of Nuclear Energy 86 (2015) 72–79

or escapes through the surface can provide a measure of the euphotic zone, while the mean square distance of travel gives a measure of the average total distance traveled by a photon before either event occurs (Kirk, 1994; Mobley, 1994). The angular moments of the emerging intensity and the spatial moments of the interior intensity also provide a different set of analytical equations for performing benchmark accuracy checks of solutions of idealized half-space albedo problems to supplement those of Ganapol (2008). An advantage of using a set of moments for computing the distribution in the interior, rather than an eigenmode expansion of the analytically-derived angular intensity, is that no singular integrals need to be numerically evaluated. Following a mathematical statement of the objectives of the paper in the next section, the eigenmodes and Chandrasekhar functions are reviewed in Section 3 before deriving the equations for the angular moments of the emerging angular intensity in Section 4. This is followed by equations for obtaining the spatial moments of the total flux in Section 5 and numerical results for both sets of moments in Section 6. For simplicity the treatment of an azimuthally-dependent incident illumination is deferred until Section 7 where the extension to anisotropic scattering also is briefly discussed. 2. Objectives

The emerging moments provide information about the emerging intensity without the need to compute the angular intensity Ið0; lÞ; 0 6 l 6 1. The outward partial current jout;0 divided by the inward partial current jin;0 , for example, is just the albedo of the half-space. From Eqs. (3), (5a) and (5b), the moments of the total angular emerging and incident distributions then are

ðl@=@x þ 1ÞIðx; lÞ ¼ ðc=2Þ

Z

1

Iðx; l0 Þdl0 ;

ð1Þ where c is the single scattering albedo with 0 6 c < 1; x is the distance into the medium in mean free paths, and l is the cosine of the polar angle with respect to the inward surface normal. Both collimated and diffuse boundary conditions

Ið0; l; l‘ Þ ¼ dðl  l‘ Þ

ð2aÞ

Ik ð0; lÞ ¼ lk ;

ð2bÞ

will be considered for l and l‘ 2 ½0; 1. Combination of these boundary conditions gives the quite general surface illumination condition

ð6bÞ

k¼1

‘¼0

where

jout;n ðl‘ Þ ¼ jout;n;k ¼

Z

Z

1

0

lnþ1 Ið0; l; l‘ Þdl

ð7aÞ

1

lnþ1 Ik ð0; lÞdl

ð7bÞ

0

Z

jin;n ðl‘ Þ ¼ Z

1

0

lnþ1 Ið0; l; l‘ Þdl

ð7cÞ

1

lnþ1 Ik ð0; lÞdl:

ð7dÞ

0

After jout;n; total and jin;n; total have been determined for a prescribed set of c‘ and dk values, we can obtain the results for

bn ¼ jout;n;total =jin;0;total ;

x P 0; 1 6 l 6 1;

1

k ¼ 0; 1; . . . ;

L k X X 1 c‘ lnþ1 þ dk ðn þ k þ 2Þ ; ‘

jin;n; total ¼

ð6aÞ

k¼0

‘¼1

jin;n;k ¼

We restrict our attention to isotropic scattering and consider the one-speed radiative transfer equation (Case and Zweifel, 1967)

L K X X c‘ jout;n ðl‘ Þ þ dk jout;n;k

jout;n; total ¼

ð8Þ

which are the fractions of the incident current (i.e., the 0th angular mode) that emerge in the nth angular mode, with b0 the surface albedo (i.e., the ratio of the outward to the inward current) of the semi-infinite medium. When computing bn values it is numerically more convenient, however, to separately obtain analytical equations for each boundary condition by computing the normalized ratios

jratio;n ðl‘ Þ ¼ jout;n ðl‘ Þ=jin;0 ðl‘ Þ

ð9aÞ

jratio;n;k ¼ jout;n;k =jin;0;k

ð9bÞ

The coefficients of the incident illumination, c‘ and dk , are assumed known and, if normalized such that

that lie in the range ½0; 1. The ratios jratio;n ðl‘ Þ and jratio;n;k are the fractions of the incident current that emerge from the surface in the nth angular current mode for the incident illumination conditions of Eqs. (2a) and (2b), with jratio;0 ðl‘ Þ and jratio;0;k the surface albedo for each boundary condition. 1 With jin;0 ðl‘ Þ ¼ l‘ and jin;0;k ¼ ðk þ 2Þ , instead of employing Eq. (6a) we evaluate jout;n;total with the ratios jratio;n ðl‘ Þ and jratio;n;k ,

L K X X c‘ þ dk =ðk þ 1Þ ¼ 1;

jout;n; total ¼

Iin; total ð0; lÞ ¼

L K X X c‘ dðl  l‘ Þ þ dk lk :

ð3Þ

k¼0

‘¼1

ð4Þ

they are the fractions in each of the incident illumination modes, respectively, as defined by Eqs. (2a) and (2b).

ð10Þ

k¼0

‘¼1

k¼0

‘¼1

L K X X 1 c‘ l‘ jratio;n ðl‘ Þ þ dk ðk þ 2Þ jratio;n;k ;

while using Eq. (6b) for jin;n;total . 2.2. Spatial moments of the interior flux

2.1. Angular moments of the emerging intensity Our first objective is to obtain analytical equations for computing the angular moments of the intensity emerging from the half-space,

jout;n ¼

Z

qn ¼

1

l

nþ1

Ið0; lÞdl;

jin;n ¼

n ¼ 0; 1; . . . ;

0

Z

ð5aÞ

0

1

xn qðxÞdx;

n ¼ 0; 1; . . . ;

ð11aÞ

0

where the interior flux is

from knowledge illumination

Z

Our second objective is to obtain the spatial moments of the interior intensity as given by

of

the

known

incident

monochromatic

qðxÞ ¼

Z

1

Iðx; lÞdl:

ð11bÞ

1 1

lnþ1 Ið0; lÞdl; n ¼ 0; 1; . . . :

ð5bÞ

The qn are moments of the distribution with respect to the penetration distance. Such moments provide information about the depth of penetration of the incident illumination that can serve as a

74

N.J. McCormick / Annals of Nuclear Energy 86 (2015) 72–79

guideline for when to terminate particle tracking in a Monte Carlo simulation of the interior distribution or where to place an image source with the method of images. From Eqs. (3), (11a) and (11b), the spatial moments of the total flux are

qn; total ¼

L X

K X c‘ qn ðl‘ Þ þ dk qn;k

ð12Þ

where, for the incident illuminations of Eqs. (2a) and (2b),

qn ðl‘ Þ ¼ qn;k ¼

Z

Z

1 n

dx x

0 1

dx x

n

0

Z

Z

an 

1

1

1 cl ¼1 HðlÞ 2

Z

1

0

Hðl0 Þ dl0 ; l þ l0

ð23Þ

dlIð0; l; l‘ Þ

ð13aÞ

Z

1

ln HðlÞdl; n P 0;

ð24Þ

0

satisfy the following two identities (Chandrasekhar, 1960),

1

dlIk ð0; lÞ:

ð13bÞ

1

1  ca0 =2 ¼ ð1  cÞ1=2 1=2

To determine qn;total we will normalize the results obtained separately for each boundary condition with the ratios defined by n

hx ðl‘ Þi ¼ qn ðl‘ Þ=q0 ðl‘ Þ

ð14aÞ

hxn ik ¼ qn;k =q0;k

ð14bÞ

in order to compute the ratio

hxn itotal ¼ qn;total =q0;total :

ð15Þ

ð1  cÞ

a2 þ ca

2 1 =4

ð25bÞ

Z

1

0

l0 Hðl0 Þ 0 dl ; l R ½1; 0: l þ l0

ð26Þ

Busbridge (1960) derived a set of equations for the higher-order even moments, 2n1 X

ð1Þk ak a2nk ¼ ð2n þ 1Þ1 ;

n P 1;

ð27Þ

k¼1

ð16Þ

k¼0

‘¼1

¼ 1=3:

1 c ¼ ð1  cÞ1=2 þ HðlÞ 2

ð1  cÞ1=2 a2n  ðc=4Þ

L K X X ¼ c‘ q0 ðl‘ Þhxn ðl‘ Þi þ dk q0;k hxn ik :

ð25aÞ

As a consequence of Eq. (25a), the HðlÞ-function also satisfies the equation (Chandrasekhar, 1960)

From Eqs. (12), (14a) and (14b) it follows that

qn;total

l R ½1; 0;

with the constraint 1=Hðm0 Þ ¼ 0. The angular moments of HðlÞ,

k¼0

‘¼1

R where we have introduced the shorthand symbol rþ dm to denote the sum of the contributions from the discrete and continuum eigenmodes. To achieve our objectives the Chandrasekhar HðlÞfunction is needed that satisfies (Chandrasekhar, 1960)

With eigenmodes /ðm; lÞ expðx=mÞ employed in the separation of variables method in Eq. (1), we obtain

that in Appendix A is obtained more directly. (Although these higher-order identities might be of value when testing the numerical consistency of an -values, they are not needed here.) To solve half-space transport problems for the expansion coefficients AðmÞ; m 2 rþ, in Eq. (22), it has been shown (McCormick and Kušcˇer, 1966; McCormick and Kušcˇer, 1973; Kušcˇer and McCormick, 1991) that /ðm; lÞ for m 2 rþ is connected to HðlÞ via a ‘‘half-range-in-l normalization condition,’’

ðm  lÞ/ðm; lÞ ¼ cm=2;

Z

A simplified form of this equation will be obtained later after equations for q0 ðl‘ Þ and q0;k have been derived; see Eq. (59). 3. Eigenmodes and Chandrasekhar functions

1 6 l 6 1;

m 2 r;

ð17Þ

after choosing the ‘‘full-range-in-l normalization condition’’ for the eigenmodes to be

Z

1

/ðm; lÞHðlÞdl ¼ 1;

m 2 r:

ð18Þ

ð1  l=mÞ/ðm; lÞ ¼ ðc=2Þ

Here the eigenvalue spectrum (Case and Zweifel, 1967) for /ðm; lÞ is defined by r  fm 2 ½1; 1 [ m0 g, where m0 is the positive root of the transcendental equation Kðm0 Þ ¼ 0 with

Kðm0 Þ ¼ 1 

c m0 2

Z

1 1

1

m0  l

1

dl ¼ 1  cm0 tanh ð1=m0 Þ:

ð19Þ

1

/ðm; l0 ÞHðl0 Þdl0 ;

m2rþ:

ð29Þ

1

l/ðm; lÞdl ¼ mð1  cÞ; m 2 r;

ð20Þ

1

so we can combine Eqs. (18) and (20) to write 1

l‘ /ðm; lÞdl ¼ m‘ ð1  cÞ‘ ; ‘ ¼ 0; 1; m 2 r:

ð21Þ

1

To analytically solve the half-space problems of interest here with boundary conditions defined for only the half-range, 0 6 l 6 1, only half the eigenmodes are needed for x P 0, i.e., only those for rþ  fm 2 ½0; 1 [ m0 g. Then the eigenmode expansion takes the form

Z

1

Iðx; lÞ ¼ Aðm0 Þ/ðm0 ; lÞ expðx=m0 Þ þ AðmÞ/ðm; lÞ expðx=mÞdm 0 Z  AðmÞ/ðm; lÞ expðx=mÞdm ð22Þ rþ

Also, by multiplying Eq. (29) by HðlÞdl, integrating over 0 6 l 6 1, and taking advantage of Eq. (25a), a second identity directly connecting /ðm; lÞ and HðlÞ is obtained that, when combined with that in Eq. (28), gives

Z

1

l‘ /ðm; lÞHðlÞdl ¼ m‘ ð1  cÞ‘=2 ; ‘ ¼ 0; 1; m 2 r þ :

ð30Þ

0

Integration of Eq. (17) over 1 6 l 6 1 gives

Z

Z 0

1

Z

ð28Þ

that is analogous to Eq. (18) so that Eq. (17) can be written as

1

/ðm; lÞdl ¼ 1;

m 2 rþ;

0

In the same way as for the derivation of the full-range orthogonality and normalization relations (Case and Zweifel, 1967), by writing Eq. (29) for values m and m0 ; m – m0 , the corresponding half-range orthogonality relations follow (McCormick and Kušcˇer, 1966; McCormick and Kušcˇer, 1973; Kušcˇer and McCormick, 1991),

Z

1

/ðm; lÞ/ðm0 ; lÞlHðlÞdl ¼ 0;

m; m0 2 rþ; m – m0 :

ð31aÞ

0

The normalization relations can be written as

Z

1

/ðm; lÞ/ðm0 ; lÞlHðlÞdl ¼ NðmÞHðmÞdðm  m0 Þ;

m; m0 2 ½0;1

ð31bÞ

0

Z

1

/2 ðm0 ; lÞlHðlÞdl ¼ Nðm0 ÞHðm0 Þ:

ð31cÞ

0

(The definitions of NðmÞ for m 2 rþ are not needed for the angular and spatial moments derived here.) The half-range orthogonality

75

N.J. McCormick / Annals of Nuclear Energy 86 (2015) 72–79

relations enable us to obtain the expansion coefficients for Eq. (22) from

AðmÞ ¼ ½NðmÞHðmÞ1

Z

1

Ið0; lÞ/ðm; lÞlHðlÞdl;

m2rþ:

ð32Þ

0

The key equation for deriving the angular moments of the reflected intensity is the half-range reflection relation (or ‘‘albedo operator’’) given by McCormick and Kušcˇer (1966, 1973), Kušcˇer and McCormick (1991)

Ið0; lÞ ¼

1 0

0

0

0

/ðl ; lÞHðlÞHðl ÞIð0; l Þdl ;

0 < l 6 1;

ð33Þ

0

which is sometimes written with the substitution /ðl0 ; lÞ ¼ ðl0 =lÞ/ðl; l0 Þ. Substitution of Eq. (33) into Eq. (5a) gives

jout;n ¼

Z

1

dllnþ1 HðlÞ Z

dl0 /ðl0 ; lÞHðl0 ÞIð0; l0 Þ

1

ð34Þ

nþ1

with the convention here and elsewhere that from Eqs. (35a) and (36) that

k P 0;

P0

j¼1 X j

ð36Þ

f1  HðlÞ½ð1  cÞ

From Eqs. (11a), (11b) and (22) the spatial moments are given

Z

dmAðmÞ

rþ

ð37bÞ

ð38aÞ

j¼1

qn ¼ n!

Z

dm

rþ

l‘ þ ðc=2Þa1   l‘ 2 jratio;2 ðl‘ Þ ¼ l‘  Hðl‘ Þ½ð1  cÞ1=2 l2‘ þ ðc=2Þða1 l‘  a2 Þ:

dx xn expðx=mÞ:

ð41Þ

l and substitution of Eq. (32) gives

mnþ1 NðmÞHðmÞ

Z

1

dl Ið0; lÞ/ðm; lÞlHðlÞ

ð42Þ

0

after use of Eq. (18). After substitution of the boundary conditions in Eqs. (2a) and (2b), respectively, it follows that

qn;k ¼ n!

Z rþ

Z

mnþ1 /ðm; l‘ Þ dm; n P 0 rþ NðmÞHðmÞ mnþ1 Mkþ1 ðmÞ dm; n P 0; k P 0: NðmÞHðmÞ

qn ðl‘ Þ ¼ n!l‘ Hðl‘ Þ

ð43aÞ ð43bÞ

Z

1

lkþ1 /ðm; lÞHðlÞdl; k P 0; m 2 rþ

ð44Þ

Mkþ1 ðmÞ ¼ m½M k ðmÞ  ðc=2Þak ;

k P 1;

m 2 rþ

ð45Þ

with the starting condition M0 ðmÞ ¼ mð1  cÞ1=2 , as given by Eq. (30). Eqs. (43a) and (43b) give the flux moments we wish to analytically express in terms of the H-function and its moments. To derive the flux spatial moments we use the half-spectrum closure relation (McCormick and Kušcˇer, 1973; Kušcˇer and McCormick, 1991; Kušcˇer and Shure, 1967) that proves the completeness of the /ðm; lÞ for m 2 rþ and non-negative l,

Z rþ

/ðm; lÞ/ðm; l‘ Þ dm; NðmÞHðmÞ

0 6 l;

l‘ 6 1:

ð46Þ

ð39aÞ ð39bÞ

Z

Examples of the ratios from Eq. (38a) for the Dirac delta incident illumination include

jratio;1 ðl‘ Þ ¼ Hðl‘ Þ½ð1  cÞ

1

(Equation (46) can be checked by multiplying it by /ðm0 ; lÞlHðlÞdl; m0 2 rþ, and using the orthogonality and normalization relations of Eqs. (31a)–(31c) after integrating the result over ½0; 1.) To obtain the equations for hxn ðl‘ Þi we begin by multiplying Eq. (46) by dl‘ and integrating over ½0; 1 to find, after using Eq. (30),

ð38bÞ

j¼1

jratio;0 ðl‘ Þ ¼ 1  Hðl‘ Þð1  cÞ1=2

Z 0

1

dðl  l‘ Þ ¼ l‘ Hðl‘ Þ

 ð1  cÞ1=2 anþkþ1

k X þ ðc=2Þ ð1Þj aj anþkþ1j :

1=2

dl/ðm; lÞ

and satisfies the recursion relation

after use of Eq. (26). The corresponding angular moments of the 1 inward intensity are just jin;n ðl‘ Þ ¼ lnþ1 and jin;n;k ¼ ðn þ k þ 2Þ ‘ which lead to the sets of ratios jout;n ðl‘ Þ=jin;0 ðl‘ Þ and jout;n;k =jin;0;k we desire for Eq. (10),

1

1

0

j¼1

jratio;n;k ¼ ð1Þk ðk þ 2Þ½ðn þ k þ 2Þ

Z

Here M kþ1 ðmÞ is defined by

0

jratio;n ðl‘ Þ ¼ ð1Þn ln‘ f1  Hðl‘ Þ½ð1  cÞ1=2 n X þ ðc=2Þ ð1Þjþ1 aj lj ‘ g

ð40fÞ

by

1=2

k X þ ðc=2Þ ð1Þjþ1 aj lj gdl

ð40eÞ

5. Spatial moments of the total flux

Mkþ1 ðmÞ ¼

1

ð40dÞ

Again we note that Eqs. (39a) and (40a) for n ¼ 0 give the equations for the albedo of the half-space for the two boundary conditions.

ð37aÞ

j¼1

l

a3 þ ðc=2Þa1 a2  ð1=4Þ jratio;1;2 ¼ 4½ð1=5Þ  ð1  cÞ1=2 a4  ðc=2Þða1 a3  a22 Þ jratio;2;1 ¼ 3½ð1  cÞ1=2 a4 þ ðc=2Þa1 a3  ð1=5Þ: jratio;1;1 ¼ 3½ð1  cÞ

ð35bÞ

 0, it follows

jout;n ðl‘ Þ ¼ ð1Þn lnþ1 f1  Hðl‘ Þ½ð1  cÞ1=2 ‘ n X þ ðc=2Þ ð1Þjþ1 aj lj ‘ g jout;n;k ¼ ð1Þ

ð40cÞ

1=2

Integration over xand

n X lnþ1 l ¼ ð1Þn ln‘ þ ð1Þnþj lj lnj n P 0; ‘ ; l þ l‘ j¼1 l þ l‘

nþkþ1

ð40bÞ

ð35aÞ

after use of Eq. (17) in addition to Eqs. (2a) in (35a) and (2b) in (35b). Because

k

ð40aÞ

jratio;0;2 ¼ 1  4ð1  cÞ1=2 a3

qn ¼

c l HðlÞ l Hðl‘ Þ dl 2 ‘ l þ l‘ 0 Z Z 1 c 1 ðl0 Þkþ1 Hðl0 Þ ¼ dllnþ1 HðlÞ dl0 ; 2 0 l0 þ l 0

Z

n P 0;

jratio;0;1 ¼ 3½ð1  cÞ1=2 a2 þ ðc=2Þa21   1

1

0

jout;n ðl‘ Þ ¼ jout;n;k

Z

0

which leads to

jratio;n;0 ¼ 2½ðn þ 2Þ1  ð1  cÞ1=2 anþ1 ;

with examples for other jratio;n;k values given by

4. Angular moments of the emerging intensity

Z

The diffuse illumination ratios of Eq. (38b) satisfy the especially simple equation for k ¼ 0,

ð39cÞ

rþ

m/ðm; lÞ 1 dm ¼ ; NðmÞHðmÞ ð1  cÞ1=2

ð47Þ

76

N.J. McCormick / Annals of Nuclear Energy 86 (2015) 72–79

so from Eq. (43a) it follows that

l‘ Hðl‘ Þ

q0 ðl‘ Þ ¼

ð1  cÞ1=2

ð48Þ

:

Next, multiply Eq. (46) by

l

n ‘

lnþ1 HðlÞdl and integrate to obtain

Z

/ðm; l‘ ÞM nþ1 ðmÞ ¼ dm; NðmÞHðmÞ rþ

n P 1;

ð49Þ

after taking Eq. (44) into account. Repeated use of Eq. (45) then leads to

M nþ1 ðmÞ ¼ mnþ1 ð1  cÞ1=2  ðc=2Þ

n X

aj mnþ1j ; n P 1; m 2 rþ

ð50Þ

j¼1

after use of Eq. (30). After substitution of the last equation into Eq. (49) and rearranging terms it follows from Eq. (43a) that n cX

1=2

aj qnj ðl‘ Þ ð1  cÞ qn ðl‘ Þ ; n P 1; ¼ ln‘ þ n!l‘ Hðl‘ Þ 2 j¼1 ðn  jÞ!l‘ Hðl‘ Þ "

ln‘ þ

c

n X

aj

2ð1  cÞ1=2

j¼1

ðn  jÞ!

#

hxnj ðl‘ Þi ;

n P 1; ð52Þ

with hx0 ðl‘ Þi  1. This recursion relation leads to

hxðl‘ Þi ¼ l‘ þ 2

2 ‘

hx ðl‘ Þi ¼ 2l þ

c

ð53aÞ (

ð1  cÞ1=2

"

a1 l‘ þ

ca1 2ð1  cÞ1=2

#

) þ a2

ð53bÞ

and so on. The quickest way to derive equations for qn;k in Eq. (43b) is to recognize that qn ðl‘ Þ, obtained using the dðl  l‘ Þ boundary condition, can be viewed as a Green’s function so that

qn;k ¼

Z 0

lk‘ qn ðl‘ Þdl‘ :

"

n! ð1  cÞ1=2

# n cX aj ; anþkþ1 þ q 2 j¼1 ðn  jÞ! nj;k

n P 1;

q0;k ¼ akþ1 =ð1  cÞ1=2

ð55Þ

ð56Þ

obtained similarly with the help of Eq. (48). From repetitive use of Eq. (55) along with Eq. (56) it now follows that

q0;0 ¼ a1 =ð1  cÞ1=2 a2 ca21 q1;0 ¼ þ 1=2

q3;0

2ð1  cÞ

ð1  cÞ 2a3

ð1  cÞ1=2 6a4

þ

2ca1

3ð1  cÞ3=2   6ca1 a3 c 1 3c2 a41 ¼ þ þ þ ca21  2 1=2 1c 16 ð1  cÞ 3 ð1  cÞ

q0;1 ¼ a2 =ð1  cÞ1=2 a3 ca1 a2 q1;1 ¼ þ 1=2 q2;1 ¼

2a3

þ

2c

a1 3ð1  cÞ   6 a4 6ca3 c 1 3c2 a41 2 hx3 i0 ¼ þ þ þ c a  1 a1 ð1  cÞ1=2 ð1  cÞ3=2 a1 3 16 a3 ca1 hxi1 ¼ þ a2 2ð1  cÞ1=2

ð58cÞ ð58dÞ

"

hx2 i1 ¼ a1 2

 # ca1 a3 c 1 ca21 : 2a4 þ þ þ 4 ð1  cÞ1=2 ð1  cÞ 3

ð58bÞ

ð58eÞ

Eqs. (52) to (53b) etc., and (58a)–(58e) provide the hxn iðl‘ Þ and hxn ik needed for computing hxn itotal in Eq. (59). Finally, Eqs. (48) and (56) provide the opportunity to express qn;total in Eq. (16) in the most direct way,

qn;total ¼ ð1  cÞ

1=2

" # L K X X n n c‘ l‘ Hðl‘ Þhx ðl‘ Þi þ dk akþ1 hx ik :

ð59Þ

k¼0

Additional equations for hxn i1 for a surface exposed to an isotropic source that gives an incident illumination of the type Ið0; lÞ ¼ l1 (i.e., for k ¼ 1 in Eq. (2b)) also can be derived, as shown in Appendix B. (The hxn ðl‘ Þi1 results for incident illumination Ið0; lÞ ¼ dðl  l‘ Þ=l‘ that correspond to those in Eq. (52) need not be derived because the constant l‘ cancels out in the ratios in hxn ðl‘ Þi ¼ qn ðl‘ Þ=q0 ðl‘ Þ.) Appendix B also contains results for an incident illumination Ið0; lÞ ¼ /ðm0 ; lÞ, the angular distribution deep within the semi-infinite medium.

ð54Þ

with the starting condition

q2;0 ¼

hx2 i0 ¼

ð58aÞ

1

Substitution of Eq. (51) followed by integration then gives the recursion relation

qn;k ¼

  3ca21 1 þ 4 3ð1  cÞ1=2 a1 1

‘¼1

ca1 2ð1  cÞ1=2

hxi0 ¼

ð51Þ

so after use of q0 ðl‘ Þ from Eq. (48) to obtain hxn ðl‘ Þi ¼ qn ðl‘ Þ=q0 ðl‘ Þ, the final result is

hxn ðl‘ Þi ¼ n!

after use of Eq. (25b) to simplify the results in Eqs. (57c), (57d) and (57g). From Eqs. (57a)–(57g) we see that the ratios hxn ik ¼ qn;k =q0;k are given by

2ð1  cÞ

ð1  cÞ 2a4 ð1  cÞ

1=2

þ

  ca1 a3 ca2 1 ca21 þ þ 3=2 3 1  c ð1  cÞ 4

ð57aÞ ð57bÞ ð57cÞ ð57dÞ ð57eÞ ð57fÞ ð57gÞ

6. Numerical results for angular and spatial moments The computation of the HðlÞ function and an moments has been examined by Busbridge (1960); Das (2008); Williams (2012) and others (Hapke, 1993; Kawabata and Limaye, 2011). All tabulated results presented here were obtained with the numerical results of Bosma and de Rooij (1983). In Tables 1 and 2 for jratio;n ðl‘ Þ ¼ jout;n ðl‘ Þ=jin;0 ðl‘ Þ and jratio;n;k ¼ jout;n;k =jin;0;k are numerical results obtained from Eqs. (39a)–(39c) and Eqs. (40a)–(40f), respectively. The results in Table 1 are consistent with physical reasoning that the angular moments jratio;n ðl‘ Þ diminish with increasing l‘ values and increase with increasing c. It is perhaps surprising, though, how insensitive the values of jratio;n;k in Table 2 are to k for a fixed c and n. The mean of the powers of the distance of travel hxn ðl‘ Þi and hxn ik monotonically increase with increasing n and c and hxn ð0Þi < hxn ik < hxn ð1Þi. Sample numerical results from Eq. (52) are given in Table 3 for hxn ðl‘ Þi for n ¼ 1 to 4 and in Table 4 for hxn ik from Eqs. (58a)–(58e), for n ¼ 1 to 3 and k ¼ 0 and 1. The results for hxn ik for a given n are only weakly dependent on k, with larger percent differences for smaller values of c. If the medium is very weakly absorbing, the numerical values in Tables 3 and 4 for n > 2 are surprisingly large until it is remembered that the asymptotic eigenvalue m0 becomes very large for such media and spatial integrals of the eigenmode Aðm0 Þ/ðm0 ; lÞ expðx=m0 Þ dominate those integrals with eigenmodes in 0 6 m 6 1.

77

N.J. McCormick / Annals of Nuclear Energy 86 (2015) 72–79 Table 1 Angular moment ratios jratio;n ðl0 Þ for isotropic scattering with Ið0; lÞ ¼ dðl  l0 Þ illumination.

Table 3 Mean distances of travel hxn ðl0 Þiin Ið0; lÞ ¼ dðl  l0 Þ illumination.

c

l0

jratio;0 ðl0 Þ

jratio;1 ðl0 Þ

jratio;2 ðl0 Þ

jratio;1 ðl0 Þ

c

l0

0.5

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.24172 0.21266 0.19115 0.17412 0.16014 0.14840 0.13835 0.12963 0.12199 0.11523

0.13762 0.12546 0.11524 0.10657 0.09912 0.09266 0.08700 0.08200 0.07754 0.07355

0.09536 0.08821 0.08183 0.07622 0.07130 0.06695 0.06309 0.05965 0.05656 0.05377

0.33265 0.13340 0.07342 0.04639 0.03175 0.02292 0.01720 0.01329 0.01051 0.00848

0.5

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.31336 0.41337 0.51337 0.61336 0.71336 0.81336 0.91336 1.01336 1.11336 1.21336

0.7

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.39037 0.35231 0.32279 0.29858 0.27813 0.26053 0.24516 0.23160 0.21952 0.20868

0.22535 0.21043 0.19686 0.18476 0.17399 0.16437 0.15574 0.14796 0.14091 0.13450

0.15722 0.14889 0.14062 0.13291 0.12585 0.11941 0.11355 0.10820 0.10331 0.09883

0.87968 0.37067 0.21181 0.13791 0.09678 0.07137 0.05455 0.04283 0.03437 0.02807

0.7

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.9

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.62934 0.59161 0.56002 0.53246 0.50794 0.48585 0.46578 0.44742 0.43054 0.41495

0.37239 0.36131 0.34873 0.33612 0.32393 0.31233 0.30136 0.29103 0.28131 0.27216

0.26312 0.25867 0.25192 0.24442 0.23677 0.22924 0.22194 0.21495 0.20828 0.20193

2.34349 1.06877 0.65099 0.44742 0.32907 0.25291 0.20053 0.16277 0.13457 0.11293

0.9

0.99

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.87751 0.86002 0.84413 0.82925 0.81514 0.80166 0.78874 0.77640 0.76431 0.75272

0.53504 0.53971 0.53930 0.53649 0.53236 0.52745 0.52206 0.51638 0.51052 0.50458

0.38423 0.39229 0.39524 0.39561 0.39445 0.39232 0.38954 0.38633 0.38282 0.37911

4.69507 2.32083 1.51745 1.11221 0.86790 0.70473 0.58824 0.50108 0.43355 0.37981

0.999

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.96077 0.95464 0.94891 0.94341 0.93806 0.93284 0.92772 0.92268 0.91773 0.91285

0.59249 0.60533 0.61220 0.61611 0.61826 0.61925 0.61944 0.61905 0.61823 0.61708

0.42819 0.44261 0.45124 0.45684 0.46057 0.46304 0.46461 0.46553 0.46594 0.46596

5.69253 2.88936 1.93641 1.45310 1.15992 0.96276 0.82095 0.71399 0.63041 0.56330

Table 2 Angular moment ratios jratio;n;k for isotropic scattering with Ið0; lÞ ¼ lk illumination. c

n

k

jratio;n;k

c

n

k

jratio;n;k

c

n

k

jratio;n;k

0.5

0 0 0 1 1 1 2 2 0 0 0 1 1 1 2 2

0 1 2 0 1 2 0 1 0 1 2 0 1 2 0 1

0.14654 0.13657 0.13133 0.09105 0.08572 0.08283 0.06567 0.06213 0.79456 0.78341 0.77694 0.52228 0.51854 0.51601 0.38847 0.38701

0.7

0 0 0 1 1 1 2 2 0 0 0 1 1 1 2 2

0 1 2 0 1 2 0 1 0 1 2 0 1 2 0 1

0.25656 0.24181 0.23390 0.16121 0.15338 0.14904 0.11695 0.11178 0.92971 0.92537 0.92280 0.61691 0.61804 0.61831 0.46140 0.46373

0.9

0 0 0 1 1 1 2 2

0 1 2 0 1 2 0 1

0.47802 0.45977 0.44962 0.30652 0.29725 0.29184 0.22481 0.21888

0.99

0.999

hxðl0 Þi

mfp

for

scattering

with

hx3 ðl0 Þi

hx4 ðl0 Þi

0.44153 0.54420 0.68688 0.86955 1.09222 1.35489 1.65757 2.00024 2.38291 2.80559

1.21069 1.40475 1.69641 2.12169 2.71656 3.51704 4.55912 5.87880 7.51209 9.49498

4.70723 5.34675 6.25865 7.61765 9.65607 12.6639 16.9885 23.0351 31.2665 42.2029

0.53368 0.63368 0.73368 0.83368 0.93368 1.03368 1.13368 1.23368 1.33368 1.43368

1.07259 1.21932 1.40608 1.63279 1.89953 2.20626 2.55300 2.93973 3.09947 3.36647

3.68627 4.09609 4.62994 5.32384 6.21378 7.33576 8.72579 9.79995 10.4199 12.4540

17.3715 19.1739 21.4530 24.4152 28.3246 33.5029 40.3290 49.2406 60.7313 75.3534

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

1.27444 1.37444 1.47444 1.57444 1.67444 1.77444 1.87444 1.97444 2.07444 2.17444

4.63421 4.92910 5.26399 5.63888 6.05377 6.50865 7.00354 7.53842 8.11332 8.72821

26.2453 27.8125 29.5927 31.6217 33.9357 36.5707 39.5625 42.9473 46.7610 51.0397

199.326 211.078 224.339 239.423 256.700 276.598 299.603 326.260 357.168 392.9871

0.99

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

5.18455 5.28455 5.38455 5.48455 5.58455 5.68455 5.78455 5.88455 5.98455 6.08455

1041.27 1059.88 1079.18 1099.19 1119.96 1141.52 1163.90 1187.15 1211.29 1236.37

24143.4 24574.7 25021.9 25485.5 25966.8 26466.5 26985.8 27525.7 28087.5 28672.3

0.999

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

35324.4 35519.9 35717.6 35917.6 36119.7 36324.2 36531.0 36740.2 36951.7 37165.7

2.5808(6) 2.5950(6) 2.6094(6) 2.6241(6) 2.6389(6) 2.6538(6) 2.6689(6) 2.6842(6) 2.6997(6) 2.7153(6)

17.6541 17.7541 17.8541 17.9541 18.0541 18.1541 18.2541 18.3541 18.4541 18.5541

hx2 ðl0 Þi

isotropic

59.8896 60.9665 62.0834 63.2403 64.4372 65.6741 66.9510 68.2679 69.6248 71.0217 644.678 648.248 651.859 655.510 659.201 662.932 666.702 670.513 674.364 678.255

Table 4 Mean distances of travel hxn ik in mfp for isotropic scattering with Ið0; lÞ ¼ lk illumination. c

hxi0

hx2 i0

hx3 i0

hxi1

hx2 i1

0.5 0.7 0.9 0.99 0.999

0.88782 1.11357 1.86442 5.78740 18.2620

1.68449 2.58603 7.05442 67.0858 667.098

4.94806 9.16414 40.1952 1166.61 36553.1

0.96791 1.19150 1.93845 5.85698 18.3299

1.91267 2.84342 7.41095 67.9668 669.652

7. Generalizations Only the case of azimuthally symmetric boundary conditions was treated in the previous sections to keep the presentation as simple as possible. For a transport application that depends on azimuthal angle u through a boundary condition

78

N.J. McCormick / Annals of Nuclear Energy 86 (2015) 72–79

Ið0; l; u; l‘ ; u0 Þ ¼ dðl  l‘ Þdðu  u0 Þ; 0 6 l 6 1, all the previous results are applicable, however, because we have derived equations only for the directionally integrated moments q and j. The generalization to anisotropic scattering is more complicated, however, than the isotropic scattering analysis presented here. The approach of using the Busbridge qn ðlÞ polynomials (Busbridge, 1960) and adjoint eigenmodes (McCormick and Kušcˇer, 1966; Kušcˇer and McCormick, 1991) to compute the emerging angular and spatial moments offers a possible solution.

jratio; n;1 ¼ ðn þ 1Þ½anþ1  ðn þ 2Þ1 ;

Ivan Kušcˇer would have been pleased to join me in honoring our friend Mike Williams whom we both have admired so much. I thank Eugene d’Eon for suggesting this problem and his early interest in numerically checking some of the identities and Ray Aronson for many suggestions that greatly improved the manuscript. Appendices

qn;1 ¼ n!

Here we illustrate a quick way of deriving Eq. (27) for isotropic scattering by generalizing Eq. (23) to the complex z-plane,

1 cz ¼1 HðzÞ 2

Z

1

0

HðlÞ dl zþl

ðA:1Þ

and then doing the power series expansion

H1 ðzÞ ¼ 1  ðc=2Þ

1 X

an =zn ; jzj > 1:

ðA:2Þ

n¼0

Z

¼

Z

cz 2

1

ðA:3Þ jzj > 1;

ðA:4Þ

n¼0

for b2n ¼ ð2n þ 1Þ1 and then matching the coefficients of z2n by employing the Wiener-Hopf factorization (McCormick and Kušcˇer, 1973; Kušcˇer and McCormick, 1991) H1 ðzÞH1 ðzÞ ¼ KðzÞ, where KðzÞ is defined in Eq. (19), to obtain 2n X

a2n  ðc=4Þ

ð1Þk ak a2nk ¼ ð2n þ 1Þ1 ;

n P 0:

ðA:5Þ

k¼0

Appendix B. Moments for other incident illuminations B.1. Angular and spatial moments for

0

Z

l1 ‘ qn ðl‘ Þ dl‘ "

n! ð1  cÞ1=2

n cX aj an þ q 2 j¼1 ðn  jÞ! nj;1

# ðB:5Þ

it follows that

a1

hxi1 ¼

ðB:6aÞ

ð1  cÞ1=2 a0   2 3ca21 ; ¼ 1þ 3ð1  cÞa0 4

ðB:6bÞ

etc., after using Eq. (25a) to derive both results. B.2. Spatial moments for /ðm0 ; lÞ incident illumination Deep within the semi-infinite medium the intensity reaches its asymptotic spatial variation Iðx; lÞ / /ðm0 ; lÞ expðx=m0 Þ that depends only on c which determines m0 , so it is interesting to examine the spatial moments for the surface illumination

Ið0; lÞ ¼ /ðm0 ; lÞ;

0 6 l 6 1:

ðB:7Þ

1

dllnþ1 HðlÞ

0

Z

1

0

dl0

Hðl0 Þ ; l0 þ l

n P 0;

ðB:1Þ

so from Eq. (28) it follows that

jout; n;1 ¼ anþ1  ðn þ 2Þ1

ðB:2Þ 1

qn ðm0 Þ ¼ n!

Z rþ

dm

mnþ1 NðmÞHðmÞ

Z

1

dl/ðm; lÞ/ðm0 ; lÞlHðlÞ:

ðB:8Þ

0

The inner integral vanishes for 0 6 m 6 1, because of the half-range orthogonality relation of Eq. (31a), and equals Nðm0 ÞHðm0 Þ for m ¼ m0 so

qn ðm0 Þ ¼ n!mnþ1 0 :

ðB:9Þ

Thus it follows that the crow’s-flight mean and mean square distances of travel in the semi-infinite medium are

hxðm0 Þi ¼ m0

ðB:10aÞ

hx2 ðm0 Þi ¼ 2m20 ;

ðB:10bÞ

References

l1 incident illumination

Instead of an incident illumination from Eq. (2b) for k P 0, we consider an isotropic external source illumination with k ¼ 1 that provides for a unit current in the inward direction, R1 lIð0; lÞ dl ¼ 1. From Eq. (35b), the equation for the outward 0 current is

c 2

ðB:4Þ

results that are consistent with the fact the surface incident illumination excites none of the eigenmodes in 0 6 m 6 1 so Iðx; lÞ / /ðm0 ; lÞ expðx=m0 Þ for x P 0 (McCormick, 1992).

Eq. (27) follows from this after utilizing Eq. (25a).

jout; n;1 ¼

n P 0:

From Eq. (42) it follows that

1 dl 1 z  l 1 X ¼ 1  ðc=2Þ b2n =z2n ;

KðzÞ ¼ 1 

dm;

1

One also can do a power series expansion of the dispersion function

Z

mnþ1

rþ NðmÞHðmÞ

But with qn ðl‘ Þ from Eq. (55),

hx2 i1

Appendix A. Alpha-moment identities for isotropic scattering

ðB:3Þ

For the spatial moments of the interior fluxes, Eqs. (28) and (43b) with k ¼ 1, lead to

qn;1 ¼

Acknowledgments

n P 0:

and, with jin;n;1 ¼ ðn þ 1Þ , the ratios of the outward to inward angular current moments become

Bosma, P.B., de Rooij, W.A., 1983. Efficient methods to calculate Chandrasekhar’s H-functions. Astron. Astrophys. 126, 283–292. Busbridge, I.W., 1960. The Mathematics of Radiative Transfer. Cambridge, Cambridge, UK. Cacuci, D.G., Goldstein, H., 1977. Neutron slowing down and transport in a medium of constant cross section. I. Spatial moments. J. Math. Phys. 18, 2436– 2444. Case, K.M., Zweifel, P.F., 1967. Linear Transport Theory. Addison-Wesley, Reading, MA. Chandrasekhar, S., 1960. Radiative Transfer. Dover, New York. Crawford, D.S., Ring, T.W., 2012a. Verification of analytic energy moments for the one-dimensional energy dependent neutron diffusion equation with MCNP5 and Attila-7.1.0. Ann. Nucl. Energy. 50, 18–32. Crawford, D.S., Ring, T.W., 2012b. Validation of energy moments from the onedimensional energy dependent neutron diffusion equation, MCNP5 and Attila7.1.0 with the GODIVA experiment. Ann. Nucl. Energy 50, 220–231. Das, R.N., 2008. Application of WienerHopf technique to linear nonhomogeneous integral equations for a new representation of Chandrasekhar’s H-function in

N.J. McCormick / Annals of Nuclear Energy 86 (2015) 72–79 radiative transfer, its existence and uniqueness. J. Math. Anal. Appl. 337, 1366– 1382. d’Eon, E., 2014a. A Dual-beam 3D Searchlight BSSRDF. ACM SIGGRAPH 2014 Talks, Association for Computing Machinery, http://dx.doi.org/10.1145/2614106. 2614140. Ganapol, B.D., 1985. Explicit expressions for the moments of the monoenergetic neutron-transport equation for isotropic scattering in an infinite plane medium. Ann. Nucl. Energy 12, 35–38. Ganapol, B.D., 2008. Analytical Benchmarks for Nuclear Engineering Applications: Case Studies in Neutron Transport Theory. Organisation for Economic Co-operation and Development Nuclear Energy Agency Report No. 6292. Hapke, B., 1993. Theory of Reflectance and Emittance Spectroscopy. Cambridge, Cambridge, UK. Kawabata, K., Limaye, S.S., 2011. Rational approximation formula for Chandrasekhar’s H-function for isotropic scattering. Astron. Astrophys. 332, 365–371. Kirk, J.T.O., 1994. Light & Photosynthesis in Aquatic Ecosystems. 2nd ed., Cambridge, Cambridge, UK. J. Krˇivánek, J., d’Eon, E., 2014b. A Zero-variance-based Sampling Scheme for Monte Carlo Subsurface Scattering. ACM SIGGRAPH 2014 Talks, Association for Computing Machinery, http://dx.doi.org/10.1145/2614106.2614138.

79

Kušcˇer, I., Shure, F., 1967. Closure relations for the eigenfunctions of one-speed transport theory. J. Math. Phys. 8, 823–826. Kušcˇer, I., McCormick, N.J., 1991. Some analytical results for radiative transfer in thick atmospheres. Transport Theory Stat. Phys. 20, 351–381. Kušcˇer, I., McCormick, N.J., Summerfield, G.C., 1964. Orthogonality of Case’s eigenfunctions in one-speed transport theory. Ann. Phys. (NY) 30, 411–421. Larsen, E.W., 1996. How rapidly and how far do neutrons migrate? Ann. Nucl. Energy 23, 341–353. McCormick, N.J., 1992. Asymptotic optical attenuation. Limnol. Oceanogr. 37, 1570– 1578. McCormick, N.J., Kušcˇer, I., 1966. Bi-orthogonality relations for solving half-space transport problems. J. Math. Phys. 7, 2036–2045. McCormick, N.J., Kušcˇer, I., 1973. Singular eigenfunction expansions in neutron transport theory. In: Henley, E.J., Lewins, J. (Eds.), Adv. Nucl. Sci. Technol., vol. 7. Academic, New York, pp. 181–282. Mobley, C.D., 1994. Light and Water Radiative Transfer in Natural Waters. Academic, New York. Sanchez, R., 2003. The moments method for the one-group transport and diffusion equations. Ann. Nucl. Energy 30, 1645–1663. Williams, M.M.R., 2012. A simple method for calculating the moments of the Chandrasekhar H function. Ann. Nucl. Energy 46, 232–233.