Three-dimensional radiative transfer with polarization in a multiple scattering medium exposed to spatially varying radiation

J. Quanr. Sj~ccrro~c. Radiaf.

Pergamon PII: SOO224073(%)0008~1

Transfer Vol. 57, No. 1, pp. 81-105, 1997 Coovriaht 0 1997 Elsevier Science Ltd Printed’k &eat Britain.All rights reserved 0022473/97 517.00 + 0.00

THREE-DIMENSIONAL RADIATIVE TRANSFER WITH POLARIZATION IN A MULTIPLE SCATTERING MEDIUM EXPOSED TO SPATIALLY VARYING RADIATION D. W. MUELLER, Thermal

JR. and A. L. CROSBIEt

Radiative Transfer Group, Department of Mechanical and Aerospace Engineering Engineering Mechanics, University of Missouri-Rolla, Rolla, MO 65401. U.S.A.

and

(Received 5 February 1996)

Abstract-The present investigation of three-dimensional radiative transfer accounts for polarization and multiple scattering by using the vector transport equation. The intensity vector is comprised of the four Stokes parameters, which are functions of three coordinates describing position and two angles specifying direction. Scattering is characterized by a general phase matrix for randomly oriented particles with a plane of symmetry. The medium is finite in the z-direction and infinite in the X- and y-directions. The incident radiation is polarized and spatially varying on the upper boundary. The vector transport equation is converted, using a double Fourier transform, to an ordinary integro-differential equation with a complex coefficient. The problem of collimated radiation incident on the upper boundary of a cylindrical medium is shown to be of fundamental importance. The analysis is greatly simplified by expressing the transformed intensity vector in a form that depends on only one transform variable. Using superposition, a linear Fredholm integral equation of the second kind is developed for the 4 x 4 generalized source matrix. The bidirectional reflectance and transmittance of a scalar analysis are replaced by 4 x 4 reflection and transmission matrices. The emergent radiation can be determined from these matrices for incident radiation of any polarization. Symmetry relationships are developed using successive orders of scattering and symmetry properties of the phase matrix. Copyright 0 1997 Elsevier Science Ltd

INTRODUCTION

The propagation of a laser through a multiple scattering medium is important to researchers in many fields including material processing, remote sensing, particle sizing, optical communications, and biomedical and ocean optics. Often the beam is polarized by the internal optics of the laser, and the finite size of the beam requires multidimensional analysis. However, most radiative transfer researchers neglect polarization because of the mathematical complexities associated with solving the vector transport equation. In multidimensional situations, the analysis becomes much more complicated, and polarization is almost always neglected. The purpose of the present investigation is to include the effects of polarization in the theoretical analysis of three-dimensional radiative transfer. In one-dimensional situations, the proper treatment of polarization in a multiple scattering medium is well documented.‘4 The radiative transfer is described by a vector transport equation for the Stokes parameters. The most common problem considered is the situation of unpolarized, collimated radiation incident uniformly over the upper surface of a homogeneous planar medium;7-‘3 although several investigations have considered linearly and circularly polarized incident radiation.14.15 Rayleigh scattering is usually assumed, and the effects of neglecting polarization in this situation have been examined.‘~‘6’8 Other researchers have compared the scalar (neglecting polarization) and vector approaches for more complicated scattering laws.‘9.‘oRecent tTo whom all correspondence

should be addressed. 81

82

D. W. Mueller, Jr. and A. L. Crosbie

studies have focused on scattering by non-spherical particles,2’-24surface reflection,2s‘27the inverse problem,**.29and the internal intensity so1ution.30.3’ Few investigations of multidimensional radiative transfer with polarization exist in the open literature. Several experimental studies32A2have examined the problem of a polarized laser beam directed normal to the upper surface of a planar, scattering medium. These experiments reveal that a linearly polarized beam can make an apparent two-dimensional problem into a three-dimensional problem. Only two theoretical investigations of two-dimensional radiative transfer have been found that include polarization. In both studies, the problem was two-dimensional because of variations in the optical properties. Jin43 formulated the problem of two-dimensional radiative transfer in an inhomogeneous, scattering and emitting planar medium. The zeroth- and first-order solutions for the brightness temperature were computed. Whitt and UlabyM investigated first-order vector radiative transfer in periodic vegetation canopies exposed to unpolarized, collimated radiation. No analytical studies of three-dimensional radiative transfer with polarization were found. An integral transform or a separation of variables technique is often applied in a multidimensional problem to reduce the governing equation to a one-dimensional form, that is more suitable for numerical calculation. After solving the one-dimensional problem in the transform space, an inverse transform is performed to construct the multidimensional solution. This technique has been applied to problems that are multidimensional because of spatially varying incident radiation,4s-65 and to those that are multidimensional because of spatially varying optical properties.6673 The classical searchlight problem of a narrow beam incident on an isotropically scattering medium has been the subject of several theoretical investigations.4s-” Crosbie and co-workers have formulated the integral equations in cylindrical coordinates describing two-dimensional anisotropic scatteringSS and three-dimensional isotropic scattering.s6 The three-dimensional rectangular formulations7 has been related to the two-dimensional cylindrical formulation.55 A recent theoretical study considered the problem of three-dimensional radiative transfer in a medium with a refractive index greater than unity.” Numerical results have been obtained for an axisymmetric Gaussian laser beam incident normally on a medium with isotropic,SM’ 1inear,62.63 Rayleigh,@ and other types of anisotropic scattering.65 In the atmospheric sciences, researchers6”73 have applied a finite transform to solve problems that are multidimensional due to spatial variations in the optical properties. None of these multidimensional theoretical investigations included polarization. This paper formulates the vector equations governing three-dimensional radiative transfer in a multiple scattering medium with the effects of polarization. Spatially varying radiation is incident on the top surface of a finite planar medium. The scattering is characterized by a general phase matrix for randomly oriented particles with a plane of symmetry. The method of analysis closely follows the three-dimensional study of Crosbie and Shieh. 58A column vector for the four Stokes parameters replaces the scalar intensity, and the 4 x 4 phase matrix replaces the phase function. Using a double Fourier transform, the vector transport equation is reduced to an ordinary integro-differential equation with a complex coefficient. By considering collimated incident radiation and rewriting the equations in cylindrical coordinates, the transformed intensity vector is related to a vector that depends on only one transform variable. The emergent radiation is described by generalized 4 x 4 reflection and transmission matrices. Symmetry relationships for the generalized matrices are developed by applying a procedure similar to that of Hovenier74 which involves the integral equation for the source matrix and symmetry properties of the phase matrix. FORMULATION

Physical model

This investigation considers three-dimensional radiative transfer with polarization. The radiation incident on the top surface is spatially varying and polarized. The medium is assumed to be infinite in the X- and y-directions and finite in positive z-direction. The medium is absorbing, anisotropically scattering, and homogeneous. The index of refraction is unity. The extinction coefficient, #L, the sum of the absorption coefficient and the scattering coefficient, and the single-scattering albedo, w, the ratio of the scattering coefficient and the extinction coefficient,

Three-dimensional

radiative

transfer

83

characterize the medium. The scattering is specified by the scattering matrix which is described later. Location within the medium is specified by three optical coordinates, z, = flex, zJ = &J, and z; = /?J, which are the product of the extinction coefficient and the corresponding physical coordinate. The optical thickness of the medium is finite and given by z. = j&L. In many situations, the cylindrical coordinate system is more convenient, in which case t, and T,.are replaced by z, and +, where zr = /$r = dm and $ = arctan(r,./r,). The global azimuthal angle 1(1is measured clockwise from the z,-axis when looking in the positive z,-direction. Figure 1 shows the rectangular and cylindrical coordinate systems. The direction at a specific location is specified by two angles 8 and C$which are defined with respect to a local coordinate system (see Fig. 1). The axes of the local coordinate system are denoted by r,,, rY,, and Q. The polar angle 8 is measured from the positive t,-axis. For notational simplification, let p = cos 8. The local azimuthal angle, 4, is measured clockwise from the r,,-axis when looking in the positive t,,-direction. In polarization studies, quantities are defined with respect to a plane of reference. The meridional plane associated with a certain direction is defined as the

top view

meridional plane ’ Fig. 1. The coordinate

system.

84

D. W. Mueller, Jr. and A. L. Crosbie

T

direction of propagation

Fig. 2. Polarization ellipse.

plane that passes through that direction and the r:.-axis. The relationship coordinate systems and angular variables is shown in Fig. 1.

between the various

Stokes parameters

In this investigation, the four Stokes parameters I, Q, U, and V are used to completely describe the radiation. The intensity of the radiation is represented by I, while the state of polarization is characterized by Q, U, and V. In radiative transfer studies, the Stokes parameters are often considered to be elements of a 4 x 1 column vector I

known as the intensity vector. The superscript T denotes the transpose. Each of the Stokes parameters is a function of both position and direction. A geometric representation is often useful to describe the state of polarization. Consider a polarized beam of radiation whose electric vector traces out an ellipse in the plane perpendicular to direction of propagation as shown in Fig. 2. Propagation is in the direction r x 1. I is the sum of the intensity in the I and f-directions; Q is the difference between the intensity in the l- and i-directions. U is the difference between the intensity components in the 45” and 135” directions measured from the I-axis. The angle that the major axis of the ellipse makes with I-axis determines the sign of U. V is related to the circular polarization. The direction that the electric vector traces the ellipse determines the sign of V. The polarization is right-handed if the electric vector moves in the clockwise sense, when looking in the direction of propagation, and left-handed if the electric vector moves in the counterclockwise sense. An important point to emphasize is that the Stokes parameters are defined with respect to a plane of reference that passes through the direction of propagation and the I-axis. For a specific beam, the direction of propagation is known, but the I-axis is chosen for mathematical convenience. The transformation between the Stokes parameters in one reference plane to the Stokes parameters in another reference plane that is rotated counterclockwise through an arbitrary angle 5, when looking in the direction of propagation, is accomplished by the matrix multiplication I,,, = L(5)L

(21

Three-dimensional

where for this choice of parameters,

L(5) =

85

radiative transfer

L(r) is simply

[;

::izt

:;:

(3)

;I.

A more detailed description of the Stokes parameters is available in Refs. 1-6. Scattering matrix

A single scattering event is described by a 4 x 4 matrix that transforms the Stokes parameters of an incident beam to those of a scattered beam, i.e. 1, = F(0)&“,.

(4)

This matrix multiplication represents a change in the direction of the radiation. For both the incident and scattered Stokes vector in Eq. (4) the plane of reference is the scattering plane, defined by the incident and scattered directions. In this study, only single-scattering matrices of the form

F(B) =

al(@)

b,(0)

0

0

by o

a2(@) o

0

0

a,(O)

b>(O)

(5)

0 0 --b*(0) a.+(@) 1 I are considered. This form of the matrix describes scattering by randomly oriented particles with a plane of symmetry.74.75Each element depends only on 0, the angle between the incoming and scattered radiation, which can be related to the polar and azimuthal angles through cos 0 = #Up’+ (1 - /2)“2(1 - @I)‘,‘*cos(f$’ - C$)

(6)

where (p’, 4’) and (p, 4) denote the incoming and scattered directions, respectively. Two important properties of the single scattering matrix are reciprocity and mirror symmetry.74~7s To illustrate these concpets consider a beam described by I,, = [I,,,, QInc,U,,,, I’inclTtraveling in the direction (p’, 4’) scattered by a single particle through the angle 0 producing the beam I, = [L, QXa, Usca,VW]’ which is traveling in the direction (CL,4). The reciprocity relation can be illustrated by simply reversing time which reverses the scattering event so the incident and scattered directions are switched. Because the direction of propagation is reversed, the angle of the major axis of the ellipse changes, and the sign of the Stokes parameter U must be changed. The mathematical statement of the reciprocity of the scattering matrix is F(O) = A&O)A,

(7)

where the tilde denotes the transpose of the matrix and A, = diag[l, 1, - 1, 11.

(8)

Pre- and post-multiplication by Aj reverses the sign of U. Mirror symmetry dictates that if the original scattering event described above occurs, then the incident beam I,,, = [I,,,, Qinc, -U,.,, - Vi,]’ will produce the scattered beam I,, = [I,,, Qscd,-UK,, - V,,]‘. Both the orientation of the major axis of the polarization ellipse and the direction of the trace on the ellipse are changed. The mathematical statement of mirror symmetry with respect to the scattering plane is F(O) = A,A,F(O)AjA.,,

(9)

where A4 = diag[l, 1, 1, - 11.

(10)

AX4= A3A4= A4A, = diag[l, 1, - 1, - 11

(11)

To simplify the notation, let

86

D. W. Mueller,

Jr. and A. L. Crosbie

which has the property that Ax4AJ4= diag[l, 1, 1, 11, the identity matrix. Mathematically, the sign inversion of U and V is accomplished by pre- and post-multiplication by A34. The single scattering matrix of the form given in Eq. (5) clearly satisfies the relationship given in Eqs. (7) and (9). Understanding the symmetry properties of the single scattering matrix is necessary to understand the symmetry relationships developed in this paper. Phase matrix

In multiple scattering studies, the convenient reference plane for the Stokes vector is not the scattering plane as in single scattering studies, but the meridional plane at a given point. Thus, in order to properly handle scattering in multiple scattering situations, it is necessary to rotate the reference plane of the incident Stokes vector to the scattering plane, perform the scattering event, and then rotate the reference plane to the emergent direction. Mathematically, these rotations are expressed by a series of matrix multiplications,‘” viz. T(PC,4) = L(n - y)F(W(-y’)I($,

C),

(12)

where (CL’,4’) denotes the incident direction and (p, C#J) denotes the scattered direction. The rotation angles are shown in Fig. 3. From spherical trigonometry, the rotation angles y and y’ can be related to the incident and scattered directions by cos y’ = b(l - P’~)“~- ~‘(1 - ,D’)“’CO~(C#J’ - 4)]/sin 0,

(13a)

sin y’ = (1 - p2)1’2sin($’ - $)/sin 0,

(13b)

cos y = b’(l - p2)lf2- ~(1 - ,B”)“’ COS(C$’ - 4)]/sin 0,

(13c)

sin y = (1 - P”)“~ sin($’ - $)/sin 0,

(13d)

and

where the limiting cases when sin 0 = 0 must be considered separately. Next by defining the phase matrix as IQ, 4; P’, 6) = L(n - y)F(W(-0,

(14)

Eq. (12) can be expressed as I(& #) = P(P1,& CL’,~‘)I($, 4’).

(15)

Actually, the phase matrix does not depend on C#J and @‘, but on the difference 4’ - 4. This choice

Fig. 3. The rotation

angles.

Three-dimensional

87

radiative transfer

of notation is used for consistency with other matrices in the paper. Finally, after performing the matrix multiplications given in Eq. (14), the explicit form of the phase matrix,5.74.76is given by 0

-h(Q)S(Y’)

al(Q)

h(QMY’) h WC(Y) adQkW(y’)

-a2(@k(rMv’)

--a3(@MybW

-a3(@MvMr’)

a2(Ws(y)W

-a2(@Mb(r’)

w/b 45 P’, 4’) =

-~z(@)s(Y)

1 1

W@MY)

+a3(@MMy’)

0

-

(16)

b2(@MY)

+a3(@MyMy’)

bz(O)s(y')

-b*(@kW)

a4(@)

where c(y) = cos 2y and s(y) = sin 2~. Using the properties of the scattering and rotation matrices, the phase matrix can be shown to satisfy useful symmetry relationships.74 Four relationships can be deduced by considering the basic transpositions that leave the scattering matrix unchanged, i.e. interchanging p and p’, interchanging 4 and +‘, simultaneously changing the sign of p and p’, or simultaneously changing the sign of C$and 4’. Each of these transpositions does not change the single scattering matrix, yet obviously affects the rotation matrix. For example, the interchange of p and $ results in interchanging y and y’ which implies that VP, 4; P’, 9’) = A&P’, 4; PC,cP’)& Interchanging

(17)

# and 4’ results in a sign change in Eqs. (13b) and (13d) to give VP, &P’,&)

=

A34W

K;P',

4)&4

(18)

From Eqs. (6) (13a)(13d), interchanging C$and 4’ has the same effect as simultaneously changing the signs of p and p’ so that VP, 4; P’, 4’) = and is mathematically

A34P(-p,

identical to simultaneously P(P, 4; $3 4’) =

A34Wc1,

4;

-P’,

(19)

@)A34

changing the sign of 4 and 4’ so that -4;

P’,

-qW34.

(20)

A fifth relationship can be obtained by combining Eqs. (17) and (18) to yield VP, 4; P’, 4’) =

A4&‘,

6’;

P,

$I&.

(21)

The physical significance of these and other similar relationships are explained in detail by Hovenier.74 Finally, another useful property evident from Eqs. (6) and (13a)-(13d) is VP, 4; CL’,4’) = WC, $J + 5; $3 4’ + 5)

(22)

where again c is an arbitrary angle. Vector transport equation

For the coordinate system shown in Fig. 1, the three-dimensional, equation can be written as

vector radiative transport

where the source vector is given as I

I

S(T.r, T,, 7:; p, 4) = g

P(P, 4; P’, ~‘NG,

ss0

ty, t:; P’, 4’) W W’,

(24)

-I

with emission being neglected. The intensity and the source vectors are functions of both location and direction. The integration in Eq. (24) is over all directions. JQSRT

57/I-D

88

D. W. Mueller, Jr. and A. L. Crosbie

Boundary conditions

In general, the upper boundary is exposed to incident radiation specified as a separable function of position and direction, i.e. I+(%, rY9r: = 0; PY4) =f(z~, while no radiation is incident on the lower boundary,

~j)Ii(P,

(25)

$1,

i.e.

I_(&, q, z, = ro;/i,4) = 0.

(26)

I+ and I- denote the intensity vectors in the positive (0 < p < 1) and negative (- 1 6 ~1< 0) rl-directions, respectively. I&, (6) = [II, Qi, U,, Vi]’ describes the angular variation of the incident radiation. Each element of the vector is a function of direction specified by (p, $), and the spatial variation of the incident radiation f(z,, 7,) is dimensionless. INTEGRAL TRANSFORM

Double Fourier transform

One approach to reduce the problem to a more tractable form is to apply a double Fourier transform.46~s*To apply this transform, multiply each term in Eq. (23) by exp[-i(/?Vr, + &7?)]/2n and integrate over the 7,7,.-plane. i is the positive square root of negative one. /I1 and flVare the variables in the transform space which correspond to 7, and 7,.. The double Fourier transform of Eqs. (23) and (24) can be written as (27)

and

where the transformed

intensity vector is defined as I(z,, 7,., 7; p, ~$1exp[-i(S,7,

+ Br7x)l d7, d7,,

(29)

with 7 replacing zI for notational simplification. The hat symbol indicates that the double integral transform has been applied in the rectangular coordinate system. Applying the transform to the boundary condition at the upper surface yields ~+(B.Yv 8~70; PL,4) =A~Y:. PyPi(P7 9)~

(30)

where (31)

Implicitly assumed throughout the development is that the integral transform of the spatial variation of the incident radiation exists. One commonly considered spatial variation is found in the searchlight problem.4S” In the searchlight problem, the radiation is incident at only one point on the r,rY-plane. The spatial variation of the incident radiation and its transform are given by f(zI, 7?) = 8(7, - 7.&5(7?. - 7?.3 =-.f(/L,

By) = exp[-Q?y7.~,

+ &J/2~

(32)

where the two Dirac-delta functions restrict the loading to an arbitrary point (7, = 7ro, 7,. = 7!,). Cylindrical coordinate system

Although the rectangular representation of the Fourier transform is suggested by the familiar form of the transport equation, in many situations the cylindrical coordinate system is more

Three-dimensional

89

radiative transfer

convenient. For example, when the spatial variation cylindrical coordinates, i.e.

of the incident radiation

f(L r?.) =f(rr cos +, rr sin $) = g(r,, $),

is specified in (33)

the integral transform, as a function of the cylindrical transform variables, can be written as

’g(r,,

cost+ - a)]z, dl// dr,.

$1 exp[-i/k

The bar indicates the double integral transform has been applied in the cylindrical coordinate system. The relationship between the transform variables in the two coordinate systems is given by /? = Jm and a = arctan (/&/fil) which implies that .QL, A) =Y(
(35)

for the searchlight problem can also be expressed in cylindrical

g(t,, II/) = 6(r, - r&W

-

+0)/z,

*

g(B,

Co

=

exp[-i/h,

cosW0

-

cO1/2~

(36)

where (fro, t,bo) designates the arbitrary point where the loading is applied. The transform g given in Eq. (36) can be shown to be identical to the transform p by noting that r,, = trOcos $o, r,., = T,, sin tjo and recalling the relationships between the transform variables. Crosbie and co-workerP5 considered an axisymmetric Gaussian beam centered at the origin. The spatial distribution and transform for this type of loading are given by g(rJ = exp (- rf/rZJ Q g(B) = (l/2H0 exp( - cf,P”/4)

(37)

where T,, is the effective optical radius of the beam. Because the beam is axisymmetric and centered at the origin, the spatial distribution of the radiation does not depend on $ and the resulting transform is not a function of a. Next, because the problem is linear, the spatial variation of the incident radiation can be removed from the boundary conditions by defining f(/L B,,

7; P,

4)

=

f(P

c-4

B

sina,

T;P,

4)

=

S(A

a)f(P,

a, 5

P,

4).

(38)

Thus, radiative transport equation, Eq. (27), can be written as i/Id-

cos(4 - a)]1 = S

(39)

and the source vector, Eq. (28) as

1’ P@,#; $3 #mP , a, r; P’, (6’)W d#‘.

(40)

The boundary conditions become I+@, a, 0; P, 4) = I&

d)

(41)

and I-@, 01,ro; P, 4) = 0.

(42)

The integral transform has converted the partial differential transport equation into an ordinary differential equation with a complex coefficient. The boundary condition depends only on direction, i.e. the spatial variation of the incident radiation has been removed. Note that transformed intensity vector is still a function of five variables: two transform variables, j? and a; two angles, p and 4; and one spatial coordinate, T. If /I = a = 0, Eqs. (39) and (40) reduce to the classical one-dimensional vector transport equation.’ If the phase function replaces the phase matrix and the intensity replaces the intensity vector, Eqs. (39) and (40) reduce to the scalar transport equation in Refs. 46 and 58.

90

D. W. Mueller, Jr. and A. L. Crosbie

Collimated loading This study considers the fundamental problem of collimated radiation incident on the upper boundary of a cylindrical medium, specified as

I+(% $30: P, 4) = gh, tiYo&P - PohvdJ- 40),

(43)

and no radiation incident on the lower boundary. b = [lo, Qo, UO, VOIT describes the polarization of the incident radiation. The angular variation of the collimated radiation is described mathematically by the two Dirac-delta functions. This problem is fundamental because the spatial variation and polarization of the incident radiation are completely general and the solution to problems with other angular distributions of incident radiation can be constructed from the solution for the collimated case. For collimated incident radiation, Eq. (39) can be written as p g + [l + i&/_112

cos($J - LX)]1 = s

(44)

and Eq. (40) as

where the incident of the transformed Next, one of the and Shieh.58 First,

direction of the collimated radiation (b, Cpo)is now included as an argument intensity and source vectors. transform variables, namely a, is eliminated following the procedure of Crosbie replacement of #J by 4 + a, & by +,, + CL,and 4’ by #J’ + a yields p g + [l + i&/Fj?

cos$]T = S

and

x r(P, a, r; P’, 4’ + or;b, $0 + a) d$ d$‘.

(47)

The boundary condition at the upper surface becomes simply

I+(/%a, 0; P, 4 + a; CLO,40 + ~1)= Io6(p - P#(+ - $J~),

(48)

where the spatial variation of the incident radiation has been removed using the procedure described in the previous section. Then, by recalling Eq. (22), the phase matrix can be shown to be independent of a. a can also be eliminated from the limits of integration by noting that both of the functions in the integrand have a period of 2n and integration over any period is equivalent. Thus, Eqs. (46)-(48) no longer explicitly depend on the transform variable a. The only place that tl now appears is as an argument of the intensity and source vectors. Therefore, by introducing two new functions

and Sk;

P, & cco,40; ro) = S(B, a, z; PL,4 + a; PO,40 + a),

(50)

Eqs. (46) and (47) can now be written as p

2 +[l + i&/C-j?

cos $]I1 = SB

(51)

Three-dimensional

91

radiative transfer

and

ss n

I

S&;p,~;Po,$o;To)=~ o _, VP, 4; P’, 4’Yd~;

P’, 4’;

PO, 40;~0)W

W’.

(52)

The transform variable, j?, is included as a subscript, and the optical thickness is explicitly included as an argument. The boundary conditions are I; (0; P, 6; PO,60; To) = IOW - PO)&4 - $0)

(53)

Ii (ro; & 4; /&, $0; ro) = 0.

(54)

and Equations (5 1)-(54) indicate that for the special case of collimated incident radiation the transform variable, cc, can be eliminated, so that the governing equations depend on only one transform variable. The incident radiation can still have an arbitrary spatial distribution. This result can be used to significantly reduce the numerical computations. Inverse integral transform Once the solution in the transform space is known, the complete three-dimensional intensity solution can be constructed for any spatial variation of the incident loading by performing the inverse double Fourier transform, i.e.

In the cylindrical coordinate I(rf,IL,z2;W#Q=z;;

system, the spatial varying intensity vector can be found by o ss

D g(B, a)W, a, r.-; ~1~6) exp[IPr, cosW -

418 da dB

(56)

or for the special case of collimated loading by - a; pi,

~#JO-

a; to) exp[ij?r, cos(+ - a)]fi da d/I.

(57) In general, the intensity vector found from Eq. (57) depends on location (G, +, t:) and direction (A 6). GENERALIZED

SOURCE MATRIX

In this section, an integral equation for the generalized source matrix is developed. This integral equation is rewritten using a linear integral operator, and the source matrix is expressed as infinite series using the successive orders of scattering approach. Four symmetry relations for the generalized source matrix are also developed. Formal solution Using an integrating factor, the formal solution to Eqs. (51)-(54) can be written as Id (r; cc, 4; clo, 40;

~0)

=

IO& - PO)S(~- 40)exp(-az) +

exp[-4z - QlSdr’; ~1,4; PO,$0; s0

and

zo)W/p)

(58)

92

D. W. Mueller, Jr. and A. L. Crosbie

where d = a(~, 4) = [l + i/3JCj2 Substitution

cos $1/p.

(60)

of Eqs. (58) and (59) into Eq. (52) and integration yields

%(r; P, 4; b, &;

70)

=

z

ew(--oo7)P(p,

4;

b,

40)Io

I +to

_,

sr

P(P,

4;

p’,

4’)

‘exp[-0’(7

-

sT

~71

x Sdz’; p’, 4’; PO,40; ~o)(d,u’/- $1 d&’ dz’ n I +z, x

o

ss

Ss(7’;

,u’,

P(P, 4; P’, 4’)

- 01

w[-f-f@

s0 $‘;

PO,

40;

zo)(dp’/p’)

d4’

d7’,

(61)

where b. = a(~~,

40)

=

[l

+

i&/iTZ

cos

4ol/po

(62)

and CT’= o(P’, 4’) = [l + iSJ1_E-z

cos $‘I/$.

(63)

Equation (61) is a rather specific vector relationship that depends on the polarization of the incident loading, IO. This specific relationship can be generalized using superposition. By introducing Jg, the 4 x 4 transformed source matrix associated with scattering medium, i.e. Ss(r; cc, 4; cco,40; TO)= z

Js(z; ~14; PO,40; roYo

(64)

and by combining the integral terms, Eq. (61) can be rewritten as

J&; ~4; PQ,40;

70)

=

exp(--ao7)&,

VP,

4;

4;

P’T

PO,

40)

-

7’,

TO 47

4’)

/A’,

4’)Ja(7’;

p’, 4’; PO, 40;

70)

d$ d$’ dz’,

(65)

s0

with

exp(-~7)/bI

d7, P, 4) =

o i

7lP

>

0

7//l

<

0.

(66)

Equation (65) is a linear Fredholm integral equation of the second kind for the dimensionless B-varying source matrix. Equation (65) involves only the properties of the medium and is more general than Eq. (61) which also depends on the polarization of the incident radiation. An operator notation can be useful in simplifying equations and illustrating mathematical relationships. The operator notation of Busbridge” has been generalized to include polarization74 and an angular varying extinction coefficient. ‘* A similar notation is adopted in this study. By defining the linear operator A,,,,, as &., {.I = t

’ e(z -

X ’ P(PY 4; $9 4’) ss0

z’,

p’,

c#J’){.}

dp’d4’dz’,

(67)

s0

-I

Eq. (65) can be written compactly as Jd7;

P, 4; PQ, 40;

70)

=

ew(-ao7)Wp,

4;

PO,

40)

+

&.+{Jd7’;

~1’~4’; b,

40;

70)>.

(68)

Three-dimensional

93

radiative transfer

The non-homogeneous term is due to the first scattering of the incident radiation. Integral equation (68) shows that the generalized source matrix at a specific location r and direction (cc, 4) depends on the generalized source matrix at all locations, i.e. (0 6 r < ro), and all directions, i.e. (-1
The solution to the linear integral equation given in Eq. (65) can be represented by the following Neumann series

assuming that the series converges. The nth-term in the summation can be related to the previous term by q(r; CL,4; po,40; ro) =

&.m{Ji;-‘CCP’, 4’; PO,40;TO)).

(70)

By representing repeated application of the operator as 4.4

= n,,,n,,,,,~n,“,,,:,.

. . . A,(. I)+(”- ‘l&J- 1))

the nth-term can be related to the zeroth term, i.e. the non-homogeneous

(71) term in Eq. (65) viz.

J;k ~1.4; po,40;zo)= n:,p.6{exp(-a’t’)P(~‘,4’; bLo, 40>>,

(72)

so that upon summation

J&; P, 4;

PO,

40;

70)

=

f

n:,,,{exp(-o’z’)P(~‘,

4’;

PO,

40)).

(73)

“=O

This expression directly relates the source matrix to the phase matrix for which several properties, Eqs. (17H22), are known. Therefore, this notation is useful in the development of symmetry relationships for the source matrix. Symmetry

Four relationships for the generalized source matrix are derived in Appendix A. The approach taken in the derivations is similar to that taken by Hovenier,74 who considered the one-dimensional case. The symmetry properties of the generalized source matrix are based on the symmetry properties of the phase matrix. The relations presented below are used in a later section to prove relationships for the generalized reflection and transmission matrices. The first two relationships are obtained directly from the linear integral equation for the source matrix, Eq. (65). From Eqs. (A3) and (A7) in Appendix A,

Jab; P, 4; PO,40; ~0)= AuJdr; P, - 4; POLO, - 40; TO)&

(74)

and

JB(T;,u,4; PO,#o; to) = &Jf(r; ~1,n - 4; &Lot x - 40;Ed&

(74)

where the asterisk denotes the complex conjugate. These two expressions are a result of the mirror symmetry of the phase matrix with respect to the scattering plane. Next, the operator notation and successive orders of scattering are used to derive two integral relations involving nth-order generalized source matrix. The derivation is based on induction. The general relationship derived in Appendix A are put in a useful form by letting C(T) = exp( - ar) and X(T) = exp( - aor) in Eqs. (AlO) and (A21) which yields

s 0

e“P(-aT)Ji;(T;

0

-Lb

4;

/ho,

$0, To) dr =

s0

‘exP(-CQr)[A&r;

-h,

40; ,% #; %,)A,]dr (76)

94

D. W. Mueller, Jr. and A. L. Crosbie

and

sr”ex~[-~o(~o

=

z)l[&ji(z;

b,

40; p, 4; TO)&] dr.

(77)

0

The specific choice of the functions i(t) and x(r) is made to obtain relations for the generalized reflection and transmission matrices. These two relationships are a result of the reciprocity of the phase matrix. Again, the tilde denotes the transpose of the matrix. GENERALIZED

REFLECTION

AND

TRANSMISSION

MATRICES

In this section, the emergent transformed intensity vectors are expressed as the integral of the source vector for a given state of polarization. More general expressions, for the generalized reflection and transmission matrices, are developed using the generalized source matrix. Finally, using the properties of the generalized source matrix, symmetry relations for the generalized reflection and transmission matrices are obtained. Definition The P-varying reflected intensity vector can be found by evaluating Eq. (59) at r = 0 which yields

I, (0; PL, 6;

Pa,

40; to) =

sr”exp(-4S~(~t; sI0 P,

4; PO,40; To)(dT’Ip).

(78)

0

By defining the generalized b-varying reflection matrix as

RB(P,&PO,@O;TO)= exp(-@J&;

-A

4;

ho,

0

+o;zo)dr,

(79)

Eq. (78) can be written as Ia (0; P, 4; PO,40; 70) = $

BP@, $; PO,40; ro)Io.

(80)

Similary, the #?-varying transmitted intensity vector leaving the lower surface can be expressed by evaluating Eq. (58) at z = ro, i.e. I;(ro;p,

&k,

40;

~~1 =

IO&cl

-

~00)~(4

-

40)exp(-~~0)

Jo

By defining the generalized b-varying transmission matrix as TB(P, $;b,

$O;TO)

=

r"ew[--a(ro

-

z)lJah P, 4; ho,40; TO) dr,

(82)

Eq. (81) can be written as Ia’@&CL9 4; CLO, 40; To) = IO&P - w@(+ - 40) exp( --(Tzo) + $

T,&, 4; &, c$~;ro)Io.

(83)

The first term represents the direct intensity vector passing through the medium, while the second term is the diffuse transmitted intensity vector. The matrices defined in Eqs. (79) and (82) are generalizations of one-dimensional matrices found

Three-dimensional

95

radiative transfer

in Refs. 1,2, and 4 and of the multidimensional scalar reflection and transmission functions in Ref. 58. These matrices are fundamental properties of the medium that depend on the albedo, phase matrix, and optical thickness. The generalized matrices can be thought of as “building blocks” that can be combined in a special way, i.e. the inverse transform, to obtain the spatially varying reflection and transmission matrices for a given distribution of incident radiation. Successive orders of scattering

The successive orders of scattering approach is introduced to relate the source matrix to the phase matrix and to facilitate the derivation of two symmetry relations. In order to use these symmetry relations for the nth-order generalized source matrix, the nth-order generalized source matrix must be related to the nth-order generalized reflection and transmission matrices. By using expressions similar to Eqs. (79) and (82), the nth-order generalized reflection matrix can be expressed as

R;(P, 4;

s” s”

expt-oz)J;tr;

40; ~0)=

PI,

4; PO,do; ~~1dr,

-L

(84)

0

and the nth-order generalized transmission matrix can be expressed as

T;;(P, 4;

PO*$0; 70)=

exp[-+o-

rllJ”8tz;p, 4; po,(bo;to) dr.

(85)

0

As was done for the source matrix, the generalized nth-order matrices can be summed to give the generalized reflection and transmission matrices, i.e.

and

TAP,

9;

PO,

Cpo; ~0)

=

f

T”B(pL,

4;

b6,b;

(87)

70)

II=0

assuming uniform convergence. Symmetry

Based on the symmetry relations developed in Appendix A and presented in Eqs. (74X77), six useful relationships for the generalized reflection and transmission matrices can be obtained. These “local” relationships do not depend on spatial coordinates; however, the six symmetry relations are valid at every point in the transform space. The following expressions are similar to the one-dimensional vector relationships obtained by Hovenier74 and to the multidimensional scalar relationships presented by Crosbie and Shieh.58 The first set of symmetry relationships follow directly from the linear integral relationship between the generalized reflection and transmission matrices and the generalized source matrix. Noting that a(/,~, -4) = a(,~, 4) and using Eq. (74) in Eq. (79) implies that R&T d; bLo,40;70) = AJW,

-9;

PO, -40;

-4;

PO, -90; ro)Aw.

70)&4r

VW

while from using Eq. (82) T&,

6 PO,40; ro) = &T&P,

(88b)

These relationships illustrate the mirror symmetry with respect to the meridional plane passing through the #J = 0 direction and reduce to the one-dimensional results of Hovenier74 for /I = 0. Similarly, recalling that a(~, K - 4) = a*(~, 4) and using Eq. (75) implies that R&L, 4; PO,40; 70)

=

&4RB(pL,

II -

4;

PO,

7~ -

40;

7oU34

WW

96

D. W. Mueller, Jr. and A. L. Crosbie

and T&T 4; &o,60; 70) = &4Tgl(~(rIL- 4; K,, II - $0; r&4.

(89b)

Equations (89a) and (89b) illustrate mirror symmetry with respect to the meridional plane passing through the 4 = n/2 direction. In the one-dimensional problem, the generalized matrices are not complex quantities, and these relationships are identical to those given in Eqs. (88a) and (88b). That these two sets of relationships must be identical in the one-dimensional case is also true for another reason. In the one-dimensional problem, the reflection and transmission matrices can be expressed as functions of 4 - $J,,. In the above relationships 4 - &,, is the same; therefore, the relationships must be the same. However, the fact that these two relationships are different from the two given in Eqs. (88a) and (88b) indicate that the generalized reflection and transmission matrices are actually functions of 4 and $J,,and cannot, in general, be reduced to a function of 4 - &,, as was done to simplify the one-dimensional problem. The scalar equivalents to the above relations were mentioned in Ref. 58. Two additional relationships are developed from the last pair of integral relations derived in Appendix A. Recalling Eq. (73) and applying Eq. (84) implies that the generalized reflection matrix satisfies R;;(pL,4; cco,40; ro) = A&(~Lo, 6; p, 6;

7dA3.

CW

A similar expression for the transmission matrix, i.e. T;@, 4; ~0, 40; ro) = A$&,

6; ~1,9; r&r

(90b)

can be obtained by recalling Eq. (73), and applying Eq. (85). These two relationships are valid for n 2 1. Therefore, the summation can be performed to yield @la)

and T&L, 4; PLO, do;

70)

=

~+t&o>

61;

P,

4;

7oW4.

@lb)

As Hovenier74 explained for the one-dimensional situation, Eq. (91a) is a reciprocity relationship that results from simply reversing the reflection event. Equation (91b) is not a reciprocity relationship but an “exchange” relationship. The physical interpretation of this mathematical relationship is that the incident and transmitted directions are switched, but the radiation is still applied to the top surface. The reciprocity relationship involving the transmission matrix requires the loading to be moved to the bottom surface. The six relationships given in Eqs. (88a)--(89b), (91a), and (91b) are useful to verify analytical expressions, check numerical solutions, or even reduce the number of computations. These equations are also used in the next section to develop symmetry relations for the spatially varying reflection and transmission matrices. SPATIALLY VARYING REFLECTION AND TRANSMISSION MATRICES In this section, the emergent intensity vectors are developed in terms of the spatially varying reflection and transmission matrices. The spatially varying matrices are expressed in terms of the inverse transforms of the generalized matrices described in the previous section. Symmetry relationships for the spatially varying matrices are derived using the properties of the generalized matrices. Definition The reflected intensity vector can be written in terms of a spatially varying reflection matrix, viz.

I-(%, $‘, 0; if, #‘) = where recalling Eq. (57),

&

R(7r,

II/;

PL, Cp; ~30, $0;

70%

(92)

Three-dimensional

97

radiative transfer

x exp[iflz, cos(+ - a)]/3 da d/I. The transmitted

(93)

intensity vector can be given by

f+(G, $3 ro; P, 4) = IO&P - PO)&4 - 90) exp(-r0/p)TO(r,,

yk P. 4; 70) + $

where the spatial variation of the direct component g(p, a) exp[ -ifir

T(r,, +; cc, +; ~0, 40; ro)Io (94)

is

tan 8 cos($ - a)] x exp[iar, cos($ - a)]/? da dfl

(95)

and the diffuse spatially varying transmission matrix is

ss

T(r,,JI;~r,~;~,#o;t,)=~ m~s(P,a)TB(P,#-a;/lo,#o-a;~O) 0

0

x exp[iflr, COS(I,+ - a)]fl da dfl.

(96)

The direct transmitted component is collimated and attenuated according to Beer’s law; the polarization is the same as the incident radiation. The spatial variation of the direct transmitted component can be expressed in terms of the spatial variation of the incident radiation by using the appropriate coordinate shifts presented by Crosbie and FarrelP6 viz. To(z,, $; p, 4; zo) =f(tr cos + + r. tan 13cos 4, r, sin I+G + ~~tan 8 sin 4) =f(r,

+ 5. tan O cos 4, r! + ~~tan

e sin 4).

(97)

Equations (93) and (96) are expressions for the diffuse reflection and transmission from a multidimensional medium exposed to a spatially varying, collimated beam. Once these 4 x 4 spatially varying reflection and transmission matrices are determined, the emergent intensity can be found for any polarization of the incident beam by using Eqs. (92) and (94). Symmetr)

Based on the symmetry relationships for the generalized reflection and transmission matrices and the linear inverse integral transform given in Eqs. (93) and (96) symmetry relations for the spatially varying reflection and transmission matrices can be derived. The derivation of reflection relationships are included in Appendix B. The derivation of transmission relationships are not presented but are similar. The first four relationships are related to the “mirror” symmetry; the second two are related to the “reciprocity”. Mirror symmetry was examined in detail by Hovenier for a plane-parallel medium74 and planetary atmospheres.” The relationships developed in Appendix B are generalizations of the one-dimensional results of Hovenier to the situation of spatially varying incident radiation. In order for the final results to by symmetric, the incident radiation must possess the same symmetry. Thus, if the spatial variation of the incident radiation is symmetric to the t,-axis, then B(r,, ti; P, 4; PO,40; ro) = AwR(r,, -ti;

P, -4;

bL0,-40; ro)&.,.

(98)

while if the spatial variation of the incident radiation is symmetric to the r,.-axis, then R(L 1(1;~1,4; PO,40; 70)

=

AMR(G,7~- $; P, n - 4; /JO, 7~- 40;~o)A34.

(99)

The physical interpretation of these mathematical relationships is similar to the interpretation of the mirror symmetry exhibited by the single scattering matrix. Figure 4 illustrates the concept of spatially varying mirror symmetry for a collimated beam at the origin. Consider a beam

D. W. Mueller,Jr. and A. L. Crosbie

98

/

/zy

mirror symmetry about the

/

2,

-axis

Zy

mirror symmetry about the zy -axis Fig. 4. Spatiallyvaryingmirror symmetry. I, = [lo, Qo, U. , Vo]OIT incident at the origin in the direction (M, 40) which produces a reflected beam I = [I, Q, U, v]’ emergent from the point (7,, +) in the direction (p, 4) as shown by the solid lines in Fig. 4. Mirror symmetry with respect to the T,-axis dictates that the beam Uo, Vo]’ incident at the origin in the direction (b, - 40) produces a reflected beam IO= PO, Qo, I = [I, Q, - U, - v]’ emergent from the point (G, -$) in the direction (cc, -4). Note that the polarization of both the incident and reflected beams, the azimuthal angle of both the incident and reflected beams, and the point of observation are appropriately shifted. The interpretation of mirror symmetry with respect to the r,-axis given by Eq. (99) follows similar arguments. Similar relationships can be obtained for the transmission matrix. If the spatial variation of the incident radiation is symmetric to the r,-axis, then T(7r,

vk

P,

4;

ho,

40;

70)

=

bT(fr,

-$;

CL, -4;

PO> -40;

70)&4,

(100)

while if the spatial variation of the incident radiation is symmetric to the t,-axis, then T(7,,

II/;

p,

4;

pa,

40;

70)

=

&T(tr,

n -

$;

/1, II -

4;

PO, K -

40;

7o)Aw.

(101)

Three-dimensional

radiative

transfer

99

These four relationships, Eqs. (98)(101) are still quite general even though the spatial variation of the incident radiation is restricted. The important case of axisymmetric loading, i.e. the incident radiation does not depend on JI, is permitted. A general discussion of reciprocity is given by van de Hulst4 in Chap. 3 of his book. To illustrate the concept of reciprocity consider an original event in which a “cause” at point 1, produces an “effect” at point 2. Reciprocity dictates that shifting the “cause” to point 2 and reversing the path of the original event will produce the same “effect” as in the original event at point 1. In the context of this study, the “cause” is the radiation at the point (r,,, $,) incident in the direction (b, +O). while the “effect” is the radiation at the point (T,,, &) emergent in the direction (p, 4). Figure 5 illustrates the original event with solid lines and the reverse event with dashed lines. The mathematical statement of the reciprocity of the spatially varying reflection matrix, derived in Appendix B, is given by R(r,,, +z; P, 4; ho, 40; ro) = A&,,,

$1; bco,40 + n; P, 4 + n; zo)A3.

(102)

The azimuthal angle describing incident direction on the left-hand side of Eq. (102) differs from the reflected direction on the right-hand side by rt due to the convention. Implicit in this expression is the fact that for the left-hand side the incident radiation was applied at (tr,, $,) and for the

original event

reverse event Fig. 5. Spatially

varying

reciprocity.

D. W. Mueller, Jr. and A. L. Crosbie

loo

right-hand side the incident radiation was applied at (rcq,tj2). The sign of U is changed because the direction of propagation is reversed. A similar procedure for the spatially varying transmission matrix yields T(r,, &; P, 4; PO,40; ro) = A,+k,,, $1; PO,40 + n; P + TI;~oN3.

(103)

On the left-hand side the incident radiation is applied in the direction (po, 40) at (r,,, Ic/,) on the upper surface, and the emergent intensity is observed in the direction (p, 4) at (rq, r,k) on the lower surface. On the right-hand side the incident radiation is applied in the direction (b, +. + R) at (rr2, &) on the upper surface, and the emergent intensity is observed in the direction (cl, 4 + n) at (TV,,$]) on the lower surface. This result is not a true reciprocity relationship. A true reciprocity relationship for the transmission matrix requires the spatial distribution of the incident radiation to be moved to the bottom surface. CONCLUDING REMARKS The effects of polarization and multiple scattering have been included in the three-dimensional radiative transport equation. The resulting vector transport equation is necessary to properly model realistic problems such as the propagation of a polarized laser beam through a multiple scattering medium. The fundamental problem considered consists of collimated incident radiation with an arbitrary spatial variation on the upper surface of a cylindrical medium that is infinite in the radial direction. The double integral transform method is applied to reduce the three-dimensional equation to a one-dimensional form. Considerable simplification is achieved by showing that the transformed intensity vector can be expressed as a function of only one transform variable. Generalized reflection and transmission matrices are introduced and shown to satisfy several symmetry relationships. The reflection and transmission for any spatial distribution of incident radiation can be found with these generalized matrices. The emergent radiation for any polarization of the incident intensity can also be determined. The present investigation should provide a firm foundation for future studies of polarization in multidimensional situations. Acknowledgemenr-This

work was supported in part by the National Science Foundation through Grant CTS-9103971.

REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26. 27.

S. Chandrasekhar, Radiative Transfer, Dover, New York (1960). J. E. Hansen and L. D. Travis, Space Sci. Rev. 16, 527 (1974). A. Ishimaru, Wave Propagation and Scattering in Random Media, Academic Press, New York (1978). H. C. van de Hulst, Multiple Light Scattering, Academic Press, New York (1980). J. W. Hovenier and C. V. M. van der Mee, Astron. Astrophys. 128, 1 (1983). L. Tsang, J. A. Kong, and R. T. Shin, Theory of Microwave Remote Sensing, Wiley, New York (1985). J. W. Hovenier, Astron. Astrophys. 13, 7 (1971). H. Domke and E. G. Yanovitskij, JQSRT 36, 175 (1986). R. D. M. Garcia and C. E. Siewert, JQSRT 36, 401 (1986). J. F. de Haan, P. B. Bosma, and J. W. Hovenier, Astron. Astrophys. 183, 371 (1987). M. I. Mishchenko, JQSRT 43, 171 (1990). K. F. Evans and G. L. Stephens, JQSRT 46, 413 (1991). W. M. F. Wauben, J. F. de Haan, and J. W. Hovenier, Astron. Astrophys. 282, 277 (1994). A. Ishimaru and R. L. T. Cheung, Ann. Telecommun. 35, 373 (1980). A. Ishimaru and R. L. T. Cheung, Appl. Opt. 21, 3792 (1982). C. N. Adams and G. W. Kattawar, JQSRT 10, 341 (1970). G. W. Kattawar and C. N. Adams, LimnoI. Oceanogr. 34, 1453 (1989). M. I. Mishchenko, A. A. Lacis, and L. D. Travis, JQSRT 51, 491 (1994). J. E. Hansen, J. Atmos. Sci. 28, 1400 (1971). C. N. Adams and G. W. Kattawar, Appl. Opt. 32, 4610 (1993). M. I. Mishchenko, JQSRT 46, 171 (1991). W. M. F. Wauben and J. W. Hovenier, JQSRT 47, 491 (1992). W. M. F. Wauben, J. F. de Haan, and J. W. Hovenier, JQSRT 50, 237 (1993). L. Tsang, J. A. Kong, and R. T. Shin, Radio. Sci. 19, 629 (1984). Q. Ma, A. Ishimaru and Y. Kuga, Radio. Sci. 25, 419 (1990). C. M. Lam and A. Ishimaru, IEEE Trans. Antennas Propugut. 41, 851 (1994). C. M. Lam and A. Ishimaru, IEEE Trans. Antennas Propagat. 42, 145 (1994).

Three-dimensional

28. 29. 30. 31. 32. 33. 34. 35. 36.

37. 38. 39. 40. 41. 42. 43. 44. 45. 46.

47. 48. 49. 50.

51. 52. 53. 54. 55. 56. 57. 58. 59. 60. 61. 62. 63. 64. 65. 66. 67. 68. 69. 70. 71. 72. 73. 74. 75. 76. 77. 78. 79.

101

radiative transfer

C. E. Siewert, JQSRT 30, 523 (1983). N. J. McCormick and R. Sanchez, JQSRT 30, 527 (1983). H. Domke and E. G. Yanovitskij, JQSRT 26, 389 (1981). P. Stammes, J. F. de Haan, and J. W. Hovenier, Astron. Astrophys. 225, 239 (1987). S. R. Pal and A. I. Carswell, Appl. Opt. 24, 3464 (1985). D. C. Look, Jr., Opt. Eng. 28, 160 (1989). D. C. Look, Jr., Proc. SPIE 1166, 518 (1989). D. C. Look, Jr. and Y. R. Chen, Proc. SPZE 1779, 130 (1992). M. Dogariu and T. Asakura, J. Optics (Paris) 24, 271 (1993). C. I(. Carniglia, J. P. Black, S. E. Watkins, and B. J. Pond, Appl. Opt. 32, 5504 (1993). D. C. Look, Jr. and Y. R. Chen, AIAA Paper 93-2727 (1993). D. C. Look, Jr. and Y. R. Chen, J. Thermophys. Heat Trans. 7, 631 (1993). M. Dogariu and T. Asakura, Wave Rand, Media 4, 429 (1994). D. C. Look, Jr. and Y. R. Chen, AIAA Paper 94-2093 (1994). D. C. Look, Jr. and Y. R. Chen, Appl. Opt. 34, 144 (1995). Y. Jin, J. Appl. Phys. 69, 7594 (1991). M. W. Whitt and F. T. Ulaby, Int. J. Rem. Sens. 15, 1813 (1994). H. G. Kaper, J. Math. Phys. 10, 286 (1969). G. B. Rybicki, JQSRT 11, 827 (1971). I. L. Katsev, Izv. Akad. Nauk. BSSR 1, 105 (1973). I. L. Katsev, Zzu. Atmos. Ocean Phys. 10, 425 (1974). M. M. R. Williams, J. Phys. A: Math. Gen. 15, 965 (1982). C. E. Siewert and W. L. Dunn, J. Appl. Math. Phys. 34, 627 (1983). W. L. Dunn and C. E. Siewert, J. Appl. Math. Phys. 36, 581 (1985). C. E. Siewert, JQSRT 33, 551 (1985). C. E. Siewert, JQSRT 41, 467 (1989). B. D. Ganapol, D. E. Kornreich, J. A. Dahl, D. W. Nigg, S. N. Jahshan, and C. A. Wemple, Nucl. Sci. Engng 118, 38 (1994). A. L. Crosbie and R. L. Dougherty, JQSRT 25, 551 (1981). A. L. Crosbie and J. B. Farrell, JQSRT 31, 397 (1984). A. L. Crosbie and L. C. Lee, JQSRT 38, 231 (1987). A. L. Crosbie and S. M. Shieh, JQSRT 44, 299 (1990). A. L. Crosbie and R. L. Dougherty, JQSRT 20, 151 (1978). A. L. Crosbie and R. L. Dougherty, JQSRT 27, 149 (1982). A. L. Crosbie and S. M. Shieh, JQSRT 49, 535 (1993). A. L. Crosbie and R. L. Dougherty, JQSRT 28, 233 (1982). A. L. Crosbie and R. L. Dougherty, JQSRT 33, 487 (1985). A. L. Crosbie and R. L. Dougherty, JQSRT 30, 255 (1983). A. L. Crosbie and L. C. Lee, JQSRT 49, 185 (1993). D. J. Diner and J. V. Martonchik, JQSRT 31, 97 (1984). D. J. Diner and J. V. Martonchik, JQSRT 32, 279 (1984). J. V. Martonchik and D. J. Diner, JQSRT 34, 133 (1985). G. L. Stephens, JQSRT 36, 51 (1986). G. L. Stephens, J. Afmos. Sci. 45, 1818 (1988). G. L. Stephens, J. Atmos. Sci. 45, 1837 (1988). T. Kobayashi, J. Atmos. Sci. 48, 22 (1991). P. M. Gabriel, S. C. Tsay, and G. L. Stephens, J. Afmos. Sci. 50, 3125 (1993). J. W. Hovenier, J. Atmos. Sci. 26, 488 (1969). H. C. van de Hulst, Light Scattering by Small Parricles, Dover, New York (1981). J. W. Hovenier and C. V. M. van der Mee, Astron. Astrophys. l%, 287 (1988). I. Busbridge, The Mathematics of Radiative Transfer, Cambridge Univ. Press (1940). B. D. Ganapol and R. B. Myneni, JQSRT 48, 321 (1992). J. W. Hovenier, Astron. Astrophys. 7, 86 (1970).

APPENDIX A The purpose of this Appendix is to prove four relationships for the generalized source matrix that are used to develop symmetry properties of the generalized reflection and transmission matrices. In the derivations below, note that from Eq. (60) a@,#)

= c(lr, -4)

= “*f&x - 4) = -c(-/.&

9),

Y, 4) =

47,

= e*(7,

-k

(Al)

which from Eq. (66) implies that 0,

cc, -4)

P, n -

4)

= 4-7,

$1,

(A2)

where the asterisk denotes the complex conjugate. The first two derivations use the integral equation for Jb. Eq. (65), while the next two derivations use the operator notation introduced in Eq. (67).

102

D. W. Mueller, Jr. and A. L. Crosbie

(i) Show that Ja(r; ~3 4; no, 40; ro) = AYJc(T; p, -$;

n’o,-40; ro)Aw.

(A3)

To obtain the right-hand side of the above equation replace 4 with -I#J and & with -&, in Eq. (65) which yields J&; n, -&no,

-40; TO)= exp(-oaor)P(n,

-4; KY,-40)

m

-4; $3 4’)

’ e(r - T’, n’. #)Jp(?‘; /I’, 4’; /IX,,-&;

ro) dy’ d4’ dr’.

(A4)

Next, replacement of 4’ with -4’. utilization of the periodicity of the azimuthal integration to rewrite the limits, and preand post-multiplications by Aw gives AwJk;

P. -$;m,

-&;

TO)AW = exp(-oor)A34P(p,

-$J; IIO,+)Aw

” e(r - T’, n’. -K)AwJr(r’; s0

AwP(p, - 4; p’. -+‘)A,

y’, -4’; /.Lo,-60; ro)Au d$ d4’ dr’

(AS)

where the integrand was rewritten with AMPJ~Aw = A~PA~AMJ~A~. Then, application of Eqs. (20) and (A2) yields AjrJ~(r; n, -4;

no, -400; ro)A,, = exp(--uor)P(y,

4; b, $0)

so

e(r - r’, n’, +‘)[AvJe(r’; /I’, - 4’; no, - 40; so)Aw]dn’ d#’ dr’.

w4&P’,4’)

(A6)

0

Since both sides of Eq. (A3) satisfy the same integral equation, Eq. (A3) is true. (ii) Show that J~(T; n,. 4; pro,&; ro) = AwJB*(T;n, x - 4; w, n - 40; TD)AIG.

(A7)

This derivation is similar to that of Eq. (A3). First, replacement of I#Jwith A - 4, $0 with n - & and 4’ with n - 4’ in Eq. (65) yields

= +g

I P(P.l-

IS-* _,

&P',n-4')

'e*(r- T',p'.+‘)J~(T’; $, n - 4’; h, R - &o;TO)dn’ d4’ dr’ J:

(A@

where Eqs. (Al) and (A2) have been applied. Using the perodicity of the azimuthal integration to rewrite the limits, preand post-multiplying by Aw, recalling Eqs. (20) and (22), and taking the complex conjugate of each term yields A~JB(T; n, n - 4; m, R - 40; %)A,, = exp(-aor)P(n,

4; ho, $0)

* I +zT

'e(r- T',p',+‘)[AMJ$(~‘; n’, a - 4’; h, n - 40; TO)AM] d$ d4’ dr’. s0

_IwG~;P',&)

I$

(A9)

Since both sides of Eq. (A7) satisfy the same integral equation, Eq. (A7) is true. (iii) Show that ’ i(r)&.,{&‘,

4’; ,4 &)X(r’)} dr =

L

” ~(r)a,[n:-,,,,,{P(~~, s0

4’; P, +)W)}]‘A, dr,

(Alo)

where the superscript T denotes the transposed matrix. The derivation follows by induction. First, the case of n = 1 is considered, i.e.

I

x

J-7 0

Interchange T with T',let p’ = -$,

s

o&T)&,.,{P(P’,

4’;

M,

l -

I

-I

'e(T

Y-PC, 4; P’. 4’)

T’, p’, &){P(P’, 4’;

/AI,

+o)X(T’)}d/i d$’ dr’.

(A1 1)

recall Eq. (A2), and apply Eqs. (19) and (21) to write

&I)X(T’)}dr =

$

0

s‘X(T) dr

0

n

I

Is

-1

A&d. 4’;

P.

+)A>

-

’ e(T

T', /l’. &){A#(-ha,

60; ~‘9 $‘)A,&‘)} dkd@

dr’.

(At2)

Three-dimensional

103

radiative transfer

where the property A, = &AM is used. Cancellation and rearrangement of terms yields

s

7"i(~)n,-,.,{P(~',~';~a,

b)x(7')}

d7

=

z

'x(7M7 s'0

0

s

r

x A1

I

[i?0

-I

P(-pa.

40; !J’. 4’)

and application of the definition of the n-operator ’ i(r)k.,{P(p’, s0

‘Oe(r s0

-

T', p'.

+‘){P(p’,

4’; /I,

TA, dfi’d+’ dr’.

4)[(7')}

(A13)

I

yields

4’; pa>40)x(7’)} dr =

’ ~(7)A,tn~.-ao.mn{P(~‘,4’; fi< d)T(~‘))l’Al dr.

L414)

s’0

Thus, relationship (AlO) is true for n = 1. Next, Eq. (AlO) is assumed true for n, and the case of n + 1 is considered. The left-hand side of Eq. (AIO) becomes ” i(rM:‘,!,,{P(lc’, I0

’ i(7) dl s0

4’; PO,Wx(r’)} dr = t

I I x A:-,,

IS0 I -I

P(p”, 4”; p’, 4’)

’ e(r” I’0

where for the linear operator &‘,I,$ = A:_,, A, .@.+.Rearrangement

x

T', p',

d’)P(p’.

9’; pob.&)X(7’) dp’ d4’ dr’

,

(Al 5)

40)W d4’ dr’,

(Al@

of terms gives

’ [s

I

~(T)A:-,.~{PW’,4~“;P’, 4’) e(s” - T’, P’, 4’)) dr VP’, 4’;

0

P+

and application of Eq. (AlO) yields

x

0 [J‘

” e(r - T’. /.t’. &)A,[~:_,,~, {P(p”, 4”; /I. +)i(r”}lTA,dr

Interchange r with T’, let p’ = -p’,

Di(r)4’,!,,{P(p’,

4~‘;PO, &)x(f))

x

1

P(p’, 9’; pa, 40) dp’ d4’ dr’,

(Al3

recall Eq. (A2), and apply Eqs. (19) and (21) to write d7 =

SC

1

[S

‘C(T - T’, Jo’, $J’)A,[LI:,<,.{P(~L”,0”; p. $J)<(T”}I~AY A,?( -pa, 40; /.I’, 4’)A, dp’ d4’ dr’.

0

(Al8)

Cancellation and rearrangement of terms yields

s

’ C(r)A:‘,!.,{P(~‘, 4’; w, &)x(7’)} d7 = w 0 x {ny.,,.6{P(g“,

Finally, application of the definition of the n-operator

1

I#J,“; p, c#~)&“)}j d7 ‘As dp’ d+’ ds:

(A19)

gives

0 $0

t(~)A::,f.~{P(p’, 4’; poo,4o)xW}dr =

where for the linear operator A:tj,,, (iv) Show that ’ ((70

-

7)A:,,{W,

’ ~(r)A,[n:~:,,,jP(~‘, $0

d’; P(.4KW)l’A~ dr.

(A201

= Ar,_r,,+.o~..p.,4. Therefore, the relationship in Eq. (AIO) is true for n 3 1.

4'; PO, d~dx(7'))

d7

=

'x(70

-

7)A.[4,,,,,{W,

d'; P,

$KW)}lT~d~.

(A.21)

104

D. W. Mueller, Jr. and A. L. Crosbie

The derivation of Eq. (A21) is similar to the previous derivation. First, the case of n = 1 is considered

r

x

I

%S -I

Let

T' =

TV -

'&a

-

5,

let

T =

7,, -

T',

w4 4; P’, 4’)

s

’ e(7 - 7’3p’, #‘){P(fi’, 4’; ~0, &)x(7’)} dp’ d# d7’. (A22) 0

and apply Eq. (21) to write

r)n,.,,,(P(~',~';~,~o)~(r')}

d7

=

g

"x(70-

r)dr

*

' &%',d~';r,(b)&

'47

-

7', $7

+')

-I

ss0

x {A&-w,

40; $3 &)&[(r’)} d$ d#’ dr’.

(A23)

Cancellation and rearrangement of terms yields

s '~70

-T)&.+

{P(r’,@;bo,

d~oo)x(r')l d7

=

g

"x(70

-

P(-rco, 40; $, 4’)

~)drb

%7

-

7'7 $9

4')

0

X {P(n’, 9’; ~9 4X(7’)}

A4d$ d4’ d7’

(A24)

T and application of the definition of the n-operator

yields

Thus, relationship (A21) is true for n = 1. Next, Eq. (A21) is assumed true for a, and the case n + 1 considered. Let T' = TO - 7 and 7 = 70and follow the procedure of the previous derivation. Thus, Eq. (A21) is true for n 2 1. APPENDIX

T',

leave p’ unchanged,

B

The purpose of this Appendix is to derive three relationships for the spatially varying reflection matrix. The derivations of the relationships for the spatially varying transmission matrix are similar. In each of these derivations, the spatial distribution of the incident radiation must possess a certain symmetry. (i) Show that R(r,, 9; p. 4; /IG 400;70) = First, pre- and post-multiplication yields &.R(sr,

-$;p.

-4;

/b,

-&~;7o)A,.

AwR(rr,

-+;

p.

-4;

by A,, and replacement of 4 by -4,

=

&

g(B,

x AwR&

pa,

-4~0;

$JOby -&

TO&.

$ by -+,

(Bl)

and a by -a in Eq. (93)

-a)

- 4 +

a; pi, - 40 + a;

TO)AW exp[ifi7, cos( - (I + a)@ da d/J

(B2)

where the negative sign has been used to reverse the limits of integration. Next, the integral can be written in a form similar to that in Eq. (93) by using the periodicity of all quantities to change the limits, noting that cosine is an even function, and recalling Eq. (88a). Thus, ” AwR(rr, -$; p, -4;

bb. --do; 7a)Aw= &

z?(B, -a) x Rp@, $J - a; b, 41,- a; 70)exp[i/?7r,co@ - a)]/3 da d/3. (B3)

The integral in the above equation would be identical to the integral in Eq. (93), if g(jJ, -a) = g(B, a), i.e. d is an even function of a. From Eq. (34), 2 will be an even function of a if g is an even function of $ or symmetric to 7,-axis. Therefore, if the spatial variation of the incident radiation is symmetric to the r,-axis, Eq. (BI) is true. (ii) Show that R(~,,II_;c~.&~,Qo; Pre- and post-multiplication (93) yields

TO) =

AwR(r.,x

-

$;p,n

-

+;wo,n

-

+0;7o)Ar.

(W

by Aa, and replacement of 4 by x - 4, QOby n - 40, + by rt - $, and a by n - a in Eq. m

AwR(ro s - +; p, n - 4; /IO,x - $00;TO)AY= $7

ss0

x AMR&

-*

f(B, x - a) - 4 + a; w, - & + a; 70)Awexp[i/37,cos( - $ + a)]fi da d/J

(85)

Three-dimensional

105

radiative transfer

where the negative sign is used to reverse the limits of integration. Next, noting that cosine is an even function, using periodicity to change the limits, and applying Eq. (88a) yields Adt(~r,

I[ - +; p, x - 4; p’o, n - 40; 7o)A14 = &

N,

n-a)

x R&, 4 -

a; h,

40 -

a; 70)

exp[ibs, cos($ -

a)]/3

da da.

(B6)

The integral in the above equation would be identical to the integral in Eq. (93) if &?, n - a) = g(/I, a). This is true if the loading is symmetric to rY-axis or tj = n/2. Therefore, if the spatial variation of the incident radiation is symmetric to the r,-axis, Eq. (B4) is true. (iii) Show that R(7,,, $2; PC,4; no, 40; 70)

=

&ft(rr,,

$1;

po,

40 +

n; p(. 4

+

tr;

(87)

ra)A~.

First, consider a spatial distribution of incident radiation that is azimuthally symmetric with respect to the origin given by g(7,) with the transform g(b). Now, consider the beam incident in the direction (b, 40) centered at the point, (T,,, +I). Using Eq. (93) the spatially varying reflection at (7,>, 1/1?) in the direction (p, 4) for this situation can be expressed as

x

RB(P, $ - a; hco,40 - a;

70) exp[iPrr,

cos(JI:

-

a)lBda dB (W

where the incident distribution is appropriately shifted from the origin to (T,,, +I). Reciprocity dictates that the propagation of the beam must be reversed, i.e. the spatial distribution must be shifted to (7?, $2), the incident direction must be changed to (p, (b + n), and the emergent direction must be changed to (~0, Cpo+ B). For this situation, the spatially varying reflection at (L!. &) in the direction @CO, 40 + n) is given by

NT,,,

$1;

PO, 90

+

n;

p.

4

+

R;

70)

=

&

g(p)

exp[ - iprr, cos($2 - a)]

x Rp(po, 40 + n - a; p, I$ + n - a; 70)exp[ijr,, cos(IJIl Next, replacement of a by

a +

n

a)]B

da

dB.

(B9)

yields

g(p)

exp [iaTr2cos($2 -

t-f)1

x Rp@, 40 - a; p, 4 - a; 70)exp[-ir?7,, COS($I- a)]B da d/I.

(B)O)

where the property of the cosine is used to change the sign in the exponentials. Finally, utilization of periodicity to rewrite the limits. transposition of both sides. pre- and post-multiplication by A,, and application of Eq. (9la) yields

x

RN@, 4 - a; pu. 40 -

a; 70)

Since, the right-hand sides of Eqs. (B8) and (BI I) are identical, Eq. (B7) is true.

exp[-i/IT,, cos($~ - r)lB da da.

(Bl 1)