Direct collocation meshless method for vector radiative transfer in scattering media

Journal of Quantitative Spectroscopy & Radiative Transfer 163 (2015) 50–62

Contents lists available at ScienceDirect

Journal of Quantitative Spectroscopy & Radiative Transfer journal homepage: www.elsevier.com/locate/jqsrt

Direct collocation meshless method for vector radiative transfer in scattering media Xun Ben a, Hong-Liang Yi a,n, Xun-Bo Yin b, He-Ping Tan a,n a b

School of Energy Science and Engineering, Harbin Institute of Technology, 92 West Dazhi Street, Harbin 150001, PR China Department of Mathematics, School of Science, Harbin Institute of Technology, 92 West Dazhi Street, Harbin 150001, PR China

a r t i c l e in f o

abstract

Article history: Received 25 January 2015 Received in revised form 29 April 2015 Accepted 3 May 2015 Available online 9 May 2015

A direct collocation meshless method based on a moving least-squares approximation is presented to solve polarized radiative transfer in scattering media. Contrasted with methods such as the finite volume and finite element methods that rely on mesh structures (e.g. elements, faces and sides), meshless methods utilize an approximation space based only on the scattered nodes, and no predefined nodal connectivity is required. Several classical cases are examined to verify the numerical performance of the method, including polarized radiative transfer in atmospheric aerosols and clouds with phase functions that are highly elongated in the forward direction. Numerical results show that the collocation meshless method is accurate, flexible and effective in solving onedimensional polarized radiative transfer in scattering media. Finally, a two-dimensional case of polarized radiative transfer is investigated and analyzed. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Radiative transfer Polarization Meshless method Direct collocation Multiple scattering

1. Introduction Solar light propagating through the atmosphere, the ocean, and vegetation is often analyzed in the framework of scalar radiative transfer theory. However, due to the electromagnetic nature of light, the effects of polarization should be considered in the treatment of scattering, reflection, and transmission [1,2]. The interaction of light with media, such as scattering by particles and reflection at phase boundaries, will generally alter the polarization of the incident beam [3,4]. Polarization is also decoupled from the frequency of the light, giving extra degrees of freedom. Accurate and efficient polarized radiative transfer calculations are essential in many applications, such as the retrieval of atmospheric and oceanic constituents from remotely sensed satellite observations [5]. For skylight

n

Corresponding authors. Tel.: þ 86 451 86412674. E-mail addresses: [email protected] (H.-L. Yi), [email protected] (H.-P. Tan). http://dx.doi.org/10.1016/j.jqsrt.2015.05.001 0022-4073/& 2015 Elsevier Ltd. All rights reserved.

detection, the additional polarimetric measurements can significantly improve the retrieval of some aerosol parameters, and several space-borne and airborne instruments have been designed to measure the polarization of skylight. This also has potential applications in other disciplines, such as biomedical optical tomography [6]. Many applications in medical diagnostics profit from polarization properties, e.g. noninvasive glucose sensing, tissue anisotropy and concentration measurements. Considering the importance of polarization information, polarized radiative transfer is an active area of research. The standard vector radiative transfer equation (VRTE), which is applied to the propagation of light in a plane-parallel layer, has been utilized for more than half a century [7]. To investigate the polarization characteristics in graded index media, a corresponding VRTE was derived [8]. Recently, for astrophysical applications, a VRTE in an orthogonal curvilinear coordinate system was derived [9]. Many numerical methods have been proposed and developed to investigate polarized radiative transfer, including Monte Carlo (MC) methods [10–14], FN methods [15], discrete ordinate (DO)

X. Ben et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 163 (2015) 50–62

Nomenclature

Z

A a a B d F H Δh I i, j, k

Greek symbols

i1, i2 L M m Ν Nc Nloc p p pj Q Rd r r S S U V Vr w

stiffness matrix coefficient matrix expansion coefficient of Legendre polynomial source matrix diameter of the support domain, m incident light irradiance, W/m2 sr height, m local discrete length, m radiation intensity, W/m2 sr unit vectors of x-axis, y-axis and z-axis directions, respectively rotation angles rotation matrix single scattering Mueller matrix total number of nodes in the computational domain nodal basis function number of the incident directions for collimated radiation sources number of the local nodes in Vr vector of Legendre polynomial order of Legendre polynomial polynomial basis function of jth order linear polarization aligned parallel or perpendicular to the z-axis, W/m2 sr reflection matrix for Lambertian surface spatial coordinate vector distance between ri and r Stokes vector matrix Stokes vector linear polarization aligned 7451 to the z-axis, W/m2 sr circular polarization, W/m2 sr spatial subdomain weight function

methods [16–19], successive order of scattering (SOS) methods [20–23], matrix operator methods [24], etc. To compare computational performance for different numerical methods, Kokhanovsky et al. [25] summarized seven codes of vector radiative transfer, including three techniques based on the DO method, two MC methods, the SOS method, and a modified doubling–adding technique. In that reference, benchmark data for one-dimensional vector atmospheric radiative transfer were presented. Recently, a new class of meshless methods has been applied to solve the scalar radiative transfer equation (SRTE) in scattering media [26,27]. Compared to traditional numerical methods based on equation discretization (such as the finite difference method, the finite volume method and the finite element method), meshless methods have an advantage, since they only use scattered nodes to discretize the domain of the problem. In addition, the meshless method is much faster than ray tracing methods such as MC and zone methods. Due to the characteristics of flexibility and convenience, meshless methods have been developed to solve

δ φ κe κe κs λ μ, η, ξ Θ θ ρ τ Ω, Ω0 Ω0 ω

51

scattering phase matrix

Dirac delta function azimuth angle extinction matrix extinction coefficient, m  1 scattering coefficient, m  1 support domain amplifying factor direction cosines of x-axis, y-axis and z-axis directions, respectively scattering angle zenith angle reflection coefficient optical thickness vector of radiation direction solid angle, sr single scattering albedo

Subscripts 0 c col d i loc row s w

inflow collimated light column diffuse light incident direction local nodes row scattered dirction wall

Superscript T

transposition

radiative transfer problems. To improve the numerical stability of meshless methods, a new second-order form of SRTE was presented by Zhao et al. [28], and Luo et al. [29] applied a kind of upwind scheme to solve the SRTE for strongly inhomogeneous media. To date, the meshless method hasn’t been used to solve VRTE. In this work, a direct collocation meshless method (DCM) is extended to solve the multidimensional VRTE. In this method, the angular discretization is based on the discrete-ordinates approach, where the spatial discretization is conducted using a collocation approach. Here, the trial function is constructed by a moving least squares approximation. Several test cases of polarized radiative transfer are taken to verify the performance of the method. The paper is organized as follows. Firstly, the theory of VRTE is introduced, and Lambertian boundary conditions are presented. Then, the formulation and solution steps of the DCM are presented in detail. In Section 3, two onedimensional cases are studied to test the stability and convergence characteristics of the present method. Finally,

52

X. Ben et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 163 (2015) 50–62

a two-dimensional case for polarized radiative transfer is investigated, and the simulation results are analyzed.

2. Theory 2.1. The vector radiative transfer equation Based on the incoherent addition principle of Stokes parameters, the steady state VRTE for polarized monochromatic radiation in a homogeneous medium without considering thermal emission can be written as Z κs Ω U∇Sðr; ΩÞ ¼ κe Sðr; ΩÞ þ Zðr; Ω0 -ΩÞSðr; Ω0 ÞdΩ0 ; ð1Þ 4π 4π where S¼(I, Q, U, V)0 is the Stokes vector, I is the radiation intensity, Q is the linear polarization aligned parallel or perpendicular to the z-axis, U is the linear polarization aligned 7451 to the z-axis, V is the circular polarization, Ω¼ μiþηj þξk is the unit direction vector of radiation, r is the spatial coordinate vector, κe is the extinction matrix, κs is the scattering coefficient, and Z is the scattering phase matrix describing the scattering properties of the medium. The paths of the incident and scattered rays are shown in Fig. 1, which are comprised of scattering planes and meridian planes. As seen in Fig. 1, OAB is the scattering plane, while OAC and OBC are the meridian planes for the incident ray Si(Ω0 ) and the scattered outgoing ray Ss(Ω), respectively. Since the definition of the polarization component of the incident and scattered rays is not aligned with the scattering plane OAB, a rotation of the Stokes parameters is needed to obtain the scattering phase matrix Z. The transform of the phase matrix Z is expressed as Z ¼ Lð  i2 ÞMðΘÞLð i1 Þ;

ð2Þ

where the matrix M is the single scattering Mueller matrix [7] with 4  4 elements used for describing the transformation of the Stokes vector of the incident wave into that of the scattered wave. For spherical particle scattering, the matrix M can be simplified as block-diagonal with eight non-zero elements, of which only four are independent.

The form is given as 0 M1 M2 0 B 0 B M2 M1 MðΘÞ ¼ B B 0 0 M 3 @ 0 0 M 4

0

1

C 0 C C: M4 C A M3

ð3Þ

where the matrix elements M1, M2, M3 and M4 are dependent on the scattering angle Θ, but independent from the optical depth of a given medium. The matrix L in Eq. (2) represents a rotation in the clockwise direction with respect to an observer looking into the direction of the ray propagation [14], where a negative sign represents the counterclockwise direction. The form is expressed as 0 1 1 0 0 0 B0 cos 2ψ sin 2ψ 0 C B C LðψÞ ¼ B ð4Þ C: @ 0  sin 2ψ cos 2ψ 0 A 0

0

0

1

To complete the transform given in Eqs. (2), i1, i2 and Θ are required. Θ is the scattering angle, and can be obtained by Eq. (5) ! ΩU Ω0  0 : Θ ¼ arccos ð5Þ jΩjΩ  To get the two rotation angles, a vector notation is presented here. As shown in Fig. 1, the unit normal vectors of the scattering plane (nOBA), incident meridian plane (nOCA), and scattering meridian plane (nOCB) can be obtained as nOBA ¼ 

Ω  Ω0 ; Ω  Ω0 

k  Ω0 ; nOCA ¼  k  Ω0 

kΩ ; nOCB ¼  k  Ω

ð6Þ

where k is the unit vector of the z direction. With the help of these vectors in Eqs. (6), i1 and i2 can be conveniently represented as i1 ¼ arccosðnOCA U nOBA Þ;

i2 ¼ arccosð  nOCB U nOBA Þ

nOBA U k 4 0;

ð7aÞ i1 ¼ arccosð  nOCA U nOBA Þ;

i2 ¼ arccosðnOCB U nOBA Þ

nOBA U k o 0:

ð7bÞ

z

If nOBA  k ¼0 or Ω  Ω0 ¼0, the two meridian planes coincide and are aligned with the scattering plane. Thus no rotation is needed, and Z ¼M. In Eq. (1), the extinction matrix κe of the medium is also a matrix with 4  4 elements. For randomly oriented particles with rotational symmetry, the extinction matrix can be simplified to a diagonal matrix as Eq. (8) 0 1 κe 0 0 0 B 0 κe 0 0 C B C κe ¼ B ð8Þ C; @ 0 0 κe 0 A

C S i (Ω')

S s (Ω) i1

i2

A

B Θ

O

x

0

y

Fig. 1. Geometry of scattering and meridian planes.

0

0

κe

where κe is the extinction coefficient of the medium, m  1. For the case with collimated radiation sources, the Stokes vector propagating in the medium consists of two parts: the collimated Stokes vector Sc and the scattered or

X. Ben et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 163 (2015) 50–62

2.2. Direct collocation meshless discretization

scattered nodes

Vr

r6

r4

d

r5

r

r2

r3

nw

r1

computional domain Fig. 2. The domain of the moving least-squares approximation.

diffuse intensity Sd, that is Sðr; ΩÞ ¼ Sc ðr; ΩÞ þ Sd ðr; ΩÞ:

ð9Þ

For the collimated case, it follows the law of Lambert  Beer, which can be expressed as Ω U∇Sc ðr; ΩÞ ¼  κe Sc ðr; ΩÞ:

ð10Þ

Substituting Eq. (9) into Eq. (1), the VRTE is changed to Ω U ∇½Sc ðr; ΩÞ þ Sd ðr; ΩÞ ¼  κe ½Sc ðr; ΩÞ þ Sd ðr; ΩÞ Z κs þ Zðr; Ω0 -ΩÞ½Sc ðr; Ω0 Þ þ Sd ðr; Ω0 ÞdΩ0 : 4π 4π

ð11aÞ Considering Eq. (10), then it can be simplified to Z κs Zðr; Ω0 -ΩÞSd ðr; Ω0 ÞdΩ0 Ω U∇Sd ðr; ΩÞ ¼ κe Sd ðr; ΩÞ þ 4π 4π þ

Nc κs X Zðr; Ωk -ΩÞSc ðr; Ωk Þ; 4π k ¼ 1

ð11bÞ

where Nc is the number of the incident directions for collimated radiation sources. In general, the boundary condition for the VRTE [Eq. (11b)], which takes into account the Lambertian (diffuse) reflection and collimated beam irradiation, can be written as Sðrw ; ΩÞ ¼ Sc ðrw ; ΩÞδðΩ Ωc Þ Z   1 R ðrw ; Ω0 -ΩÞSd ðrw ; Ω0 Þnw U Ω0 dΩ0 þ π nw U Ω0 o 0 d Nc X

1 þ π

53

  Rd ðrw ; Ωk -ΩÞSc ðrw ; Ωk Þnw U Ωk ;

For the application of direct collocation meshless discretization, the moving least-squares approximation is employed to construct the trial functions using collocated nodes, and the DCM is adopted to construct the linear equations. Tan and Zhao [26–28] have applied the DCM to solving the first order and two order scalar radiative transfer equations. In this paper, we extend this method to solve VRTE. For the differential term in Eq. (11b) corresponding to a scattered node of the domain of interest, a moving leastsquares approximation is used to build connections with the other scattered nodes around the node. It should be noted that for first-order linear partial differential equations such as Eq. (11b), the interpolation approximation is only applied for the differential terms. The detailed steps of completing moving least-squares approximations for element I of the Stokes vector are introduced as follows. As shown in Fig. 2, a spatial subdomain Vr, in the neighborhood of a node denoted by r, is considered as the domain of the definition of the moving least-squares approximation for the trial function at r. To approximate the distribution of a generic trial function I(r,Ω) in Vr over a number of local nodes {ri}, i¼ 1, 2, …, Nloc, the moving ̃ least-squares approximation I(r,Ω) of I(r, Ω), 8 rA Vr, is given by ~ ΩÞ ¼ Iðr;

p X

pj ðrÞaj ðr; ΩÞ ¼ pT ðrÞaðr; ΩÞ;

where p is the order of a Legendre polynomial, pj is the jth polynomial basis function, and aj is the corresponding expansion coefficient. The orthogonal vector pT can be chosen to be either linear or quadratic: pT ðrÞ ¼ ½1; x; y;

ð15aÞ

pT ðrÞ ¼ ½1; x; y; x2 ; xy; y2 :

ð15bÞ

In this article, the complete second order polynomial basis in Eq. (15b) is used for the moving least-squares approximation. The least squares values can then be calculated with the help of the local inner product as ~ ΩÞ  Iðr; ΩÞ‖2 min J r ðaÞ ¼ min‖Iðr; aj

aj

¼ min o aj

k¼1

p X

pj ðrÞaj  Iðr; ΩÞ;

j¼0

p X

pk ðrÞak Iðr; ΩÞ 4 r :

k¼0

ð16Þ

nw UΩk o 0 nw U Ω4 0;

ð12Þ

where δ is the Dirac delta function, and Rd is the reflection matrix for Lambertian surfaces with 4  4 elements. For an ideal Lambertian surface, this has the form as 0 1 ρ 0 0 0 B0 0 0 0C B C Rd ¼ B ð13Þ C; @0 0 0 0A 0

0

0

0

where ρ is the reflection coefficient of the wall.

ð14Þ

j¼0

The inner product 〈f, g〉 is defined as 〈f ; g〉r ¼

Nloc X

f ðri Þgðri Þwðjri rjÞ;

ð17Þ

i¼1

where w is a weight function with compact support. The basic idea of the moving least-squares approximation is to locally project the unknown fields into a prescribed function space using the least squares approach. Hence, it will yield different projections for different positions. By tuning the function space, the moving least-squares approximation can

54

X. Ben et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 163 (2015) 50–62

be very accurate. The most frequently used function space is the polynomial space. The local projection can be achieved through weight functions with compact support. In the following exposition, the function space spanned by the complete 2-D polynomial basis is second order, as Eq. (15b), so a weight function with more than second order derivatives is needed. In this paper, the fourth order Wendland's function [30] is chosen, which is compactly supported and continuous up to the fourth derivative. It is defined as ( ð1 r=dÞ6 ð35ðr=dÞ2 þ18ðr=dÞ þ3Þ ðr=dÞ r 1 wðr=dÞ ¼ ; 0 ðr=dÞ 41 ð18Þ where r is the distance between ri and r, and d is the diameter of the support domain given by d ¼ λðp þ1ÞΔh;

ð19Þ

where Δh is the local discrete length, and λ is the support domain amplifying factor. If the scattered nodes in the entire computational domain have an inhomogeneous distribution, an amplifying factor can improve the accuracy of the moving least-squares approximation. In this paper, we consider only scattered nodes with a uniform distribution, with a choice of λ¼1.0. A least squares procedure is employed to solve Eq. (16), and a solution to the minimization problem (using the extreme value condition ∂Jr(a)/∂aj ¼0), yields p X

o pk ðrÞ; pj ðrÞ 4 r aj ¼ opk ðrÞ; Iðr; ΩÞ 4 r ;

ð20Þ

j¼0

which can be written in matrix form as Φa ¼ f;

ð21Þ

where ΦðrÞ ¼ ½Φkj k ¼ 1…p;j ¼ 1…p ¼ ½〈pk ðrÞ; pj ðrÞ〉r j ¼ 1…p;k ¼ 1…p ; ð22aÞ aðr; ΩÞ ¼ ½aj j ¼ 1…p ¼ Φ

1

ðrÞfðr; ΩÞ;

fðr; ΩÞ ¼ ½f k k ¼ 1…p ¼ ½〈pk ðrÞ;

Iðr; ΩÞ〉r k ¼ 1…p ;

ð22bÞ ð22cÞ

Substituting Eq. (22b) into Eq. (14), the locally approximated I-element can be expressed as ~ ΩÞ ¼ pT ðrÞaðr; ΩÞ ¼ pT ðrÞΦ  1 ðrÞfðr; ΩÞ Iðr; ¼

Nloc X

pT ðrÞΦ  1 ðrÞpðri Þwðjri  rjÞIðri ; ΩÞ:

ð23Þ

i¼1

where Nloc is the number of local nodes in Vr. A nodal basis function (at node i) from the moving least-squares approximation can then be written as N i ðrÞ ¼ pT ðrÞΦ  1 ðrÞpðri Þwðjri rjÞ:

ð24Þ

The local polynomial approximation [Eq. (14)] can be written in the equivalent nodal basis expansion as ~ ΩÞ ¼ Iðr;

N loc X i¼1

Ni ðrÞIðri ; ΩÞ:

ð25Þ

The first-order partial derivatives of the I-element can be written as ~ ΩÞ ¼ ∇Iðr;

Nloc X

∇Ni ðrÞIðri ; ΩÞ

i¼1

¼

N loc X

∇pT ðrÞΦ  1 ðrÞpðri Þwðjri  rjÞIðri ; ΩÞ:

ð26Þ

i¼1

For the other three parameters of the Stokes vector, the nodal basis functions Ni(r) are the same as the I-element. Thus, the Stokes vector from the moving least-squares approximation has the form: ~ ΩÞ ¼ Sðr;

N loc X

Ni ðrÞSðri ; ΩÞ:

ð27Þ

i¼1

Substituting Eq. (26) into Eq. (11b), the VRTE from DCM using the moving least-squares approximation can be expressed as ΩU

N loc X

∇Ni ðrÞSd ðri ; ΩÞ ¼  κe Sd ðr; ΩÞ þ

i¼1

þ

κs 4π

Z 4π

Zðr; Ω0 -ΩÞSd ðr; Ω0 ÞdΩ0

Nc κs X Zðr; Ωk -ΩÞSc ðr; Ωk Þ: 4π k ¼ 1

ð28Þ

The Stokes vector of every node in the computational domain can be obtained by solving Eq. (28). If the total number of the nodes in the computational domain is m, Eq. (28) can be rewritten in matrix form as ASðΩÞ ¼ BðΩÞ

ð29Þ

where A is the stiffness matrix with m  m elements, which is determined by the nodal basis functions and the extinction coefficient of the medium. It is expressed as 0 1 a11 a12 ⋯ a1m B C a2m C B a21 a22 C ; ð30Þ A¼B B ⋮ ⋱ ⋮ C @ A am1 am2 … amm mm

where the specific form of ajk is written as ( jak ΩU ∇N k ðrj Þ ajk ¼ : κe þ Ω U∇Nk ðrj Þ j ¼ k

ð31Þ

Here, SðΩÞ is the Stokes vector matrix with m  4 elements and has the form 0 1 I d ðr1 ; ΩÞ Q d ðr1 ; ΩÞ U d ðr1 ; ΩÞ V d ðr1 ; ΩÞ B C B I d ðr2 ; ΩÞ Q d ðr2 ; ΩÞ U d ðr2 ; ΩÞ V d ðr2 ; ΩÞ C C SðΩÞ ¼ B B C ⋮ ⋮ ⋮ ⋮ @ A I d ðrm ; ΩÞ Q d ðrm ; ΩÞ U d ðrm ; ΩÞ V d ðrm ; ΩÞ m4

ð32Þ and B(Ω) is the source matrix with m  4 elements, which is determined by scattering sources, expressed as 0 1 b11 b12 b13 b14 B C B b21 b22 b23 b24 C C BðΩÞ ¼ B ; ð32Þ B ⋮ ⋮ ⋮ ⋮ C @ A bm1 bm2 bm3 bm4 m4

X. Ben et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 163 (2015) 50–62

where the specific form of bij is written as κs bij ¼ 4π

"

R

þ

4π

0

0

0

0

#

zj1 ðri ; Ω -ΩÞI d ðri ; Ω Þ þ zj2 ðri ; Ω -ΩÞQ d ðri ; Ω Þ þ dΩ zj3 ðri ; Ω0 -ΩÞU d ðri ; Ω0 Þ þ zj4 ðri ; Ω0 -ΩÞV d ðri ; Ω0 Þ

Nc zj1 ðri ; Ωk -ΩÞI c ðri ; Ωk Þ þ zj2 ðri ; Ωk -ΩÞQ c ðri ; Ωk Þ þ κs X 4π k ¼ 1 zj3 ðri ; Ωk -ΩÞU c ðri ; Ωk Þ þ zj4 ðri ; Ωk -ΩÞV c ðri ; Ωk Þ

! :

ð33Þ The Eq. (29) is a system of linear equations, and it is solved by the codes based on IMSL 7.0 (IMSL numerical libraries can be downloaded from http://www.roguewave. com/products-services/imsl-numerical-libraries). 2.3. Solution To solve the discretized VRTE Eq. (28) combined with boundary condition Eq. (12), the implementation of the DCM solution is given as follows. Step 1: Set the initial parameters in the simulation, and set an appropriate number of scattered nodes based on the requirements of the solution domain. Step 2: Choose an appropriate support domain amplifying factor λ, and then calculate the nodal basis functions of all nodes from the moving least-squares approximation. Step 3: Loop over each angular direction, and compute the system of linear equationsASðΩÞ ¼ BðΩÞ. (a) According to the nodal basis functions of all nodes, assemble a stiffness matrix A for the whole solution domain. (b) Impose boundary conditions, and replace the corresponding elements in the stiffness matrix A. If the node located at the boundary is the jth node in the computational domain, then the jth row of A is replaced by 0 elements, except for ajj ¼1. (c) Update the source matrix B(Ω). The new B(Ω) is assembled by the Stokes vector obtained from the previous iteration as Eq. (33). Then compute the system of linear equations to obtain the Stokes vector of each node. (d) Terminate the iteration process if the stopping criterion (defined as the maximum relative error being smaller than a given tolerance) is satisfied. Otherwise, go back to step (a). The convergence criterion for the iterations is given as ) (R  S dΩ  R S  n  4π n  1 dΩ max  4π R   100%o error   Sn dΩ

55

model is shown in Fig. 3, where a parallel solar beam penetrates into a scattering atmospheric layer from above, with a direction of μ0 ¼cos θ0 ¼0.2. The irradiation flux intensity of the beam is π. Along the z-axis, the optical thickness τ based on the height H of the atmosphere is 1.0, and the single scattering albedo ω is 0.99. The phase matrix used is that of Mie scattering at a wavelength of 0.951 μm from a gamma distribution of particles with effective radius of 0.2 μm, effective variance of 0.07, and refractive index n ¼1.44. The Legendre coefficients for the four unique elements of the phase matrix have been obtained by Evens and Stephens [31]. The bottom of the atmospheric layer has diffuse reflection with ρ¼ 0.1. For the one-dimensional model, the VRTE is simplified as μz

N loc X ∂N i ðrÞSd ðri ; ΩÞ

∂z

i¼1

¼  κe Sd ðr; ΩÞ þ þ

κs 4π

Z 4π

Zðr; Ω0 -ΩÞSd ðr; Ω0 ÞdΩ0

Nc κs X Zðr; Ωk -ΩÞSc ðr; Ωk Þ; 4π k ¼ 1

ð35Þ

where μz is the cosine of the direction with respect to the z-axis. To solve the one-dimensional problem, two columns of secondary nodes are needed to support the DCM solution, as shown in Fig. 3. The primary nodes have a uniform distribution along the z-axis. The simulation is run to obtain the angular distribution of Stokes parameters along μ¼cos θ at φ  φ0 ¼π⧸2, counterclockwise. At z/H¼0.5, due to the multiple scattering, all parameters of the Stokes vector are non-zero for both upward and downward directions. This case was first presented in Ref. [15] using the FN method, and it was compared by using MC methods In Ref. [32]. To investigate the performance of DCM, the Stokes vector curves obtained by the three methods are presented in Fig. 4. As displayed in Fig. 4, the Stokes parameters at different directions obtained by the DCM of this paper agree very well with those of Garcia and Siewert [15] using the FN method. It can be also seen that the precision of DCM is higher than that from the MC method [32], especially for the parameter V. As seen in Table 1, the CPU time when using DCM is only a few minutes, while it is more than ten hours

z S0 θ0

ð34Þ

4π

primary node

where n denotes the nth iteration.

3. Results and discussions 3.1. Validation of the model 3.1.1. Case 1: One-dimensional model of polarization by Mie scattering In this case, a classical ‘L¼13 problem’ tested by Garcia and Siewert [15] is investigated. The schematic of the

H secondary node

Fig. 3. Schematic of the oblique incident solar beam irradiating into a layer of atmosphere.

56

X. Ben et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 163 (2015) 50–62

1.2

×10−2

6.0 5.0

DCM FN

1.0

DCM FN

4.0

MC

MC

I

Q

0.8

×10−2

3.0

0.6 2.0 0.4

z/H=0.5

0.2 -1.0

0.0

z/H=0.5

1.0

-0.5

0.0 μ

0.5

0.0 -1.0

1.0

×10−2

3.0

-0.5

0.0 μ

V

U

-1.0 -1.5 -2.0 -2.5 -1.0

0.0 μ

0.5

DCM FN

0.0

MC

-1.0

-3.0

MC

-0.5

1.0

-2.0

DCM FN

z/H=0.5

1.0

×10−2

2.0 -0.5

0.5

-4.0 1.0

-5.0 -1.0

z/H=0.5

-0.5

0.0 μ

0.5

1.0

Fig. 4. Comparisons of Stokes parameters against data using the FN method [15].

Table 1 Computational time for different numbers of scattered nodes. The number of nodes (Ncol  Nrow)

Computational time

The number of solid angles (Nθ  Nφ)

35 3  11 3  21 3  31 3  41

2 min 3 min 1 min 2 min 2 min

40  40 40  40 40  40 40  40 40  40

24 s 5s 32 s 10 s 38 s

when using the MC method [32]. It can be concluded that DCM is much more efficient than the MC method. For DCM, the angular discretization is given as Nθ  Nφ ¼40  40, and the number of nodes is taken as Ncol  Nrow ¼ 3  41. The angular and spatial discretization values are demonstrated to be sufficient to obtain convergent results with good accuracy. To observe and analyze the influence of the nodes on the results, a node independence test is carried out. The number of primary nodes is taken to be 5, 11, 21, 31, and 41, respectively, and the angular discretization is taken to be Nθ  Nφ ¼40  40. Fig. 5 shows the curves for the four parameters I, Q, U and V of the Stokes vector at z/H ¼0.5 position. The curves for 21, 31 and 41 primary nodes agree well. The four parameters of the Stokes vector tend to

converge as the number of primary nodes increases. The maximum relative error between results using 31 and 41 nodes is less than 0.1%. To guarantee the precision of the results, 41 primary nodes are used for the next cases. In addition, Table 1 shows the computational times for different numbers of primary nodes. The computing environment we used is an Intel(R) Core™ i5-2320 CPU @ 3.00 Hz. It can be seen that the computational time increases with the increase of the number of nodes when the number of primary nodes is greater than 21. When the number of the nodes is too small, more iteration is needed to obtain convergent results, and consequently more computation time is required. 3.1.2. Case 2: One-dimensional model of polarization by Rayleigh scattering In this case, DCM is validated against the results by the MC method developed in Ref. [33] for a purely Rayleigh scattering atmospheric layer. A homogeneous planeparallel slab subjected to a collimated irradiation is considered. The optical thickness of the slab is τ¼ 0.5, and the scattering albedo ω is 1.0. The collimated radiation Stokes vector S0 ¼(π, 0, 0, 0)T covers the bottom of the slab with cosθ0 ¼0.4 and φ0 ¼01. A uniform angular discretization of Nθ  Nφ ¼60  24 is considered, and the number of nodes is taken to be Ncol  Nrow ¼3  30. The Mueller matrix for

X. Ben et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 163 (2015) 50–62

×10−2

6.0

Number of nodes 5×3 11× 3 21× 3 31× 3 41× 3

1.0

I

0.8 0.6

Number of nodes 5 ×3 11× 3 21× 3 31× 3 41× 3

4.0 3.0 2.0

0.4

1.0

z/H=0.5

0.2 -1.0

0.0

×10−2

5.0

Q

1.2

57

-0.5

0.0 μ

0.5

z/H=0.5

0.0 -1.0

1.0

×10−2

3.0

-0.5

0.0 μ

Number of nodes 5×3 11× 3 21× 3 31× 3 41× 3

1.0 0.0

U

Number of nodes 5×3 11× 3 21× 3 31× 3 41× 3

-1.5 -2.0

z/H=0.5

-2.5 -1.0

-0.5

0.0 μ

0.5

1.0

V

-1.0

1.0

×10−2

2.0 -0.5

0.5

-1.0 -2.0 -3.0 -4.0 -5.0 -1.0

z/H=0.5 -0.5

0.0 μ

0.5

1.0

Fig. 5. Comparisons of Stokes parameters with different numbers of nodes.

Rayleigh scattering is given by 0 cos 2 Θ þ 1 cos 2 Θ 1 B 3B cos 2 Θ  1 cos 2 Θ þ1 MðΘÞ ¼ B 4@ 0 0

2 cos Θ

0

0

0

2 cos Θ

0

0

0

0

0

1 C C C: A ð36Þ

Fig. 6(a) and (b) plots the values of I and Q Stokes parameters for the radiation exiting the boundaries at z¼ 0 and zmax. Three viewing azimuth angles 301, 901, and 1501 are chosen. It can be seen that the DCM results agree well with those obtained with the MC method. For this case, the CPU time taken by DCM is only 53 s. 3.1.3. Case 3 One-dimensional model of vector atmospheric radiative transfer To investigate the performance of the DCM code for solving VRTE under extreme scattering conditions, the two cases of vector radiative transfer in aerosols, and cloudy multiple scattering media, are considered. In these two cases, the phase functions are highly elongated in the forward direction, which are typical for atmospheric aerosols and clouds. Other elements of the scattering matrix oscillate considerably as a function of the scattering angle. The scattering matrices of atmospheric aerosols and clouds

appear in such an extreme situation that they can pose some difficulties to radiative transfer codes. The optical thicknesses of the aerosol and cloud layers are set to be 0.3262 and 5.0, respectively. The incident zenith angles for the two cases are both 601, and the viewing azimuth angles are 901, clockwise. Benchmark results have been provided by Ref. [25] and the elements of the scattering matrices can be downloaded at www.iup.physik.uni-bre men.de/  alexk. The benchmark data were obtained by the code SCIATRAN, which was based on a discrete ordinate method. In this work, the normalized Stokes vector is defined as S¼πSd/μ0F0, where F0 is the incident light irradiance at the area perpendicular to the incident light beam. For the two cases, the accuracy of the results is highly dependent on the number of streams (i.e. discretized directions). To obtain accurate results, the angular discretization for the aerosol scattering is taken to be Nθ  Nφ ¼180  360, while Nθ  Nφ ¼200  400 is adopted for the cloud scattering. A spatial discretization of Ncol  Nrow ¼3  41 is used for both cases. Fig. 7 shows the parameters of the normalized Stokes vector for the aerosol scattering and cloud scattering cases compared to the data from Ref. [25], and the curves for both reflected and transmitted lights are plotted in Fig. 7(a) and (b), respectively. It can be seen that the DCM results agree well

58

X. Ben et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 163 (2015) 50–62

4.5

2.0

4.0 1.5

3.5

1.0

2.5

Q

I

3.0

0.5

2.0 1.5

0.0

1.0

2.8 2.6 2.4 2.2 2.0 1.8 1.6 1.4 1.2 1.0 0.8 0.6 0.0

-0.8

-0.6

-0.4

-0.2

-0.5 -1.0

0.0

-0.8

-0.6

-0.4

-0.2

0.0

0.2

0.4

0.6

0.8

1.0

1.2 1.0 0.8 0.6 0.4 Q

I

0.5 -1.0

0.2 0.0 -0.2 0.2

0.4

0.6

0.8

1.0

-0.4 0.0

Fig. 6. Comparisons of I and Q Stokes parameters to data from Ref. [33]: (a) z¼ 0; (b) z¼ zmax.

with the benchmark data from Ref. [25], especially for the aerosol scattering case. The relative error of the I-element in the cloud scattering case is less than 3%. Under the same computing environment, the computational time for the aerosol scattering case is 28 h and 8 min, while the time for the cloud scattering case is 222 h and 6 min. The physical accuracy of DCM for different optical thicknesses was demonstrated in this case. The optical thickness for the aerosol layer is much bigger than that of the cloud layer. From a comparison of the two cases, it can be concluded that to improve accuracy, more streams (discretized directions) are needed for optically thick media. Meanwhile, much more CPU time is required for optically dense media. The two cases are both highly dependent on the number of streams, especially on the azimuth angles. Fig. 8 demonstrates the influence of azimuth discretization on the accuracy of the computational results. The number of discrete zenith angles is fixed to 180, while the numbers of the discrete azimuthal angles are set to 30, 60, 90, 180 and 360, respectively. The relative error of the Stokes parameters increases as Nφ decreases, and the V-element is more sensitive to Nφ. As Nφ decreases to a certain number such as 30, the V-element values will be positive, while the true values are negative. Due to the high dependence on azimuth

angles, sufficiently fine angular discretization is needed to obtain precise simulation results. 3.2. DCM applied to a two-dimensional model In this case, vector radiative transfers in a twodimensional system subjected to collimated irradiation are studied. As shown in Fig. 9, a collimated irradiation penetrates into a cuboid-shaped medium from the top. The length of the medium along the x-axis is over 10 times longer than that along the other two axes. The incident direction of the beam is cosθ0 ¼0.2 and φ0 ¼01 (parallel to x-axis). The beam with irradiation flux intensity π covers the entire top surface of the medium. The bottom surface is of diffuse reflection with ρ¼ 0.1. The optical thicknesses of the y-axis and z-axis directions are taken as τy ¼τz ¼1.0. For the two-dimensional model, the VRTE is expressed as

μy

N loc X ∂Ni ðrÞSd ðri ; ΩÞ i¼1

þ

κs 4π

Z 4π

∂y

þ μz

N loc X ∂Ni ðrÞSd ðri ; ΩÞ i¼1

Zðr; Ω0 -ΩÞSd ðr; Ω0 ÞdΩ0 þ

∂z

¼  κe Sd ðr; ΩÞ

Nc κs X Zðr; Ωk -ΩÞSc ðr; Ωk Þ: 4π k ¼ 1

ð37Þ

X. Ben et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 163 (2015) 50–62

5.0

59

2.0

DCM Ref.[25]

4.0

1.5

DCM Ref.[25]

Cloud 1.0

3.0

I

Q

Cloud 0.5

2.0 0.0

Aerosol 1.0 0.0

z/H=0 0

10

20

30

40

50

60

70

80

-1.0

90

0.8

0.2

0.6

0.0

0.4

20

30

40

50

60

70

80

90

Aerosol

-0.4

V

-0.2

Cloud

-0.6

Aerosol

-0.4

-0.8

-0.6

DCM Ref.[25]

-0.8

0

10

20

30

40

50

60

70

DCM Ref.[25]

-1.0

z/H=0 80

-1.2

90

0

0.0 4.0

Cloud

-0.2

DCM Ref.[25]

I 2.0

40

50

60

70

80

90

DCM Ref.[25] Aerosol

-0.6

-1.0

1.0

-1.2

z/H=1.0

-1.4 90

100 110 120 130 140 150 160 170 180

z/H=1.0

Cloud

100 110 120 130 140 150 160 170 180

1.0

0.0

0.0

-0.2

Aerosol -1.0

-0.6

-2.0

V

-0.4

Cloud

-0.8

-1.2 90

30

-0.8

Aerosol

-1.0

20

-0.4

Q

3.0

0.0 90

10

z/H=0

0.2

5.0

U

10

Cloud

0.0

U

0

Aerosol

-0.2

0.2

-1.0

z/H=0

-0.5

z/H=1.0

DCM Ref.[25]

100 110 120 130 140 150 160 170 180

Cloud

Aerosol -3.0 -4.0 -5.0 90

z/H=1.0

DCM Ref.[25]

100 110 120 130 140 150 160 170 180

Fig. 7. Comparisons of normalized Stokes vector elements for the aerosol scattering and cloud scattering cases to the data from Ref. [25]. (a) Reflected light. (b) Transmitted light.

As opposed to the one-dimensional model, all the scattered nodes are primary nodes with a uniform distribution in the domain, and the number of nodes is taken

to be Ncol  Nrow ¼41  41, which is sufficient to obtain convergent results. The angular discretization of Nμ  Nφ ¼ 40  40 is adopted.

60

X. Ben et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 163 (2015) 50–62

1.8

0.0

1.4

N =60

1.2

N =90

1.0

N =180

0.8

N =360

-0.2 -0.4

Q

N =30

I

1.6

-0.6

0.6

N =30

0.4

-0.8

0.2 0.0

90

-1.0

100 110 120 130 140 150 160 170 180

90

N =60

N =180

N =90

N =360

100 110 120 130 140 150 160 170 180

5.0

0.0

4.0 -0.2

3.0 2.0

N =30

-0.6

N =60

-0.8 -1.0 -1.2

90

1.0

V

U

-0.4

0.0 -1.0

N =90

-2.0

N =180

-3.0

N =360 100 110 120 130 140 150 160 170 180

N =30

-4.0 -5.0

90

N =60

N =180

N =90

N =360

100 110 120 130 140 150 160 170 180

Fig. 8. Comparisons of Stokes parameters with different number of discrete azimuth angles for the transmitted light of aerosol layer.

z

z θ0 S0

(1,1)

(0,1) y

y

x

(1,0)

(0,0) primary node

Fig. 9. Schematic of the oblique incident beam irradiating into a two-dimensional system.

Fig. 10 plots the curves of the Stokes vectors at three different positions (0.5,0.75), (0.5,0.5) and (0.5,0.25). The viewing azimuthal angle φ is 901, counterclockwise. The

trends of these curves are similar to those shown in Fig. 4 but the corresponding values are smaller. Due to the effect of the diffuse reflection at the bottom, the peak values are

X. Ben et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 163 (2015) 50–62

6.0

×10−2

61

×10−2 3.5

(y, z)=(0.5,0.75) (y, z)=(0.5,0.5) (y, z)=(0.5,0.25)

5.0 4.0

(y, z)=(0.5,0.75) (y, z)=(0.5,0.5) (y, z)=(0.5,0.25)

3.0 2.5

Q

I

2.0 3.0

1.5 2.0

1.0

1.0

0.5

0.0 -1.0

0.0 μ

0.5

0.0 -1.0

1.0

×10−2

1.5

-0.2

1.0

-0.4

0.5

-0.6

0.0

V

U

0.0

-0.5

-0.8

-1.2 -1.0

-0.5

0.0 μ

0.5

0.0 μ

0.5

1.0

×10−5

(y, z)=(0.5,0.75) (y, z)=(0.5,0.5) (y, z)=(0.5,0.25)

-0.5

(y, z)=(0.5,0.75) (y, z)=(0.5,0.5) (y, z)=(0.5,0.25)

-1.0

-0.5

-1.0

1.0

-1.5 -1.0

-0.5

0.0 μ

0.5

1.0

Fig. 10. Stokes parameters at three different points in the two-dimensional model.

shifted to the left (vertically down) at the positions away from the top of the two-dimensional system. The peak values also gradually decrease due to multiple scattering, except the V-element. The V-element curve shows a different tendency than the other three Stokes elements, which can be observed by noticing that the peak value increases at the positions away from the incident surface. This is because multiple scattering can strengthen the circular polarization. 4. Conclusion A direct collocation meshless method (DCM) was developed for solving vector radiative transfer in a scattering medium. Both one-dimensional and two-dimensional problems were examined. The numerical results showed that the present method is flexible and capable of solving for vector radiative transfers in multidimensional scattering media. With no fixed meshes, the scattered nodes can be flexibly configured according to the geometrical shape of the enclosure filled with a scattering medium, yet maintains good computational accuracy. The numerical model by DCM was validated by comparing the results for Mie scattering medium and in an atmospheric medium. The results compare very well with the data from previous papers. For the cases with large angular dependence, the

simulation results also show good agreement with the benchmark results. The results of the two-dimensional model demonstrate that the multiple scattering effects on the four Stokes parameters are different. With the enhancement of multiple scattering, the intensity (I element) and the linear polarization (Q and U elements) are weakened, while the circular polarization (V element) is strengthened.

Acknowledgement This work was supported by the National Natural Science Foundation of China (Grant nos. 51176040 and 51422602), and the Fundamental Research Funds for the Central Universities (Grant no. HIT. IBRSEM. A. 201413). A very special acknowledgement is made to the editors and referees who make important comments to improve this paper. References [1] Bryan AB, Ping Y, Andrew JH, Aaron B, Benjamin HC, Aronne M, et al. Ice cloud single-scattering property models with the full phase matrix at wavelengths from 0.2 to 100 μm. J Quant Spectrosc Radiat Transf 2014;46:123–39.

62

X. Ben et al. / Journal of Quantitative Spectroscopy & Radiative Transfer 163 (2015) 50–62

[2] Schepers D, Brugh JMJ, Hahne P, Butz A, Hasekamp OP, Landgraf J. LINTRAN v2.0: a linearized vector radiative transfer model for efficient simulation of satellite-bornnadir-viewing reflection measurements of cloudy atmospheres. J Quant Spectrosc Radiat Transf 2014;149:347–59. [3] Li L, Li Z, Wendisch M. Simulation of the influence of aerosol particles on stokes parameters of polarized skylight. Proceedings of the IOP conference series: earth and environmental science, vol. 17. IOP Publishing; 2014. [4] Hohmann A, Voit F, Schäfer J, Kienle K. Multiple scattering of polarized light: influence of absorption. Phys Med Biol 2014;59: 2583–97. [5] Li L, Li Z, Li K, Blarel K, Wendisch M. A method to calculate Stokes parameters and angle of polarization of skylight from polarized CIMEL sun/sky radiometers. J Quant Spectrosc Radiat Transf 2014;149:334–46. [6] Tuchin VV, Wang LV, Zimnyakov. Optical polarization in biomedical applications. Berlin: Springer; 2006. [7] Chandrasekhar S. Radiative transfer. London: Oxford University Press; 1950. [8] Zhao JM, Tan JY, Liu LH. On the derivation of vector radiative transfer equation for polarized radiative transport in graded index media. J Quant Spectrosc Radiat Transf 2012;113:239–50. [9] Freimanis J. Polarized radiative transfer equation in some curvilinear coordinate systems. J Quant Spectrosc Radiat Transf 2014;146: 250–69. [10] Wang XD, Wang LV, Sun CW, Yang CC. Polarized light propagation through scattering media: time-resolved Monte Carlo simulations and experiments. J Biomed Opt 2003;8:608–17. [11] Vaillona R, Wongb BT, Mengüç MP. Polarized radiative transfer in a particle-laden semi-transparent medium via a vector Monte Carlo method. J Quant Spectrosc Radiat Transf 2004;84:383–94. [12] Ramella-Roman JC, Prahl SA, Jacques SL. Three Monte Carlo programs of polarized light transport into scattering medium: part I. Opt Exp 2005;13:4420–38. [13] Li XY, Sun B, Yu YY. Reverse electric field Monte Carlo simulation for vector radiative transfer in the atmosphere. Chin Phys B 2014;23: 064219–1–064219-7. [14] Tynes HH, Kattawar GW, Zege EP, Katsev IL, Prikhach AS, Chaikovskaya LI. Monte Carlo and multicomponent approximation methods for vector radiative transfer by use of effective Mueller matrix calculations. Appl Opt 2001;40:400–12. [15] Garcia RDM, Siewert CE. The FN method for radiative transfer models that in include polarization effects. J Quant Spectrosc Radiat Transf 1989;41:117–45. [16] Sommersten ER, Lotsberg JK, Stamnes K, Stamnes JJ. Discrete ordinate and Monte Carlo simulations for polarized radiative transfer in a coupled system consisting of two medium with different refractive indices. J Quant Spectrosc Radiat Transf 2010;111:616–33. [17] RDM Garcia. Radiative transfer with polarization in a multi-layer medium subject to Fresnel boundary and interface conditions. J Quant Spectrosc Radiat Transf 2013;115:28–45.

[18] RDM Garcia. Fresnel boundary and interface conditions for polarized radiative transfer in a multilayer medium. J Quant Spectrosc Radiat Transf 2012;113:306–17. [19] Doicu A, Efremenko D, Trautmann T. A multi-dimensional vector spherical harmonics discrete ordinate method for atmospheric radiative transfer. J Quant Spectrosc Radiat Transf 2013;118:121–31. [20] Zhai PW, Hu YX, Charles RT, Lucker Patricia L. A vector radiative transfer model for coupled atmosphere and ocean systems based on successive order of scattering method. Opt Exp 2009;17:2057–79. [21] Zhai PW, Hu YX, Josset DB, Trepte CR, Lucker PL, Lin B. Advanced angular interpolation in the vector radiative transfer for coupled atmosphere and ocean systems. J Quant Spectrosc Radiat Transf 2013;115:19–27. [22] Zhai P.W., Hu, Y.X., Josset, D.B., Trepte, C.R., Lucker P.L., Lin B. Exact first order scattering correction for vector radiative transfer in coupled atmosphere and ocean systems. In: Proceedings of the conference on polarization-measurement, analysis, and remote sensing X, Baltimore; 2012. [23] Duan MZ, Min QL, Lu DR. A polarized radiative transfer model based on successive order of scattering. Adv Atmos Sci 2010;27:891–900. [24] Suniti S, Anthony BD, Annmarie E. vSmartMOM: A vector matrix operator method-based radiative transfer model linearized with respect to aerosol properties. J Quant Spectrosc Radiat Transf 2014;133:412–33. [25] Kokhanovsky AA, Budak VP, Cornet C, Duan MZ, Emde C, Katsev IL, et al. Benchmark results in vector atmospheric radiative transfer. J Quant Spectrosc Radiat Transf 2010;111:1931–46. [26] Tan JY, Zhao JM, Liu LH, Wang YY. Comparative study on accuracy and solution cost of the first/second-order radiative transfer equations using the meshless method. Numer Heat Transf, Part B 2009;55:324–37. [27] Tan JY, Zhao JM, Liu LH. Geometric optimization of a radiationconduction heating device using meshless method. Int J Therm Sci 2011;50:1820–31. [28] Zhao JM, Tan JY, Liu LH. A second order radiative transfer equation and its solution by meshless method with application to strongly inhomogeneous media. J Comput Phys 2013;232:431–55. [29] Luo K, Cao ZH, Yi HL, Tan HP. A direct collocation meshless approach with upwind scheme for radiative transfer in strongly inhomogeneous media. J Quant Spectrosc Radiat Transf 2014;135:66–80. [30] Wendland H. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Adv Comput Math 1995;4:389–96. [31] Evans KF, Stephens GL. A new polarized atmospheric radiative transfer model. J Quant Spectrosc Radiat Transf 1991;46:413–23. [32] Ben X, Yi HL, Tan HP. Polarized radiative transfer considering thermal emission in semitransparent media. Chin Phys B 2014;23: 099501–1–099501-9. [33] Vaillona R, Wongb BT, Mengüç MP. Polarized radiative transfer in a particle-laden semi-transparent medium via a vector Monte Carlo method. J Quant Spectrosc Radiat Transf 2004;84:383–94.