A mathematical method for solving mixed problems in multislab radiative transfer

TECHNICAL PAPERS

A mathematical method for solving mixed problems in multislab radiative transfer

M. P. de Abreu

Departamento de Modelagem Computacional; Instituto Politécnico; Universidade do Estado do Rio de Janeiro; Caixa Postal 97282; 28601-970 Nova Friburgo, RJ. Brazil; [e-mail]

ABSTRACT

In this article, we describe a mathematical method for solving both conservative and non-conservative radiative heat transfer problems defined on a multislab domain, which is irradiated from one side with a beam of radiation. We assume here that the incident beam may have a monodirectional (singular) component and a continuously distributed (regular) component in angle. The key to the method is a Chandrasekhar decomposition of the (mathematical) multislab problem into an uncollided transport problem with singular boundary conditions and a diffusive transport problem with regular boundary conditions. Solution to the uncollided problem is straightforward, but solution to the diffusive problem is not so. For then we make use of a recently developed discrete ordinates method to get an angularly continuous approximation to the solution of the diffusive problem. We suitably compose uncollided and diffuse solutions, and the task of generating an approximate solution to the original multislab radiative transfer problem is complete. We illustrate the accuracy of the proposed method with numerical results for a test problem in shortwave atmospheric radiation, and we conclude this article with a discussion.

Keywords: Radiative transfer, multislab problems, mixed beams, conservative scattering, discrete ordinates

Introduction

We have been long working on the development of analytical and numerical methods for solving basic and important problems in neutron transport theory. These methods are intended to yield accurate solutions to discrete ordinates (S_N) versions of the integrodifferential neutron transport equation (Davison, 1957; Duderstadt and Martin, 1979; Lewis and Miller Jr., 1993). Roughly speaking, an S_N version of a transport equation consists of taking the integrodifferential transport equation in a finite set of angular directions (discrete ordinates), and replacing therewith the integral source term by a suitable quadrature formula. As a result, we obtain a coupled system of linear differential equations of the first order whose unknowns are approximations to the angular density of neutrons in the discrete directions considered. This system of linear differential equations is referred to as the S_N equations (Lewis and Miller Jr., 1993). Discrete ordinates approximations to the boundary conditions can be likewise derived, and an S_N problem is said to be stated (S_N equations plus S_N boundary conditions).

The neutron transport problems that can be solved with our current S_N methods can be divided into two major classes: multiplying and non-multiplying problems. The first class incorporates neutron transport problems defined in critical systems with respect to neutron fission chain reactions (de Abreu, Alves Filho and Barros, 1993, 1996; de Abreu and Barros, 1994; de Abreu, 1995, 1997). A typical problem in this class is that of computing thermal power distribution and neutron multiplication factor in nuclear reactor design and analysis (Duderstadt and Hamilton, 1976). The second class addresses those problems where neutron regeneration is only due to scattering events, and the boundary conditions are modelled by smooth functions of the angular variable (Barros and Larsen, 1990, 1991; de Abreu et al., 1998; de Abreu, 1998, 2001a, 2001b, 2002). A classical problem here is that of computing the emerging space-energy-angle distribution of neutrons due to a highly active neutron source in radiation shielding design (Shultis and Faw, 1996).

As noticed by S. Chandrasekhar (1950) more than a half century ago, the mathematical problems that arise in neutron transport theory and radiative heat transfer are essentially the same. For example, some problems in the second class above and some problems in radiative transfer in planetary atmospheres (Chandrasekhar, 1950; Liou, 2002) and in vegetation canopies (Myneni and Ross, 1991; Ganapol et al., 1998, 1999) are formulated by essentially the same transport equation, with the differences lying only in the integral source terms. These matters have recently attracted our attention to the point that we have been dedicating part of our work with transport methods to the extension of some of our slab-geometry S_N analytical and numerical methods to solve basic and important problems in the theory of radiative transfer. Unfortunately, things are not so straightforward. For example, the boundary conditions for radiative transfer problems are often more complicated than those in neutron transport due to the presence of interacting (absorbing-emitting-scattering) boundaries (Stamnes et al., 1988; Myneni and Ross, 1991; Liou, 2002). Also, basic and important problems in radiative transfer are characterized by highly anisotropic scattering (Garcia and Siewert, 1985; Liou, 2002). In this case, the phase function of the scattering angle is often approximated by a high order polynomial expansion. And, particularly in atmospheric radiative transfer, we have most often to solve conservative problems, i.e., those associated to the conservative case of perfect scattering of radiation (Chandrasekhar, 1950). Conservative scattering is known to introduce additional difficulties into the mathematical analysis of radiative transfer problems (Shultis and Hill, 1976; Benassi et al., 1984; Siewert, 2000) and, for the lucky neutron transport analyst, it is much less frequent in neutron transport problems. On one hand, these differences make most radiative transfer problems more difficult to analyse and solve than similar neutron transport ones but, on another hand, they indicate the room left and directions for improvements in our S_N analytical and numerical methods for solving radiative transfer problems of theoretical and practical interest. Accordingly, we have been working on improvements in our S_N methods on a stepwise and timely basis. As part of our continued effort on these matters, we have recently developed a mathematical method based on improvements in our previous S_N analytical and numerical methods to yield approximate solutions to both conservative and non-conservative problems in radiative transfer defined on a multislab domain with transparent (non-interacting) boundaries, which is irradiated from one side with a monodirectional beam of radiation (de Abreu, 2004). To show the method's accuracy, we solved basic and important problems in atmospheric radiative transfer such as the four azimuthally symmetric test problems posed in 1977 by the Radiation Commission of the International Association of Meteorology and Atmospheric Physics (Lenoble, 1985), as well as some multislab problems adapted from a six-layer model atmosphere discussed in a work of Devaux et al. (1979) on the facile (F_N) method.

We find it appropriate at this point to make things clear as to what we mean when we speak of a multislab conservative problem and of a multislab non-conservative problem in radiative transfer. A standing point in the theory of radiative transfer is the formulation (and subsequent analysis) of conservative and non-conservative problems defined in a homogeneous plane-parallel medium. These nomenclatures are direct reference to the conservative case of perfect scattering of a pencil of radiation by a mass element in the medium (Chandrasekhar, 1950; Thomas and Stamnes, 1999). In the conservative case of perfect scattering, the radiant energy removed from a pencil of radiation traversing a mass element in a given direction of propagation is entirely transferred to other directions as scattered radiation of the same energy. Otherwise, part of the radiant energy removed is transformed into other forms of energy or even into scattered radiation of other energies. The fraction of the energy removed from the pencil of radiation, which appears as scattered radiation of the same energy in all directions, is the single scattering albedo (v) of the homogeneous medium. From the above discussion, it is apparent that if v = 1 then, the homogeneous medium is conservative in the sense described above, and a radiative transfer problem defined therein is said to be conservative. If otherwise 0 < v < 1, then the homogeneous medium truly absorbs (Chandrasekhar, 1950) radiant energy and a radiative transfer problem defined therein is now said to be non-conservative. These definitions are appropriate to problems defined in homogeneous media but we do not think the same for the problems dealt with in this article – multislab problems. For it may well happen that v = 1 for a given layer in the multislab domain and that 0 < v < 1 for another one, and according to the Principle of Non-Contradiction (Russell, 1976), a radiative transfer problem cannot be conservative and non-conservative at the same time. We shall therefore have to give a definition for a conservative problem defined on a multislab domain. Since a multislab is a stratified plane-parallel non-homogeneous region, it seems to be appropriate to define a multislab conservative problem as follows: a radiative transfer problem defined on a multislab domain W is conservative if v(t) = 1 for some t on W. And denying our condition for conservativeness, i.e. 0 < v(t) < 1 everywhere in W provides a definition for a multislab non-conservative problem. Though unexpected from a physical standpoint, these definitions proved to be useful as regarding computational efficiency. We shall return to this point in a section ahead.

In this article, we describe a mathematical method for approximately solving both conservative and non-conservative problems in radiative transfer defined on a multislab domain. We assume here that the multislab domain has transparent boundaries. We assume further that the multislab domain is irradiated from one side with a beam of radiation having a monodirectional (singular) component and/or a continuously distributed (regular) component in angle. This makes the method here a bit more general, for it enables us to solve problems covered by the method recently developed by the present author (de Abreu, 2004) and problems not covered so far (with mixed beams). To our knowledge, a mathematical method with the features presented here is not found in the specialized literature.

An outline of the remainder of this article follows. Right after the Nomenclature section, we formulate and perform an analysis of the target problem that represents the class of radiative transfer problems dealt with in this article. We then develop a mathematical method for accurately solving the target problem. We next discuss computational aspects of the method, and we present numerical results for a test problem. In the last section, we conclude this article with closing remarks, and we report ongoing research.

Nomenclature

a = angular component, dimensionless

A = auxiliary operator, dimensionless

c = source vector, W/(m² sr)

d = source vector, W/(m² sr)

E = real square matrix, dimensionless

F = real square matrix, dimensionless

f = constant in the particular solution component of the intensity, W/(m² sr)

g = Chandrasekhar polynomial, dimensionless, and coefficient in the ESGF equations, W/(m² sr)

H = unit step function, dimensionless

I = identity matrix, dimensionless

I = frequency-integrated intensity of the radiation field, W/(m² sr)

K = positive real number, dimensionless

L = order of Legendre expansion, dimensionless

M = diagonal matrix, dimensionless

N = order of quadrature set, dimensionless

P = Legendre polynomial, dimensionless

q = radiative heat flux, W/m2

R = number of layers, dimensionless

S = scattering source, W/(m² sr)

T = real square matrix, dimensionless

Greek Symbols

a = expansion coefficient in the homogeneous solution component of the intensity, W/(m² sr)

b = factor in a Legendre component of the scattering phase function, dimensionless

g = boundary function for the intensity, W/(m² sr)

Dt = optical thickness, dimensionless

d = Dirac distribution, dimensionless, and Kronecker discrete function, dimensionless

q = coefficient in the ESGF equations, dimensionless

µ = cosine of the polar angle, dimensionless

n = separation constant, dimensionless

r = reflectance of an interacting boundary, dimensionless

t = optical depth, dimensionless

f = angular moment of the intensity, W/m2

W = multislab domain, dimensionless

w = angular weight, dimensionless

v = single scattering albedo, dimensionless

Subscripts

d relative to diffuse reflection

i relative to ordered set

j relative to layer edge and ordered set

l relative to a component in a Legendre expansion

m relative to discrete direction

N relative to order of quadrature set

n relative to discrete direction

p relative to particular solution

R relative to right boundary and rightmost layer

r relative to layer number and layer edge

s relative to specular reflection

t relative to discrete direction

u relative to discrete direction

Superscripts

d relative to diffusive problem

r relative to layer number

T relative to transpose

u relative to uncollided

0 relative to left boundary and zeroth order

+ relative to positive directions and upward heat flux

– relative to negative directions and downward heat flux

Formulation and Analysis of the Target Problem

We start this section with a mathematical formulation of the target problem representing a class of radiative transfer problems with anisotropic scattering defined on a multislab domain irradiated from one side with a beam of radiation. Let us consider the equation of transfer with azimuthal symmetry and anisotropic scattering of the form

where t is the optical variable defined on a multislab domain W with transparent boundaries denoted by t₀ (left) and t_R (right), respectively; µ is the cosine of the polar angle defined by the direction of the propagating radiation and the positive t-axis. The quantity I(t,µ) is the frequency-integrated intensity of the radiation field in direction m at optical depth t and S(t,µ) is the scattering source function given by

The quantity v(t) is the single scattering albedo at depth t; (2l+1)b_l(t) is the lth-order component of the Legendre expansion of the scattering phase function and P_l(µ) denotes the lth-degree Legendre polynomial. We assume that the multislab domain W consists of R contiguous and disjoint layers of homogeneous material each, i.e. the quantities v(t) and b_l(t) for all l are piecewise constant functions of t on W. In accordance with our definitions for multislab conservative and non-conservative problems given in the introductory section, if v(t) = 1 somewhere in W, then the multislab problem is conservative. Otherwise, the problem is non-conservative. The transfer Eq. (1) is subject to the boundary conditions

where I₀ is a nonnegative real; µ₀ is the cosine of the polar angle defining the direction of incidence of the monodirectional component of the beam of radiation upon the left boundary of the multislab domain W; the symbol d is to denote a Dirac distribution and g₀(µ), µ > 0, is a nonnegative function of µ representing the angularly continuous component of the incident beam of radiation. Equations (1-3) define the (mathematical) target problem representing the class of radiative transfer problems dealt with in this article.

Following a decomposition technique considered by Chandrasekhar (1950) in solving a basic problem in radiative transfer in planetary atmospheres, we decompose the target problem (1-3) into the uncollided problem

with the left singular boundary conditions

and the diffusive problem

with the regular boundary conditions

so that

The quantity

in Eq. (6) is a depth-dependent anisotropic source given in terms of the solution I^u(t,m) to the uncollided problem (4-5).

We perform an analysis of both the uncollided problem and the diffusive problem in order to get analytical results important to define and support our mathematical method. We begin with the uncollided problem (4-5). Since the uncollided problem (4-5) represents an auxiliary problem defined in a purely absorbing domain with transparent boundaries, with no source of radiation and with an incident beam upon the left (t₀) boundary only, we must have I^u(t,µ) = 0, for µ < 0 and for all t Î W (Case and Zweifel, 1967; Liou, 2002). For m > 0 and t Î W, we have I^u(t,µ) > 0, and we can write the uncollided Eq. (4) in the integral form

We solve Eq. (9) for I^u(t,µ) and we successively obtain

and

We substitute the left singular boundary condition in (5) into Eq. (12) and the uncollided problem (4-5) has the closed form solution

At this point, we may substitute the closed form solution (13) into the depth-dependent anisotropic source (8) to completely define the more challenging problem – the diffusive problem (6-7). So, we substitute solution (13) into the source (8) to obtain

We now perform an analysis of the diffusive problem (6-7). We decompose the multislab domain W into its R contiguous and disjoint homogeneous subdomains (layers) and we define the local (layer-level) diffusive equations

with

and with intensity continuity conditions at layer interfaces, i.e.

where t_j , j = 1 : R-1, is to denote the jth layer interface. We remark that if R = 1, then the multislab domain W consists of one single layer. We note further that

If v(t) happens to be equal to 1 for some t in the interval [t_r-1 ,t_r] then, v(t) = 1 everywhere in [t_r-1 ,t_r], and therefore v_r = 1. In this sense, we may think of a "conservative layer", and we are giving an alternative definition (perhaps more adequate at this point) for a multislab conservative problem: a multislab problem is conservative if at least one layer is conservative. It is non-conservative otherwise. Continuing, the depth-dependent anisotropic source that appears in Eq. (15) can be expressed as

Result (17) can be written in the shorter form

where

We follow by considering standard S_N approximations (Lewis and Miller Jr., 1993) to the local diffusive equations (15) of the form

where

and {µ_m}, m = 1 : N, is a finite set of angular directions on the interval [-1,1]. In this article, we consider even-order Gauss-Legendre quadrature sets (Lewis and Miller Jr., 1993). That is to say, the directions µ_m, m = 1 : N, are the roots of the Nth-degree Legendre polynomial P_N(µ), and the angular weights w_n , n = 1 : N, are so that the quadrature formula on the right side of the equal sign in Eqs. (20) integrates exactly Legendre polynomials from P₀(µ) to P_2N-1(µ). We remark that we have ordered the directions m_m so that µ_m > 0 holds for m = 1 : N/2, µ_m < 0 holds for m = N/2+1 : N, where µ_m-1 < µ_m, m = 2 : N/2 and µ_m+N/2 = –µ_m, m = 1 : N/2. We remark also that the nonnegative integers L_r in Eqs. (20), r = 1 : R, indicate that the Legendre expansions of the scattering phase functions in corresponding layers have been truncated after (L_r + 1) terms. We note further that intensity continuity conditions still hold at layer interfaces, and that we have used the standard discrete boundary conditions (Lewis and Miller Jr., 1993)

where the subscript –m is to denote the angular direction –µ_m.

Solution to the S_N system (20) can be expressed for each r in terms of the solution to the homogeneous version of Eqs. (20), and a particular solution in the vector form

where

a_r,i , i = 1 : N, are scalars depending upon the discrete boundary conditions (21);

are the elements of a vector basis for the null space of the local S_N radiative transfer operator

where (·) denotes appropriate operations on the entries of vector (23), and

Owing to the translational invariance (Case and Zweifel, 1967) of the homogeneous version of the S_N Eqs. (20), exponential functions are the natural choice for the entries of vector (23). These are the well-known elementary solutions (Chandrasekhar, 1950; Case and Zweifel, 1967) to the homogeneous version of the S_N Eqs. (20), and they are expressed here in the form

where t_r,i, i = 1 : N, are appropriate depths. The quantities n_r,i and a_r,m(n_r,i), are the separation constants and the angular components (de Abreu, 1998; Siewert, 2000) of the elementary solutions (26), respectively.

We now show that the natural choice (26) implies that the vector basis with elements (23) is a basis of eigenvectors of an auxiliary local S_N operator. For so we substitute the exponential solutions (26) into the homogeneous version of the S_N Eqs. (20), we first order operate on the left side of the equal sign and we obtain

Equations (27) can be reformulated further and be cast in the matrix form

or, equivalently,

where the operator A_r on the left side of the equal sign in Eq. (29) is an auxiliary S_N (matrix) operator with entries

where d_mn is the Kronecker delta and

Equation (29) is an eigenvalue problem defined by the auxiliary local S_N operator A_r. The eigenvalues (n_r,i)^-1, i = 1 : N, of operator A_r are the reciprocals of the separation constants n_r,i and the corresponding eigenvectors are the basis vectors (23) with entries given by (26). Therefore, the vector basis for the null space of operator (24) with elements (23) and entries (26) is a basis of eigenvectors of the auxiliary local S_N operator A_r. But nothing has been said until now about the existence of such basis of eigenvectors. We discuss this very important point next.

For an oriented discussion of the existence for the null space of operator (24) of a basis of eigenvectors of the auxiliary operator A_r, we shall step back a bit to state some preliminary results and to make a relevant assumption. We begin by dividing Eqs. (27) by the function exp [(t – t_r,i) / n_r,i], and we notice that the resulting -dependent, N-term series on the right side of the equal sign is just the zeroth-order th-degree Chandrasekhar polynomial (1950) evaluated at n_r,i. Then, we may write

Now, with the parity relations

for Legendre (Lewis and Miller Jr., 1993) and Chandrasekhar (1950) polynomials, respectively, it is not difficult to draw that the separation constants appear in ± pairs of numbers and that the angular components satisfy the relation a_r,m(n_r,i) = a_r,–m(–n_r,i) for all r, m and i. In this article, we assume that the ± pairs of separation constants are distinct and real. The relationship for the angular components will be brought to further discussion only in the final section of this article. Let us make use of the former result (± pairs of numbers) and of our assumption (distinct real) to discuss the existence of a basis of eigenvectors of operator A_r. If the separation constants are bounded, i.e. there exists a positive real number K so that |n_r,i| < K for all i, then zero is not at all an eigenvalue of A_r, viz Eq. (29). Since the separation constants are distinct pairs of ± real numbers and zero is not an eigenvalue of A_r, the eigenvalues of A_r form a distinct set of real numbers. This implies that there exists a basis of eigenvectors of A_r (Schneider and Barker, 1989). Numerical evidence (Shultis and Hill, 1976; Siewert, 1978, 2000) indicates that this is true of non-conservative layers (0 < v_r < 1). On the other hand, it is known (Chandrasekhar, 1950; Case and Zweifel, 1967; Benassi et al., 1984; Siewert, 2000) that, for conservative layers, there exists a ± pair of unbounded constants satisfying for example Eqs. (27). Therefore, the auxiliary local S_N operator A_r has a zero eigenvalue with algebraic multiplicity (Schneider and Barker, 1989) equal to 2 (a double zero eigenvalue), viz Eq. (29). Unfortunately, the operator A_r has no special properties (to our knowledge) to guarantee that the geometric multiplicity (Schneider and Barker, 1989) of the double zero eigenvalue be also equal to 2. So, for the conservative case v_r = 1, the operator A_r may not provide a basis of eigenvectors for the null space of the local S_N operator (24). But even so we can do the following: when generating the eigenvalues and corresponding eigenvectors of the auxiliary local S_N operator A_r and for v_r = 1, we simply discard the eigenvector(s) corresponding to the double zero eigenvalue and replace it (them) with a pair of vectors of the form (23) belonging to the null space of the local S_N operator (24). These replacing vectors have entries given by the elementary solutions

and

respectively, with Dt_rº t_r – t_r-1 and |b_1,r| < 1 for all r. The elementary solutions (33-34) are S_N variants to elementary solutions of the homogeneous integrodifferential transfer equation reported on page 14 of the book of Chandrasekhar (1950) for the conservative case (the reader may check solutions (33-34) against the homogeneous version of the S_N Eqs. (20)). Similar solutions are also reported in the Appendix F of Case and Zweifel (1967) for neutron transport problems but for isotropic scattering (b_1,r = 0 in this case). Elementary solutions very similar to (33-34) can also be found in a recent work by Siewert (2000).

In this article, the entries of the particular solution vector (25) are given by the exponential functions

The exponential functions (35) are S_N versions of the particular solution used by Chandrasekhar (1950) to express the solution to the classical albedo problem in radiative transfer, and f_r,m , m = 1 : N, are constants to be determined upon substitution of functions (35) into Eqs. (20). As noted by Siewert (2000) in a recent work, the exponential functions (35) are not valid in the (unlikely) event that µ₀ be equal to one of the separation constants n_r,i , i = 1 : N. So, we are assuming here that µ₀ does not match any n_r,i for all i. The constants f_r,m, m = 1 : N, in the exponential functions (35) can be efficiently computed with the help of the parity relation for the Legendre polynomials and through a not-so-straightforward use of matrix algebra for splitting and handling the matrix equation resulting from the substitution of the exponential functions (35) into the S_N Eqs. (20). An outline of this technique follows. For fixed r, the constants f_r,m , m = 1 : N, are the entries of the (N/2)-dimensional column matrices

and

where

and

respectively. The quantities

are (N/2)-dimensional real square matrices. The symbol I is to denote the (N/2)-dimensional identity matrix and M is the (N/2)-dimensional diagonal matrix whose entries are M_ii = µ_i, i = 1 : N/2. The quantities and denote (N/2)-dimensional real square matrices with entries

and

respectively. The quantities c_r and d_r in results (36-37) are (N/2)-dimensional column matrices given by

and

respectively.

From the foregoing analysis, we may draw a number of conclusions important to define and substantiate the method to be described in the next section. Firstly, if I₀ = 0 and the boundary function g₀(µ) ¹ 0, then we draw from result (13) that I^u(t,µ) = 0 everywhere in W and for all µ, from result (14) that s^u(t,µ) = 0 everywhere in W and for all µ and from results (36) through (43) that the constants f_r,m = 0, for all r and for m = 1 : N. Therefore, the exponential functions (35) are zero everywhere in the real interval [t_r-1 ,t_r] for all r and the particular solution vector (25) is the null vector of the null space of operator (24). In this case, solution to the target problem (1-3) in each of the discrete directions in the finite set {µ_m}, m = 1 : N, is approximately given by the open form solution (22). This should not be surprising, for the resulting target problem is just the diffusive problem (6-7) with s^u(t,µ) = 0, which has been formulated and approximately solved with our previous S_N analytical and numerical methods. On another hand, if I₀¹ 0 and the boundary function g₀(µ) = 0 then, formulation and analysis provided here are mathematically equivalent to those recently reported by us (de Abreu, 2004) for problems with only a monodirectional beam. Therefore, a mathematical method founded in the results derived in the foregoing analysis may be expected to work with both one-type (singular or regular) and mixed boundary conditions. Secondly, the scalars a_r,i are hitherto unknown quantities. This is a call for one of the open questions in our preceding analysis – the determination of the scalars a_r,i in all layers of the multislab domain W, and we intend to address an appropriate reply in the coming section. Thirdly, even with the scalars a_r,i on hand, the best we can do with available information is to close up the open form solution (22) and add it to the closed solution (13) to yield an approximate solution to the target problem (1-3) only in the discrete directions considered in the S_N Eqs. (20). In the preceding analysis, nothing has been said about a way of obtaining approximate solutions to the local diffusive equations (15) that are continuous in t and µ. This gap is bridged with the equivalence between the S_N formulation considered in this article and an angularly continuous formulation classical in radiation transport theory.

At this point, we proceed to describe a mathematical method for an approximate solution to the target problem (1-3).

A Mathematical Method for Solving the Target Problem

The mathematical method described in this section is a conjugation of basic relations from more general results in the theory of radiation transport and a two-component method based on recently developed S_N analytical and numerical methods for an approximate solution to the target problem (1-3). The approximate solution we seek is a distribution on t and µ of the form

where I^u(t,µ) means that we have directly incorporated the closed form (13) into our approximate solution, and the second term on the right side of the equal sign in (44) denotes the well-known spherical harmonics (P_N-1) approximation (Duderstadt and Martin, 1979; Lewis and Miller Jr., 1993) to the solution of the diffusive Eqs. (15) given by

where

The quantities

are the P_N-1 angular moments (Duderstadt and Martin, 1979; Lewis and Miller Jr., 1993) of the diffuse component of the intensity of the radiation field. Results (45) and (46) can be shown to come up from two different and equivalent formulations of the diffusive problem (6-7) – the S_N formulation composed by the S_N Eqs. (20) and the discrete boundary conditions (21) and the classical P_N-1 formulation with the corresponding boundary conditions due to Mark (Davison, 1957; Duderstadt and Martin, 1979). Expressions (44-46) constitute the basic relations of the mathematical method described in this section. We should note that the diffuse component (45) of our approximate solution (44) is continuous in t and µ on W. We note further that the only unknown quantities in expressions (44-46) are the entries of the S_N solution vector (22), which will be provided by a two-component method based on recently developed S_N analytical and numerical methods.

As the name implies, our two-component method has two ingredients: a numerical component and an analytical component. The numerical component is to provide layer-average

and layer-edge values for the entries of the S_N solution vector (22) without having to determine the scalars a_r,i, r = 1 : R, i = 1 : N, in the open form (22). The numerical component is thus suited to radiative transfer problems where the quantities of interest are, for example, the angular distribution of radiation leaving the multislab domain and angle-integrated layer-edge quantities such as radiative heat fluxes (Chandrasekhar, 1950; Thomas and Stamnes, 1999). For the angular distribution of leaving radiation, we make direct use of results (46) through (44) at t = t₀ for µ < 0 and at t = t_R for µ > 0. For layer-edge radiative heat fluxes, we may use expression (46) at t_j, j = 0 : R, to get approximations for the downward (+) and upward (-) partial fluxes

where is the unit step function , as well as for the net fluxes

The analytical component of our two-component method is to reconstruct the approximate solution (44) by solving a system of linear algebraic equations for the scalars a_r,i in the S_N solution (22). Inputs to the system are layer-edge values for supplied by the numerical component. The analytical component is to be applied when the intensity of the radiation field I_N(t,µ) at any depth t and direction µ is sought, for then we make use of results (22) and (46) back to (44). We remark that, with the scalars a_r,i determined, the radiative heat fluxes (48-49) can also be obtained at any depth t through results (22) and (46). We describe next either component.

The numerical component of our two-component method is a numerical method designed for solving the S_N diffusive problem (20-21) with no optical truncation error. That is to say, the numerical solution of the S_N diffusive problem (20-21) generated by our numerical method agrees to the analytical solution (22) on corresponding layer-edge points, apart from computational finite arithmetic considerations and regardless of the optical thicknesses of the layers forming the multislab domain. Our numerical method is an extension to anisotropic scattering of arbitrary (Legendre) order and to depth-dependent anisotropic sources of the former spectral Green's function (SGF) method for neutron transport problems (Barros and Larsen, 1990). For this reason, we refer to our numerical method as the extended spectral Green's function (ESGF) method.

The ESGF method has two main ingredients: one is standard and the other is nonstandard. The standard ingredient is the derivation of the local radiative balance equations

where ,

are layer-edge intensities;

and, , is given by definition (47). The local balance Eqs. (50), together with the discrete boundary conditions (21), form a system of linear algebraic equations whose unknowns are layer-average and layer-edge intensities. Such a system does not have a unique solution because there are more unknowns than equations. A simple count gives us N(R+1) layer-edge intensities, NR layer-average intensities, NR balance equations and N boundary equations. The net score is N(2R+1) unknowns against N(R+1) equations. Therefore, N(2R+1) – N(R+1) = NR additional equations are needed or, equivalently, N additional equations per layer relating layer-edge and layer-average intensities.

The nonstandard ingredient is to provide NR equations to the system of N(R+1) equations referred to in the preceding paragraph. These NR additional equations are the ESGF auxiliary equations

where the layer-dependent coefficients q_r,m,u and g_r,m are determined in a layer-by-layer process in order that the analytical solution (22) do satisfy the ESGF auxiliary Eqs. (52), for arbitrary scalars a_r,i and for the entries of vector (25) given by the exponential functions (35). Since the basis vectors (23) have entries whose character depends on whether the layer is non-conservative (all entries are exponentials) or conservative (two vectors have polynomial entries and the remaining vectors have exponential entries), our presentation will be oriented by layer type. So, let us consider a layer in the multislab domain W and let us calculate the layer-average intensities in terms of analytical results derived in the preceding section. If the target layer is non-conservative, then we consider the S_N analytical solution (22) with the entries of all basis vectors (23) given by (26), we substitute the entries of (22) into definition (47) and, after some Calculus, we obtain

for m = 1 : N. However, if the target layer is conservative, we shall consider that two basis vectors (23) have entries given by the elementary solutions (33) and (34), respectively, and the remaining vectors have entries given by (26). If we tag these two basis vectors with i = 1 and i = 2, and if we proceed as in the non-conservative case, then we can obtain

Since the ESGF auxiliary Eqs. (52) are to be satisfied by the S_N analytical solution (22) for arbitrary scalars a_r,i, for the entries of the particular solution vector (25) given by the exponential functions (35), and for every layer in the multislab domain, we can substitute result (53) and the entries of the analytical solution (22) (with entries (26) for the basis vectors (23)) evaluated at t_r-1 and t_r into the ESGF auxiliary Eqs. (52) to obtain for a non-conservative layer the S_N equations

for m = 1 : N. For a conservative layer, we proceed similarly to get

Since Eqs. (55-56) hold for arbitrary scalars a_r,i , they hold for a_r,i set equal to zero for all i from 1 to N. This yields the coefficients

in either case (conservative and non-conservative). We substitute result (57) into Eqs. (55-56), we cancel equivalent terms out on both sides of the equal sign and we arrive at

for a non-conservative layer, and

for a conservative one, with m = 1 : N. Since Eqs. (58-59) hold for arbitrary scalars a_r,i , they must hold for the N ordered sets of scalars (d_1j, d_2j, d_3j, , d_Nj), j = 1 : N. These ordered sets of scalars yield, for fixed m, the system of N linear algebraic equations

for a non-conservative layer, and

and

for a conservative one. The unknowns in either system are the N coefficients q_r,m,u , m fixed, u = 1 : N.

At this point, we select appropriate depths t_r,j, j = 1 : N. Once we have assumed that the layer-dependent separation constants are pairs of ± real numbers, we set t_r,j = t_r for n_r,j > 0 and t_r,j = t_r-1 for n_r,j < 0. In so doing, the arguments in the exponentials in Eqs. (60) and (61c) are all non-positive real numbers. This makes all exponential evaluations in Eqs. (60) and (61c) fall in the interval [0,1], preventing this way computer overflow exceptions when solving the systems (60) and (61) on a digital computer. With our choice, Eqs. (60) and (61c) become

and

For a fixed layer, we may run over all m and, for each m, obtain a corresponding system of N linear algebraic equations in the N unknowns q_r,m,u, u = 1 : N. If the layer is non-conservative, the appropriate equations are just (62-63). If the layer is otherwise conservative, the appropriate equations are (61a), (61b), (62) and (63). With the coefficients q_r,m,u, m = 1 : N, u = 1 : N, determined, we may go back to expression (57) to calculate the coefficients g_r,m, m = 1 : N. We run over all r, calculating the layer-dependent coefficients q_r,m,u and g_r,m and defining completely the ESGF auxiliary Eqs. (52). Equations (50), (52) and the discrete boundary conditions (21) constitute the ESGF equations of the S_N diffusive problem defined by the S_N Eqs. (20) and the above boundary equations. For the numerical solution of the ESGF equations, we use the one-cell block inversion (CBI) iterative scheme (Barros and Larsen, 1990; de Abreu et al., 1998) and we determine the layer-average and the layer-edge intensities for all r and m with no optical truncation error.

The analytical component of our two-component method is a local (layer-level) analytical reconstruction scheme of the approximate solution (44). It is based upon solving a local system of N linear algebraic equations whose unknowns are the scalars a_r,i , r fixed, i = 1 : N. Inputs to the system are the layer-edge intensities that are incident upon the layer of interest (de Abreu and Barros, 1994). These layer-edge intensities are supplied by the numerical component. Since the numerical solution of the S_N diffusive problem (20-21) generated by the ESGF method is free from optical truncation error, the numerical values for layer-edge intensities generated by the ESGF method are exactly the same as those generated by the S_N analytical solution (22) in closed form on corresponding layer-edge points. Accordingly, the system of N linear algebraic equations

and

in the N unknowns a_r,i , i = 1 : N, must hold for a non-conservative layer. For a conservative one, we have the system of N linear algebraic equations

and

With the system (64-65) or (66-67) solved and the scalars a_r,I determined, we can make use of the analytical results derived in the preceding section and of results (46) back to (44) to reconstruct the angular distribution I_N(t,µ) at any depth t on the layer of interest. We can also obtain the radiative heat fluxes (48-49) at any depth t with results (22) and (46). We remark that our reconstruction scheme is local in the sense that calculation of the scalars a_r,i and later reconstruction of the angular distribution I_N(t,µ) are confined to the layer of interest. We may reconstruct over as many layers as we wish in this layer-by-layer process. It is important to note here that if r = 1, then

and that if r = R, then

in Eqs. (64) through (67).

We next discuss computational aspects of the mathematical method described in this section and we present numerical results for a test problem.

A Test Problem

We have codified our mathematical method using a standard version of the FORTRAN language (Press et al., 1996). Thus, our computer version of the method is a FORTRAN program consisting of routines written by the present author and routines picked up from mathematical software packages available elsewhere. It has been oriented by the two-component method described in the preceding section – given the geometrical, material and boundary datasets specifying the target problem (1-3), and the appropriate Gauss-Legendre quadrature set, we solve firstly the mathematical problems defined in the numerical component and we only then solve the mathematical problems defined in the analytical component. The most relevant problems and corresponding routines in the numerical component are: i) the algebraic eigenvalue problem (28), solved with the driver routine RG from the EISPACK system (Smith et al., 1976) for the separation constants n_r,i and angular components a_r,m(n_r,i); ii) the matrix problems (36) and (37), solved with our own routines PARTMAT (to construct and to store the matrices involved), INVMAT (a routine based on the Gram-Schmidt orthogonalization procedure for matrix inversion), and PARTSOL for the matrix operations to yield the constants f_r,m; iii) the system (62-63) for non-conservative layers or (61a-b),(62) and (63) for conservative ones, solved with the routines DGEFA and DGESL from the LINPACK collection (Dongarra et al., 1979) for the coefficients q_r,m,u; iv) the algebraic problem (57), solved with our own routine COEFGRM for the coefficients g_r,m and v) the ESGF equations, solved with our own routine ESGF for the layer-average and layer-edge intensities. Additional routines have been written in order to generate the diffuse component of the angular distribution of leaving radiation (RLEAVE) from results (45) and (46), as well as heat fluxes (HEATF) from results (46) and (48-49). In the analytical component, we consider the problem of solving the system (64-65) or (66-67) for the scalars a_r,i. For then we used over again the routines DGEFA and DGESL from LINPACK. We have written an additional routine (RECON) to locally reconstruct the diffuse component of the angular distribution of radiation at any t from results (22), (45) and (46). The routine HEATF may also be called to yield heat fluxes at any t.

We find it appropriate at this point to discuss the usefulness of our definitions for multislab conservative and non-conservative problems from a computational perspective. First of all, we acknowledge that a more physically acceptable definition is as follows: a radiative transfer problem defined on a multislab domain W is conservative if v(t) = 1 everywhere in W. It is otherwise non-conservative (0 < v(t) < 1 for some t on W). However, in our opinion, the definitions set in the introductory section allow us to construct more computationally efficient algorithms for the determination of the coefficients q_r,m,u in step iii above, as well as for determining the scalars a_r,I and the diffuse component (45) in the analytical component. According to our definitions, if a multislab problem is non-conservative, then there are no conservative layers in the multislab domain and therefore, we do not have to perform layer-type tests on a layer-by-layer basis to decide which systems are to be solved in step iii and in the analytical component. In our program, we have used a flag to indicate whether the current problem is conservative or not.

We illustrate the numerical accuracy of our mathematical method with numerical results for a test problem relevant to shortwave atmospheric radiation. We should notice that the numerical results reported here come from the execution of our FORTRAN program on an IBM-compatible PC (1.4 GHz-clock Intel Pentium 4 processor and 256 Mbytes of RAM) running on GNU/Linux, version 0.2. The executable file of our program has been generated with the g77 GNU Fortran package, release 2.95. The execution (CPU) time for the test problem was 146.4 seconds. This time was generated with the TIME GNU internal routine, option –S.

Our test problem is a multislab conservative problem based on a six-layer model for a stratified atmosphere described in a work of Devaux et al. (1979). Each of the six layers has the same scattering law but the single scattering albedo is allowed to be different in each layer. The optical thickness Dt_r and single scattering albedo v_r for each layer are provided in Table 1. The scattering law is approximated by the scattering phase function data given in Table 2. The multislab domain is irradiated with a mixed beam having a normally incident monodirectional component and a linearly anisotropic regular component. The boundary data for our test problem are t₀ = 0, t₆ = 21, I₀ = 0.5, µ₀ = 1 and g₀(µ) = µ, µ > 0.

In Tables 3 and 4, we present layer-edge results for the diffuse component of our approximate solution I₂₀₀(t,µ) and for the radiative heat fluxes (48-49) for N = 200, respectively. We notice that the directions ±µ_m, m Î C º {14,27,42,60} in Table 3 are discrete directions in the S₂₀₀ Gauss-Legendre quadrature set. Because of the form of the boundary function (g₀(µ) = µ) and of the equivalence between our S_N formulation (20-21) and the P_N-1 formulation (15) with Mark boundary conditions, the numerical values of µ_m, m Î C, are exactly those listed in the lower half of the t₀ column in Table 3. The sources of discrepancies between the other positive directions (those multiples of tenth) in the µ column in Table 3 and the corresponding values in the lower half of the t₀ column are: i) the intensity I^d(t,µ) is likely to be angularly discontinuous at (t₀,0) and ii) the P_N-1 component (45) is a polynomial fit of I^d(t,µ) in the full-range µ-interval [-1,1]. The numerical values for , j = 0 : 6, m Î C, in Table 3 and the values for radiative heat fluxes in Table 4 agree to within ±1 in the last figures given with corresponding S₂₀₀ results generated on a very fine mesh grid (40 cells per unit optical length per layer) with an up-to-date version of a fine-mesh S_N program developed by us some years ago (de Abreu and Barros, 1994) to check our coarse-mesh S_N analytical and numerical methods for accuracy.

In the next section, we share our major thoughts concerning the work on the mathematical method described in this article and we report ongoing research.

Discussion

We have described a mathematical method for approximately solving multislab radiative transfer problems with anisotropic scattering of arbitrary (Legendre) order. We have considered that the multislab domain is a union of conservative and/or non-conservative layers of homogeneous material each. We assumed further that the multislab domain is irradiated from the left side with a beam of radiation having a monodirectional (singular) component and/or a continuously distributed (regular) component in angle. The method described here is thus suited to radiative transfer problems characterized by different layer types and with a left boundary condition expressed as a sum of (rather) arbitrary singular and regular functions of the angular variable.

From a conceptual viewpoint, the mathematical method described in this article is a threefold method in the following sense: i) the multislab target problem (1-3) is decomposed into two basic radiative transfer problems – the uncollided problem (4-5) with left singular boundary conditions and the diffusive problem (6-7) with regular boundary conditions; ii) the uncollided problem is considered as is and it is solved rather easily and straightforwardly, whilst the diffusive problem is considered approximately through a standard S_N formulation, and the approximate S_N diffusive problem is solved exactly with our two-component method and iii) solutions to the two problems are composed through basic relations in order to yield an approximate solution to the multislab target problem. From a practical standpoint, the mathematical method described here is a constructive method designed to conform to the type of incident beam and to the quantities we are looking at. On one hand, our method is able to yield approximate solutions to radiative transfer problems defined on a multislab domain irradiated from one side with a beam having a monodirectional component and/or a continuously distributed component in angle.

On another hand, our method can be used for generating spatially localized quantities, such as the angular distribution of radiation leaving the multislab domain, without having to worry about the detailed distribution of radiation in space. This is accomplished by the numerical component of our mathematical method. When the detailed space-angle distribution of radiation on a specific layer is required, the layer-edge values generated by the numerical component for the layer of interest are the inputs to a local reconstruction scheme of the angular distribution of radiation at any position on the layer of interest. This is accomplished by the analytical component of our mathematical method.

We are considering an extension of the mathematical method described in this article to radiative transfer problems defined on a multislab domain with interacting boundaries (Chandrasekhar, 1950; Thomas and Stamnes, 1999; Liou, 2002). This would enlarge the applicability of our method to problems of theoretical and practical interest not covered thus far, e.g. radiative transfer problems with specularly and/or diffusely reflective boundary conditions. Target problems with reflecting boundaries are in general more complicated to formulate and solve than those with transparent boundaries, for the implicit boundary conditions induced. To illustrate the complications introduced by reflecting boundaries into multislab problems with mixed boundary conditions, let us consider a radiative transfer problem defined on a multislab domain W with reflecting boundaries. Let us assume that the boundaries are specularly and diffusely reflective. For then the boundary conditions (3) are replaced by

and

where r_s,0, r_d,0, r_s,R and r_d,R denote specular (s) and diffuse (d) reflectances on the left (0) and right (R) boundaries. Equations (1), (2) and (68) define the target problem representing a class of radiative transfer problems with specularly and diffusely reflecting boundaries. In order to construct separate singular and regular transfer problems, we Chandrasekhar decompose the target problem (1), (2), (68) into an uncollided problem defined by Eq. (4) with the implicit boundary conditions

and a diffusive problem defined by Eq. (6) with the implicit boundary conditions

so that

We shall give a comprehensive account of these matters in a forthcoming article.

Paper accepted July, 2005.

technical Editor: Atila P. Silva Freire.

• Barros, R.C. and Larsen, E.W., 1990, "A Numerical Method for One-Group Slab-Geometry Discrete Ordinates Problems with No Spatial Truncation Error", Nuclear Science & Engineering, Vol. 104, pp. 199-208.
• Barros, R.C. and Larsen, E.W., 1991, "A Numerical Method for Multigroup Slab-Geometry Discrete Ordinates Problems with No Spatial Truncation Error", Transport Theory & Statistical Physics, Vol. 20, pp. 441-462.
• Benassi, M., Garcia, R.D.M., Karp, A.H. and Siewert, C.E., 1984, "A High-Order Spherical Harmonics Solution to the Standard Problem in Radiative Transfer", The Astrophysical Journal, Vol. 280, pp. 853-864.
• Case, K.M. and Zweifel, P.F., 1967, "Linear Transport Theory", Addison-Wesley Publishing Co., Reading MA, 342 p.
• Chandrasekhar, S., 1950, "Radiative Transfer", Oxford University Press, London, 393 p.
• Davison, B., 1957, "Neutron Transport Theory", Clarendon Press, London, 450 p.
• de Abreu, M.P., Alves Filho, H. and Barros, R.C., 1993, "An accelerated numerical method for one-speed one-dimensional eigenvalue problems in neutron transport theory with no spatial truncation error", Bulletin of the Brazilian Society of Computational and Applied Mathematics, Vol. 4, No. 1, pp. 1-9.
• de Abreu, M.P. and Barros, R.C., 1994, "An Analytical Reconstruction Scheme for the Dominant Solution of One-Speed Slab-Geometry S_N Eigenvalue Problems", Proceedings of the 5^th Brazilian Nuclear Society Meeting, Vol. 1, Rio de Janeiro, Brazil, pp. 201-205.
• de Abreu, M.P., 1995, "Implicit Transmission Conditions for Discrete Ordinates Eigenvalue Problems in Neutron Transport Theory" (In Portuguese), Proceedings of the 10^th National Meeting on Reactor Physics and Thermal-Hydraulics, Vol. 1, Águas de Lindóia, Brazil, pp. 516-521.
• de Abreu, M.P., Alves Filho, H. and Barros, R.C., 1996, "A Numerical Method for Multigroup Slab-Geometry Eigenvalue Problems in Transport Theory with No Spatial Truncation Error", Transport Theory & Statistical Physics, Vol. 25, pp. 61-85.
• de Abreu, M.P., 1997, "One-Speed Slab-Geometry Discrete Ordinates Neutron Transport Problems with Anisotropic Scattering of High Order", Proceedings of the 11^th National Meeting on Reactor Physics and Thermal-Hydraulics (CD-ROM), Águas de Lindóia, Brazil.
• de Abreu, M.P., Barros, R.C., Yavuz, M., Alves Filho, H. and Mello, J.A.M., 1998, "Progress in Spectral Nodal Methods Applied to Discrete Ordinates Transport Problems", Progress in Nuclear Energy, Vol. 33, pp. 117-154.
• de Abreu, M.P., 1998, "On the Spectrum of the One-Speed Slab-Geometry Discrete Ordinates Operator in Neutron Transport Theory", Annals of Nuclear Energy, Vol. 25, pp. 1209-1219.
• de Abreu, M.P., 2001a, "A Spectral Nodal Method Compatible with Iteration on the Scattering Source", Proceedings of the 5^th Meeting on Computational Modelling, Vol. 1, Nova Friburgo, Brazil, pp. 120-129.
• de Abreu, M.P., 2001b, "A Numerical Nodal Method for One-Speed Slab-Geometry Discrete Ordinates Equations in an Integral Form", Proceedings of the 5^th Meeting on Computational Modelling, Vol. 1, Nova Friburgo, Brazil, pp. 130-139.
• de Abreu, M.P., 2002, "Numerical Methods for the Generation of the Spectrum of the Multigroup Slab-Geometry Discrete Ordinates Operator in Neutron Transport Theory", Annals of Nuclear Energy, Vol. 29, pp. 1837-1853.
• de Abreu, M.P., 2004, "A Two-Component Method for Solving Multislab Problems in Radiative Transfer", Journal of Quantitative Spectroscopy & Radiative Transfer, Vol. 85, pp. 311-336.
• Devaux, C., Grandjean, P., Ishiguro, Y. and Siewert, C.E., 1979, "On Multi-Region Problems in Radiative Transfer", Astrophysics & Space Science, Vol. 62, pp. 225-.233.
• Dongarra, J.J., Bunch, J.R., Moler, C.B. and Stewart, G.W., 1979, "LINPACK User's Guide", Society for Industrial and Applied Mathematics (SIAM), Philadelphia PA, USA, 320 p.
• Duderstadt, J.J. and Hamilton, L.J., 1976, "Nuclear Reactor Analysis", John Wiley & Sons, New York, 672 p.
• Duderstadt, J.J. and Martin, W.R., 1979, "Transport Theory", John Wiley & Sons, New York, 613 p.
• Ganapol, B.D., Johnson, L.F., Hammer, P.D., Hlavka, C.A. and Peterson, D.L., 1998, "LEAFMOD: A New Within-Leaf Radiative Transfer Model", Remote Sensing of Environment, Vol. 63, pp. 182-193.
• Ganapol, B.D., Johnson, L.F., Hlavka, C.A., Peterson, D.L. and Bond, B., 1999, "LCM2: A Coupled Leaf/Canopy Radiative Transfer Model", Remote Sensing of Environment, Vol. 70, pp. 153-166.
• Garcia, R.D.M. and Siewert, C.E., 1985, "Benchmark Results in Radiative Transfer", Transport Theory & Statistical Physics, Vol. 14, pp. 437-483.
• Lenoble, J. (editor), 1985, "Radiative Transfer in Scattering and Absorbing Atmospheres: Standard Computational Procedures", Deepak Publishing, Hampton VA, USA, 300 p.
• Lewis, E.E. and Miller Jr., W.F., 1993, "Computational Methods of Neutron Transport", American Nuclear Society, La Grange Park IL, USA, 401 p.
• Liou, K.N., 2002, "An Introduction to Atmospheric Radiation", 2^nd edition, Academic Press, New York, 597 p.
• Myneni, R.B. and Ross, J. (editors), 1991, "Photon-Vegetation Interactions: Applications in Optical Remote Sensing and Plant Ecology", Springer-Verlag, New York, 565 p.
• Press, W.H., Teukolsky, S.A., Vetterling, W.T. and Flannery, B.P., 1996, "Numerical Recipes in FORTRAN 77: the Art of Scientific Computing", Cambridge University Press, New York, 933 p.
• Russell, B., 1976, "The Problems of Philosophy", Oxford University Press, London, 225 p.
• Schneider, H. and Barker, G.P., 1989, "Matrices and Linear Algebra", Dover Publications, New York, 413 p.
• Shultis, J.K. and Hill, T.R., 1976, "The Discrete Eigenvalue Problem for Azimuthally Dependent Transport Theory", Nuclear Science and Engineering, Vol. 59, pp. 53-56.
• Shultis, J.K. and Faw, R.E., 1996, "Radiation Shielding", Prentice-Hall, Englewood Cliffs NJ, USA, 533 p.
• Siewert, C.E., 1978, "The F_N Method for Solving Radiative Transfer Problems in Plane Geometry", Astrophysics & Space Science, Vol. 58, pp. 131-137.
• Siewert, C.E., 2000, "A Concise and Accurate Solution to Chandrasekhar's Basic Problem in Radiative Transfer", Journal of Quantitative Spectroscopy & Radiative Transfer, Vol. 64, pp. 109-130.
• Smith, B.T., Boyle, J.M., Dongarra, J.J., Garbow, B.S., Ikebe, Y., Klema, V.C. and Moler, C.B., 1976, "Matrix Eigensystem Routines EISPACK Guide", Springer-Verlag, Berlin, 551 p.
• Stamnes, K., Tsay, S., Wiscombe, W. and Jayaweera, K., 1988, "Numerically Stable Algorithm for Discrete-Ordinates-Method Radiative Transfer in Multiple Scattering and Emitting Layered Media", Applied Optics, Vol. 27, pp. 2502-2509.
• Thomas, G.E. and Stamnes, K., 1999, "Radiative Transfer in the Atmosphere and Ocean", Cambridge University Press, New York, 517 p.

Publication Dates

• Publication in this collection
  02 Jan 2006
• Date of issue
  Dec 2005

History

• Accepted
  July 2005
• Received
  July 2005