Просмотреть файл в отдельном окне: 9b7c1fda64a54db5a28ba810de782f24.pdf

Home

Search

Collections

Journals

About

Contact us

My IOPscience

Efficiency of algorithm for solution of vector radiative transfer equation in turbid medium slab

This article has been downloaded from IOPscience. Please scroll down to see the full text article.

2012 J. Phys.: Conf. Ser. 369 012021

(http://iopscience.iop.org/1742-6596/369/1/012021)

View the table of contents for this issue, or go to the journal homepage for more

Download details:

IP Address: 80.249.156.162

The article was downloaded on 03/07/2012 at 17:09

Please note that terms and conditions apply.

Eurotherm Conference No. 95: Computational Thermal Radiation in Participating Media IV

IOP Publishing

Journal of Physics: Conference Series 369 (2012) 012021

doi:10.1088/1742-6596/369/1/012021

Efficiency of algorithm for solution of vector radiative

transfer equation in turbid medium slab

V.P.Budak1,3, D.S.Efremenko2, O.V.Shagalov1

1

National research university “MPEI”, Krasnokazarmennaya, 14, Moscow, 111250,

Russia

2

Remote Sensing Technology Institute, German Aerospace Centre, Oberpfaffenhofen,

Wessling, Germany

BudakVP@mpei.ru

Abstract. The numerical solution of the vectorial radiative transfer equation (VRTE) is

possible only by its discretization, which requires elimination of the solution anisotropic part

including all the singularities. Discretized VRTE for the turbid medium slab has the unique

analytical solution in the matrix form. Modern packages of matrix (linear) algebra allow only

one possible algorithm of VRTE solution by computer. Various realizations of such an

algorithm differ by the method of the elimination of the solution anisotropic part. Methods of

the solution anisotropic part elimination are analysed in the paper. The codes created by the

authors of these methods are analysed in simple situations in order to define its influence on

the code efficiency. It is shown that the most effective method is based on the small angle

modification of the spherical harmonics method (MSH). The code based on MSH is

investigated in details by the influence of different properties of hard and software.

1. Boundary value problem

The comparison of different solution codes of scalar and vector radiative transfer equations (RTE,

VRTE) has shown that results of calculations are almost the same [1]. This suggests that different

computational methods of polarization fields in a turbid medium slab are in fact variants of the

uniform VRTE solution method [2, 3]. In this paper we present VRTE general solution, analysis and

optimization of the solution algorithm, and its effective computer implementation on the basis of the

general expression.

We have VRTE boundary problem for the turbid medium slab of optical depth ?0, on which a plane

wave of the arbitrary polarized light incidents in the direction ?l [3]:

0

r

? t t?? t

? ? r ? r ?

?

?

µ

?

+

?

=

?

?

?, ?l?) d?l?,

L(

,

l

)

L(

,

l

)

R(

)

x

(

l

,

l

)

R(

)L(

?? ??

4? ?

?r

r

r

?l ? ?l ), L(?, ?l )

?L(?, ?l )

=

L

?

(

= 0;

0

0

?= 0, ( z? , l? ) > 0

??

?=?0 , ( z? , ?l ) < 0

3

To whom any correspondence should be addressed.

Published under licence by IOP Publishing Ltd

1

(1)

Eurotherm Conference No. 95: Computational Thermal Radiation in Participating Media IV

IOP Publishing

Journal of Physics: Conference Series 369 (2012) 012021

doi:10.1088/1742-6596/369/1/012021

r

where L(?, ?l ) is the vector of Stokes parameters of the light field in the medium at the optical depth ?

in the direction ?l . We use Cartesian coordinate system OXYZ where axis OZ is directed down

perpendicularly to the slab border: ?l =

{ 1? µ

2

cos ?,

}

1 ? µ 2 sin ?, µ , ?l 0 =

{ 1? µ ,

2

0

}

0, µ0 ,

t

µ = (?l , z? ) , µ0 = (?l 0 , z? ) , z? is the unit vector along OZ. x (?l , ?l?) is the phase scattering matrix, ? is the

t

single scattering albedo. R(?) is the change matrix (rotator) of Stokes parameters at the reference

plane rotation. ? is the dihedral angle between two planes (z? ? ?l ) and (?l ? ?l?) , and ?' is the angle

between planes (?l ? ?l?) and (?l? ? z? ) . The unit vectors have a symbol “^”, the columns have the right

arrow, the rows have the left arrow, the matrices have the double arrow above the symbol.

For a numerical VRTE solution the integrals should be replaced with the finite sums that could be

done by one of two methods: the method of spherical harmonics (SH) and the discrete ordinates

method (DOM). The DOM is the best method for implementation since VRTE gains a clear ray

interpretation that essentially simplifies the definition of the composite boundary conditions. In this

case VRTE and boundary conditions take a form of interacting streams of radiation for the fixed

directions in space.

However in the initial boundary value problem (1) the replacement of the integral by the Gaussian

quadrature is impossible due to the solution singularity. The singularities are an inherent feature of the

ray approximation. Really, any break in the boundary conditions propagates inside the medium slab in

the form of the edge between light and shadow that gives the singularity in the radiance angular

distribution. Therefore, it is necessary to search the solution as the generalized function. To consider

the solution singularities, the desired vector of Stokes parameters is represented as the sum of two

parts: the anisotropic one that contains all the solution singularities, and the regular one – the smooth

function of angular variables:

r

r

r

(2)

L(?, ?l ) = L a (?, ?l ) + L r (?, ?l ) .

It is assumed that the anisotropic part can be found analytically. This allows expressing its integral

also analytically – a singularity elimination method. The substitution of (2) in (1) changes the

boundary value problem to the equivalent one

r

r

r

? ? L r (?, ?l ) r

? t t?? t

+ L r (?, ?l ) =

R(?) x (l , l?) R(??) L r (?, ?l? ) d?l? + ? (?, ?l );

??µ

?

??

4?

(3)

?r

r r

r

?

? L r (?, ?l )

= 0; L r (?, l )

= 0,

?= 0, µ> 0

?=?0 , µ< 0

??

with the source function on the right side of the equation

r

r

r

? L a (?, ?l ) r

? t t?? t

?

?

?

(4)

? (?, l ) =

R(?) x (l , l?) R(??) L a (?, l? ) dl? ? µ

? L a (?, ?l ) .

??

4? ?

2. Discretization of VRTE

The effectiveness of DOM application in the case of VRTE is strongly reduced since in the scalar type

it is determined by the possibility of the phase function presented by the series of the surface

harmonics. This allows reducing the double integral on the basis of the addition theorem to the single

t

one [4]. In the case of polarization the scattering matrix is surrounded by the rotator matrices R that

disturbs the azimuthal symmetry and does not allow using the addition theorem for the surface

harmonics. Ku??er-Ribari? [5] uses the circular basis to determine polarization

? L+2 ?

? Q ? iU ?

? 0 1 ?i 0 ? ? I ?

?

?

?

?

?

?? ?

r

t

t ?1

L

1 I ? V ? 1 ?1 0 0 ?1? ? Q ? t r

.

(5)

LCP = ? +0 ? = ?

=

= TCS L SP , TSC ? TCS

? L?0 ? 2 ? I + V ? 2 ?1 0 0 1 ? ?U ?

? ?

?

?

?

?? ?

?Q + iU ?

? 0 1 i 0 ? ?V ?

? L?2 ?

2

Eurotherm Conference No. 95: Computational Thermal Radiation in Participating Media IV

IOP Publishing

Journal of Physics: Conference Series 369 (2012) 012021

doi:10.1088/1742-6596/369/1/012021

k

It changes the rotator form and allows using generalizes surface functions Pmn

(cos ?) for the

scattering matrix representation that obey to the special form of the addition theorem

k ? ?

e ? im? Pmn

(l ? l?) e? in?? =

k

? (?1)

q =? k

q

k ?

Pmq

(l ? z? ) Pqnk (z? ? ?l?) eiq?? , ?? = ? ? ?? .

(6)

As the results, all the coefficients in VRTE become complex numbers, and this makes use of the

effective numerical methods for equation system solution difficult. Therefore, after applying the

addition theorem for the generalized spherical function one should return to the Stokes presentation

(SP).

Let’s consider the local transformation matrix of the ray at scattering in details:

t

t t

t

t t

t

t

t t

t

t

(7)

S(?l?, ?l ) ? R(?) x (?l?, ?l ) R(??) = TSC R CP (?) xCP (?l , ?l?) R CP (??) TCS ? TSC SCP (?l , ?l?) TCS ,

t t

t

t

t t

t

t ??

where x (l , l?) ? T x (?l , ?l?) T is the scattering matrix, R (?) = T R(?) T is the rotator in the

CP

CS

SP

SC

CP

CS

SC

circular polarization (CP) presentation.

Then we present the elements of scattering matrix in the form of the expansion on the generalized

spherical functions

K

t

(8)

[ xCP (cos ? )]rs = ? (2k + 1) xrsk Prk, s (cos ?) ,

k = max( r , s )

where r and s have only values ±0 and ±2; K is the number of harmonics in expansion of the scattering

matrix in the generalized spherical functions.

Taking into account the addition theorem (6), one can get every element of the local transformation

matrix in CP as

k

t

?K

?

?SCP (?l, ?l?) ? = ? ? (2k + 1) ? (?1)q eiq ( ???? ) xrsk Prqk (µ) Pqsk (µ?) ? .

(9)

?

? r ,s

q =? k

? k =0

?

t

With the use of diagonal matrix Ylq (µ) = Diag ( P?l 2,q (µ), P?l 0, q (µ), P+l 0,q (µ), P+l 2,q (µ) ) and taking into

account the reciprocal relation Pmk , n (µ) = Pnk,m (µ) one can get the expression (7) as

K

k

K

k

t

t

t t

t

t t

S(?l , ?l?) = ? (2k + 1) ? (?1) q eiq ( ???? ) TSC Ykq (µ) xk Ykq (µ?) TCS ? ? (2k + 1) ? eiq ( ???? ) Sqk (µ, µ?) ,

k =0

t

where xk ? ?? xrsk ?? .

Then

q =? k

k =0

(10)

q =? k

t

t t

t

t

t t

t t

Skq (µ, µ?) = (?1)q TSC Ykq (µ) xk Ykq (µ?) TCS = Pnk (µ)? k Pnk (µ?) ,

(11)

t t t

t

where ? k = TSC xk TCS is the so-called «Greek matrix»,

t

t t

t

t

t

(12)

Pnl (µ) = (?i ) n TSC Yln (µ) TCS ? PR + i PI ;

t

t

tl

PR and PI are the real parts of matrix Pn (µ) [6].

Finally, VRTE (1) could be written as

k t

r

r

r

t t

? r

? K

(13)

k

µ L r (?, ?l ) + L r (?, ?l ) =

(2

+

1)

Pnk (µ) ? eiq?? ? k Pnk (µ?) L r (?, ?l?)d?l? + ? (?, ?l ) .

?

?

??

4? k =0

q =? k

The scattering integral contains complex values though the whole expression (13) is real. With the

t

t

t

t

symmetry of the introduced functions PRn (µ) = PR? n (µ), PIn (µ) = ? PI? n (µ) one can eliminate the

imaginary part [6] in the equation (13). The obtained expression is correct for all scattering matrices.

However, in case of the aerosol block-diagonal matrix it could be simplified [6] more. One could

present the solution of problem (3) in the form [6]:

M

r

r

(14)

L r (?, ?l ) = ? (2 ? ? 0, m ) ???1 (m?)Lm1 (?, µ) + ?2 (m?)Lm2 (?, µ) ?? ,

m=0

that reduces VRTE to the following form

3

Eurotherm Conference No. 95: Computational Thermal Radiation in Participating Media IV

IOP Publishing

Journal of Physics: Conference Series 369 (2012) 012021

doi:10.1088/1742-6596/369/1/012021

r

1

r

r

t

? Lmc (?, µ) rm

t t

? K

µ

+ Lc (?, µ) = ? (2k + 1)? mk (µ)? k ? ? mk (µ?) Lmc (?, µ?)d µ? + ? (?, µ), c = 1, 2 ,

??

2 k =0

?1

where M is the term number in the Fourier series on the azimuth, m ? 0, M ,

? Q mk (µ)

0

0

0 ?

?

?

m

m

tk

0

R k (µ) ? Tk (µ)

0 ? ?1 (?) = diag ( cos ?, cos ?, sin ?, sin ? ) ,

?

? m (µ) =

,

? 0

? Tkm (µ) R mk (µ)

0 ? ?2 (?) = diag ( ? sin ?, ? sin ?, cos ?, cos ? ) ,

?

?

0

0

Q mk (µ) ??

?? 0

(15)

(16)

(k ? m)! m

Pk (µ) . (17)

(k + m)!

In order to get the solution of the obtained equation using DOM, we present the integrals entering

into the equation (15) taking into account the slab symmetry in the double Gaussian quadrature form.

Accordingly we have the system of ordinary differential equations for the fixed m and c (therefore we

omit them in the equation)

K

r

r

r

r

t

t

t t

d r±

?N2

(18)

µi±

Li (?) + L±i (?) = ? w j ? (2k + 1)? mk (µi± )? k ? mk (µ +j ) L+i (?) + ? km (µ ?j ) L?i (?) + ? i± (?) ,

4 j =1 k = 0

d?

r

r

r

r

where L±i (?) ? Lmc (?, µi± ), ? i± ( ?) ? ? cm (?, µi± ) , µ ±j = 0.5(? j ± 1) , ?j, wj are the zeroes and weights of the

R ln (µ) = 0.5i m ( Pmk ,2 (µ) + Pmk ,?2 (µ) ) , Tln (µ) = 0.5i m ( Pmk ,2 (µ) ? Pmk ,?2 (µ) ) , Qln (µ) =

(

)

Gaussian quadrature of the N/2 order.

Now we can introduce the following matrices and vectors

r

r

r

r

r

r

t

t ?

T

T

(19)

L ? ?? L + (?) L ? (?) ?? , ? ? ?? ? + (?) ? ? (?) ?? , M ? diag ( µi+ , ?µi+ ) , W ? diag ( wi , wi ) ,

4

t ?K

t

t t

?

(20)

A ? ? ? (2k + 1)? mk (µi+ )?k ? km (µ ±j ) ? ,

? k =0

?

r

T

where L ± (?) = ?? I (?, µ1± ), Q (?, µ1± ), U ( ?, µ1± ), V (?, µ1± ), K, I (?, µ ±N / 2 ), Q(?, µ ±N /2 ), U (?, µ ±N / 2 ), V (?, µ ±N /2 ) ?? .

Other arrays are arranged in the same way.

Since the regular part is a smooth function of angle, all the matrices are finite. This allows

rewriting the system (18) in the matrix form

r

tr

t r

t t t tt

d L(?)

= ? BL(?) + M ?1 ? (?), B ? M ?1 (1 ? AW) .

(21)

d?

3. Solution regular part

The solution of the obtained equation system has the analytic view [4, 7] in the matrix form

?0 t

t r

r

t r

B ?0

? L(0) + e L(?0 ) = ? e B ? M ?1 ? (?, µ 0 )d ? .

(22)

0

This expression (22) is equivalent to the solution representation that is the sum of the general

solution of homogeneous equation and the particular solution of the inhomogeneous one [8]. A matrix

representation is more convenient for the analytic transformation and implementation of solutions for

the computer. The solution of the homogeneous equation

t

t r

r

t

r

tt

d P(?)

? B ?0

(23)

L(?0 ) = e L(0) ? P(?0 ) L(0) :

= ? BP(?) ,

d?

connects the radiation at the slab lower boundary with an expression on the top and is called the propropagator [7].

4

Eurotherm Conference No. 95: Computational Thermal Radiation in Participating Media IV

IOP Publishing

Journal of Physics: Conference Series 369 (2012) 012021

doi:10.1088/1742-6596/369/1/012021

The problem of the solution (22) is connected with negative and positive exponents in the

expression that leads to fast worsening of the system matrix conditions. In order to eliminate this

effect, we used the scale transformation [4].

The matrix

exponent is represented in the form [4]

t

t r

r ??

B ?0

?1

0

(24)

e = Ue U ,

r

t

t

t t

where U is the matrix of eigenvectors of the matrix B and ? = diag(? ? , ? + ) is the diagonal matrix

t

t

of eigenvalues sorted by ascending ordering so that ? ? = ?? + .

Consequently, the equation (22) can be rewritten as

t

t

r

r

r

t

t

u12 ? ? L + (0) ? ? e?? ?0 ut 11 e?? ?0 ut 12 ? ? L + (?0 ) ? ? J ? ?

? u11

t

r

? ? ??t ? t

(25)

??r

? = ?r ? ,

?+? t

t

?? + ?0 t ? ?

+ 0

L

(

?

)

L

(0)

e

u

e

u

u

u

?

0

?

?

?

? ? J+ ?

?

21

22 ?

? 21

22

??

r

t

t

t t

r ? J? ? t ?0 ??

t ?1 r

t ?1 ? u11 u12 ?

?1

where J ? ? r ? = S ? e U M ? (?, µ 0 )d ? , U ? ? t

t ?.

? u 21 u 22 ?

0

? J+ ?

Note that expression (25) contains the exponent only with the positive power. Expressing in

r

r

equation (25) the emerging radiation L + (0), L ? (?0 ) from a layer through the incident radiation

r

r

L ? (0), L + (?0 ) we find a smooth regular part of the solution in the form of the so-called scatterers [7]:

r

r

r

t

t

? L ? (0) ? ? F? ? ? R ? T? ? ? L + (0) ?

t ? ?r

(26)

?r

? = ?r ? + ? t

?,

? L + (?0 ) ? ? F+ ? ? T+ R + ? ? L ? (?0 ) ?

t

t

r

r

t

t

t

t

t

t ?1

? F? ? t ? J ? ? ? R ? T? ? t ? u11

? e?? ?0 u12 ? t ? ? u12

e?? ?0 u11 ?

t ? = H? t t

where ? r ? = H ? r ? , ? t

? , H ? ? ??t ? t

? .

t

t

?? ?

? u 22 ??

u 21 ??

??e + 0 u 21

?? ? e + 0 u 22

? F+ ?

? J + ? ? T+ R + ?

The obtained system (26) is the desired solution that makes definition of the radiance distribution

of the reflected and transmitted radiations possible. Note that the expression (26) is the rigorous

analytical solution of the boundary value problem for VRTE discretized by DOM. The transfer from

VRTE (1) to the equation system (26) is possible due to the singularity elimination (2) and to the

change of the scattering integral (discretization) by the Gaussian quadrature (18).

The solution obtained in the form of scatterers (26) has a functional type and allows entering a

photometric concept of the layer radiance factor by the reflection or transmission. This immediately

gives the possibility to formulate a matrix-operator method to replace the two adjacent layers by the

equivalent one described by the expression identical to (32), but with effective parameters [3].

Consequently, the solution in the form of scatterers possesses the property of invariance.

Equation (26) can be reorganized in a form similar to the propagator (23):

t t t t

t t

r

r

? T+ ? R ? T??1R + R ? T??1 ? r

t ?1 t

t ?1 ? L(0) + F .

(27)

L(?0 ) = ?

T? ?

? ? T? R +

This kind of transformation has been called the “star product” [7]. Since one knows the differential

equation (23) for the propagator, one is able to get the one-point problem with initial conditions of the

matrix Riccati equation for the matrix elements of the scatterers [7], e.g. the reflection:

t

t

t

t

t t t

t t t t t

t t t t t

t ? b1

b2 ?

d R? t t t

d R+

= b 2 + b1R ? + R ? b1 + R ? b 2 R ? ;

= ? b 2 + b1R + ? R + b1 ? R + b 2 R + ; B ? ? t

t ? , (28)

d?

d?

?? ? b 2 ? b1 ??

t

and similarly for the transmission T± .

It is easy to see that (28) is an equivalent to the equations of Ambartsumyan-Chandrasekhar [8],

derived from the principle of invariant embedding. Thus, after eliminating the solution anisotropic part

and the equation discretization we reach a unique analytical solution in the matrix form (26). In this

sense the method of successive orders of scattering is an iterative method for finding the inverse

t

matrix A in (26), and the Monte Carlo is a stochastic method.

5

Eurotherm Conference No. 95: Computational Thermal Radiation in Participating Media IV

IOP Publishing

Journal of Physics: Conference Series 369 (2012) 012021

doi:10.1088/1742-6596/369/1/012021

Implementing the suggested solution (26) as the algorithm includes the calculation of the following

components: zeroes and weights of the quadrature formula for the VRTE discretization (18), source

function (4), eigenvectors and values of the system matrix (24), and products of matrices in (26).

The practical implementation of the indicated algorithm depends mainly on the size of matrices in

(26). Note that (26) is obtained for the regular part, therefore, its size is determined mostly not by the

medium optical parameters and the boundary conditions, but by the elimination of the solution

anisotropic part. In general, M?N?K. However, in case of correct elimination of the anisotropic part

solution, i.e. the smooth regular part is near to the isotropic angular distribution, it is possible that

M

Search

Collections

Journals

About

Contact us

My IOPscience

Efficiency of algorithm for solution of vector radiative transfer equation in turbid medium slab

This article has been downloaded from IOPscience. Please scroll down to see the full text article.

2012 J. Phys.: Conf. Ser. 369 012021

(http://iopscience.iop.org/1742-6596/369/1/012021)

View the table of contents for this issue, or go to the journal homepage for more

Download details:

IP Address: 80.249.156.162

The article was downloaded on 03/07/2012 at 17:09

Please note that terms and conditions apply.

Eurotherm Conference No. 95: Computational Thermal Radiation in Participating Media IV

IOP Publishing

Journal of Physics: Conference Series 369 (2012) 012021

doi:10.1088/1742-6596/369/1/012021

Efficiency of algorithm for solution of vector radiative

transfer equation in turbid medium slab

V.P.Budak1,3, D.S.Efremenko2, O.V.Shagalov1

1

National research university “MPEI”, Krasnokazarmennaya, 14, Moscow, 111250,

Russia

2

Remote Sensing Technology Institute, German Aerospace Centre, Oberpfaffenhofen,

Wessling, Germany

BudakVP@mpei.ru

Abstract. The numerical solution of the vectorial radiative transfer equation (VRTE) is

possible only by its discretization, which requires elimination of the solution anisotropic part

including all the singularities. Discretized VRTE for the turbid medium slab has the unique

analytical solution in the matrix form. Modern packages of matrix (linear) algebra allow only

one possible algorithm of VRTE solution by computer. Various realizations of such an

algorithm differ by the method of the elimination of the solution anisotropic part. Methods of

the solution anisotropic part elimination are analysed in the paper. The codes created by the

authors of these methods are analysed in simple situations in order to define its influence on

the code efficiency. It is shown that the most effective method is based on the small angle

modification of the spherical harmonics method (MSH). The code based on MSH is

investigated in details by the influence of different properties of hard and software.

1. Boundary value problem

The comparison of different solution codes of scalar and vector radiative transfer equations (RTE,

VRTE) has shown that results of calculations are almost the same [1]. This suggests that different

computational methods of polarization fields in a turbid medium slab are in fact variants of the

uniform VRTE solution method [2, 3]. In this paper we present VRTE general solution, analysis and

optimization of the solution algorithm, and its effective computer implementation on the basis of the

general expression.

We have VRTE boundary problem for the turbid medium slab of optical depth ?0, on which a plane

wave of the arbitrary polarized light incidents in the direction ?l [3]:

0

r

? t t?? t

? ? r ? r ?

?

?

µ

?

+

?

=

?

?

?, ?l?) d?l?,

L(

,

l

)

L(

,

l

)

R(

)

x

(

l

,

l

)

R(

)L(

?? ??

4? ?

?r

r

r

?l ? ?l ), L(?, ?l )

?L(?, ?l )

=

L

?

(

= 0;

0

0

?= 0, ( z? , l? ) > 0

??

?=?0 , ( z? , ?l ) < 0

3

To whom any correspondence should be addressed.

Published under licence by IOP Publishing Ltd

1

(1)

Eurotherm Conference No. 95: Computational Thermal Radiation in Participating Media IV

IOP Publishing

Journal of Physics: Conference Series 369 (2012) 012021

doi:10.1088/1742-6596/369/1/012021

r

where L(?, ?l ) is the vector of Stokes parameters of the light field in the medium at the optical depth ?

in the direction ?l . We use Cartesian coordinate system OXYZ where axis OZ is directed down

perpendicularly to the slab border: ?l =

{ 1? µ

2

cos ?,

}

1 ? µ 2 sin ?, µ , ?l 0 =

{ 1? µ ,

2

0

}

0, µ0 ,

t

µ = (?l , z? ) , µ0 = (?l 0 , z? ) , z? is the unit vector along OZ. x (?l , ?l?) is the phase scattering matrix, ? is the

t

single scattering albedo. R(?) is the change matrix (rotator) of Stokes parameters at the reference

plane rotation. ? is the dihedral angle between two planes (z? ? ?l ) and (?l ? ?l?) , and ?' is the angle

between planes (?l ? ?l?) and (?l? ? z? ) . The unit vectors have a symbol “^”, the columns have the right

arrow, the rows have the left arrow, the matrices have the double arrow above the symbol.

For a numerical VRTE solution the integrals should be replaced with the finite sums that could be

done by one of two methods: the method of spherical harmonics (SH) and the discrete ordinates

method (DOM). The DOM is the best method for implementation since VRTE gains a clear ray

interpretation that essentially simplifies the definition of the composite boundary conditions. In this

case VRTE and boundary conditions take a form of interacting streams of radiation for the fixed

directions in space.

However in the initial boundary value problem (1) the replacement of the integral by the Gaussian

quadrature is impossible due to the solution singularity. The singularities are an inherent feature of the

ray approximation. Really, any break in the boundary conditions propagates inside the medium slab in

the form of the edge between light and shadow that gives the singularity in the radiance angular

distribution. Therefore, it is necessary to search the solution as the generalized function. To consider

the solution singularities, the desired vector of Stokes parameters is represented as the sum of two

parts: the anisotropic one that contains all the solution singularities, and the regular one – the smooth

function of angular variables:

r

r

r

(2)

L(?, ?l ) = L a (?, ?l ) + L r (?, ?l ) .

It is assumed that the anisotropic part can be found analytically. This allows expressing its integral

also analytically – a singularity elimination method. The substitution of (2) in (1) changes the

boundary value problem to the equivalent one

r

r

r

? ? L r (?, ?l ) r

? t t?? t

+ L r (?, ?l ) =

R(?) x (l , l?) R(??) L r (?, ?l? ) d?l? + ? (?, ?l );

??µ

?

??

4?

(3)

?r

r r

r

?

? L r (?, ?l )

= 0; L r (?, l )

= 0,

?= 0, µ> 0

?=?0 , µ< 0

??

with the source function on the right side of the equation

r

r

r

? L a (?, ?l ) r

? t t?? t

?

?

?

(4)

? (?, l ) =

R(?) x (l , l?) R(??) L a (?, l? ) dl? ? µ

? L a (?, ?l ) .

??

4? ?

2. Discretization of VRTE

The effectiveness of DOM application in the case of VRTE is strongly reduced since in the scalar type

it is determined by the possibility of the phase function presented by the series of the surface

harmonics. This allows reducing the double integral on the basis of the addition theorem to the single

t

one [4]. In the case of polarization the scattering matrix is surrounded by the rotator matrices R that

disturbs the azimuthal symmetry and does not allow using the addition theorem for the surface

harmonics. Ku??er-Ribari? [5] uses the circular basis to determine polarization

? L+2 ?

? Q ? iU ?

? 0 1 ?i 0 ? ? I ?

?

?

?

?

?

?? ?

r

t

t ?1

L

1 I ? V ? 1 ?1 0 0 ?1? ? Q ? t r

.

(5)

LCP = ? +0 ? = ?

=

= TCS L SP , TSC ? TCS

? L?0 ? 2 ? I + V ? 2 ?1 0 0 1 ? ?U ?

? ?

?

?

?

?? ?

?Q + iU ?

? 0 1 i 0 ? ?V ?

? L?2 ?

2

Eurotherm Conference No. 95: Computational Thermal Radiation in Participating Media IV

IOP Publishing

Journal of Physics: Conference Series 369 (2012) 012021

doi:10.1088/1742-6596/369/1/012021

k

It changes the rotator form and allows using generalizes surface functions Pmn

(cos ?) for the

scattering matrix representation that obey to the special form of the addition theorem

k ? ?

e ? im? Pmn

(l ? l?) e? in?? =

k

? (?1)

q =? k

q

k ?

Pmq

(l ? z? ) Pqnk (z? ? ?l?) eiq?? , ?? = ? ? ?? .

(6)

As the results, all the coefficients in VRTE become complex numbers, and this makes use of the

effective numerical methods for equation system solution difficult. Therefore, after applying the

addition theorem for the generalized spherical function one should return to the Stokes presentation

(SP).

Let’s consider the local transformation matrix of the ray at scattering in details:

t

t t

t

t t

t

t

t t

t

t

(7)

S(?l?, ?l ) ? R(?) x (?l?, ?l ) R(??) = TSC R CP (?) xCP (?l , ?l?) R CP (??) TCS ? TSC SCP (?l , ?l?) TCS ,

t t

t

t

t t

t

t ??

where x (l , l?) ? T x (?l , ?l?) T is the scattering matrix, R (?) = T R(?) T is the rotator in the

CP

CS

SP

SC

CP

CS

SC

circular polarization (CP) presentation.

Then we present the elements of scattering matrix in the form of the expansion on the generalized

spherical functions

K

t

(8)

[ xCP (cos ? )]rs = ? (2k + 1) xrsk Prk, s (cos ?) ,

k = max( r , s )

where r and s have only values ±0 and ±2; K is the number of harmonics in expansion of the scattering

matrix in the generalized spherical functions.

Taking into account the addition theorem (6), one can get every element of the local transformation

matrix in CP as

k

t

?K

?

?SCP (?l, ?l?) ? = ? ? (2k + 1) ? (?1)q eiq ( ???? ) xrsk Prqk (µ) Pqsk (µ?) ? .

(9)

?

? r ,s

q =? k

? k =0

?

t

With the use of diagonal matrix Ylq (µ) = Diag ( P?l 2,q (µ), P?l 0, q (µ), P+l 0,q (µ), P+l 2,q (µ) ) and taking into

account the reciprocal relation Pmk , n (µ) = Pnk,m (µ) one can get the expression (7) as

K

k

K

k

t

t

t t

t

t t

S(?l , ?l?) = ? (2k + 1) ? (?1) q eiq ( ???? ) TSC Ykq (µ) xk Ykq (µ?) TCS ? ? (2k + 1) ? eiq ( ???? ) Sqk (µ, µ?) ,

k =0

t

where xk ? ?? xrsk ?? .

Then

q =? k

k =0

(10)

q =? k

t

t t

t

t

t t

t t

Skq (µ, µ?) = (?1)q TSC Ykq (µ) xk Ykq (µ?) TCS = Pnk (µ)? k Pnk (µ?) ,

(11)

t t t

t

where ? k = TSC xk TCS is the so-called «Greek matrix»,

t

t t

t

t

t

(12)

Pnl (µ) = (?i ) n TSC Yln (µ) TCS ? PR + i PI ;

t

t

tl

PR and PI are the real parts of matrix Pn (µ) [6].

Finally, VRTE (1) could be written as

k t

r

r

r

t t

? r

? K

(13)

k

µ L r (?, ?l ) + L r (?, ?l ) =

(2

+

1)

Pnk (µ) ? eiq?? ? k Pnk (µ?) L r (?, ?l?)d?l? + ? (?, ?l ) .

?

?

??

4? k =0

q =? k

The scattering integral contains complex values though the whole expression (13) is real. With the

t

t

t

t

symmetry of the introduced functions PRn (µ) = PR? n (µ), PIn (µ) = ? PI? n (µ) one can eliminate the

imaginary part [6] in the equation (13). The obtained expression is correct for all scattering matrices.

However, in case of the aerosol block-diagonal matrix it could be simplified [6] more. One could

present the solution of problem (3) in the form [6]:

M

r

r

(14)

L r (?, ?l ) = ? (2 ? ? 0, m ) ???1 (m?)Lm1 (?, µ) + ?2 (m?)Lm2 (?, µ) ?? ,

m=0

that reduces VRTE to the following form

3

Eurotherm Conference No. 95: Computational Thermal Radiation in Participating Media IV

IOP Publishing

Journal of Physics: Conference Series 369 (2012) 012021

doi:10.1088/1742-6596/369/1/012021

r

1

r

r

t

? Lmc (?, µ) rm

t t

? K

µ

+ Lc (?, µ) = ? (2k + 1)? mk (µ)? k ? ? mk (µ?) Lmc (?, µ?)d µ? + ? (?, µ), c = 1, 2 ,

??

2 k =0

?1

where M is the term number in the Fourier series on the azimuth, m ? 0, M ,

? Q mk (µ)

0

0

0 ?

?

?

m

m

tk

0

R k (µ) ? Tk (µ)

0 ? ?1 (?) = diag ( cos ?, cos ?, sin ?, sin ? ) ,

?

? m (µ) =

,

? 0

? Tkm (µ) R mk (µ)

0 ? ?2 (?) = diag ( ? sin ?, ? sin ?, cos ?, cos ? ) ,

?

?

0

0

Q mk (µ) ??

?? 0

(15)

(16)

(k ? m)! m

Pk (µ) . (17)

(k + m)!

In order to get the solution of the obtained equation using DOM, we present the integrals entering

into the equation (15) taking into account the slab symmetry in the double Gaussian quadrature form.

Accordingly we have the system of ordinary differential equations for the fixed m and c (therefore we

omit them in the equation)

K

r

r

r

r

t

t

t t

d r±

?N2

(18)

µi±

Li (?) + L±i (?) = ? w j ? (2k + 1)? mk (µi± )? k ? mk (µ +j ) L+i (?) + ? km (µ ?j ) L?i (?) + ? i± (?) ,

4 j =1 k = 0

d?

r

r

r

r

where L±i (?) ? Lmc (?, µi± ), ? i± ( ?) ? ? cm (?, µi± ) , µ ±j = 0.5(? j ± 1) , ?j, wj are the zeroes and weights of the

R ln (µ) = 0.5i m ( Pmk ,2 (µ) + Pmk ,?2 (µ) ) , Tln (µ) = 0.5i m ( Pmk ,2 (µ) ? Pmk ,?2 (µ) ) , Qln (µ) =

(

)

Gaussian quadrature of the N/2 order.

Now we can introduce the following matrices and vectors

r

r

r

r

r

r

t

t ?

T

T

(19)

L ? ?? L + (?) L ? (?) ?? , ? ? ?? ? + (?) ? ? (?) ?? , M ? diag ( µi+ , ?µi+ ) , W ? diag ( wi , wi ) ,

4

t ?K

t

t t

?

(20)

A ? ? ? (2k + 1)? mk (µi+ )?k ? km (µ ±j ) ? ,

? k =0

?

r

T

where L ± (?) = ?? I (?, µ1± ), Q (?, µ1± ), U ( ?, µ1± ), V (?, µ1± ), K, I (?, µ ±N / 2 ), Q(?, µ ±N /2 ), U (?, µ ±N / 2 ), V (?, µ ±N /2 ) ?? .

Other arrays are arranged in the same way.

Since the regular part is a smooth function of angle, all the matrices are finite. This allows

rewriting the system (18) in the matrix form

r

tr

t r

t t t tt

d L(?)

= ? BL(?) + M ?1 ? (?), B ? M ?1 (1 ? AW) .

(21)

d?

3. Solution regular part

The solution of the obtained equation system has the analytic view [4, 7] in the matrix form

?0 t

t r

r

t r

B ?0

? L(0) + e L(?0 ) = ? e B ? M ?1 ? (?, µ 0 )d ? .

(22)

0

This expression (22) is equivalent to the solution representation that is the sum of the general

solution of homogeneous equation and the particular solution of the inhomogeneous one [8]. A matrix

representation is more convenient for the analytic transformation and implementation of solutions for

the computer. The solution of the homogeneous equation

t

t r

r

t

r

tt

d P(?)

? B ?0

(23)

L(?0 ) = e L(0) ? P(?0 ) L(0) :

= ? BP(?) ,

d?

connects the radiation at the slab lower boundary with an expression on the top and is called the propropagator [7].

4

Eurotherm Conference No. 95: Computational Thermal Radiation in Participating Media IV

IOP Publishing

Journal of Physics: Conference Series 369 (2012) 012021

doi:10.1088/1742-6596/369/1/012021

The problem of the solution (22) is connected with negative and positive exponents in the

expression that leads to fast worsening of the system matrix conditions. In order to eliminate this

effect, we used the scale transformation [4].

The matrix

exponent is represented in the form [4]

t

t r

r ??

B ?0

?1

0

(24)

e = Ue U ,

r

t

t

t t

where U is the matrix of eigenvectors of the matrix B and ? = diag(? ? , ? + ) is the diagonal matrix

t

t

of eigenvalues sorted by ascending ordering so that ? ? = ?? + .

Consequently, the equation (22) can be rewritten as

t

t

r

r

r

t

t

u12 ? ? L + (0) ? ? e?? ?0 ut 11 e?? ?0 ut 12 ? ? L + (?0 ) ? ? J ? ?

? u11

t

r

? ? ??t ? t

(25)

??r

? = ?r ? ,

?+? t

t

?? + ?0 t ? ?

+ 0

L

(

?

)

L

(0)

e

u

e

u

u

u

?

0

?

?

?

? ? J+ ?

?

21

22 ?

? 21

22

??

r

t

t

t t

r ? J? ? t ?0 ??

t ?1 r

t ?1 ? u11 u12 ?

?1

where J ? ? r ? = S ? e U M ? (?, µ 0 )d ? , U ? ? t

t ?.

? u 21 u 22 ?

0

? J+ ?

Note that expression (25) contains the exponent only with the positive power. Expressing in

r

r

equation (25) the emerging radiation L + (0), L ? (?0 ) from a layer through the incident radiation

r

r

L ? (0), L + (?0 ) we find a smooth regular part of the solution in the form of the so-called scatterers [7]:

r

r

r

t

t

? L ? (0) ? ? F? ? ? R ? T? ? ? L + (0) ?

t ? ?r

(26)

?r

? = ?r ? + ? t

?,

? L + (?0 ) ? ? F+ ? ? T+ R + ? ? L ? (?0 ) ?

t

t

r

r

t

t

t

t

t

t ?1

? F? ? t ? J ? ? ? R ? T? ? t ? u11

? e?? ?0 u12 ? t ? ? u12

e?? ?0 u11 ?

t ? = H? t t

where ? r ? = H ? r ? , ? t

? , H ? ? ??t ? t

? .

t

t

?? ?

? u 22 ??

u 21 ??

??e + 0 u 21

?? ? e + 0 u 22

? F+ ?

? J + ? ? T+ R + ?

The obtained system (26) is the desired solution that makes definition of the radiance distribution

of the reflected and transmitted radiations possible. Note that the expression (26) is the rigorous

analytical solution of the boundary value problem for VRTE discretized by DOM. The transfer from

VRTE (1) to the equation system (26) is possible due to the singularity elimination (2) and to the

change of the scattering integral (discretization) by the Gaussian quadrature (18).

The solution obtained in the form of scatterers (26) has a functional type and allows entering a

photometric concept of the layer radiance factor by the reflection or transmission. This immediately

gives the possibility to formulate a matrix-operator method to replace the two adjacent layers by the

equivalent one described by the expression identical to (32), but with effective parameters [3].

Consequently, the solution in the form of scatterers possesses the property of invariance.

Equation (26) can be reorganized in a form similar to the propagator (23):

t t t t

t t

r

r

? T+ ? R ? T??1R + R ? T??1 ? r

t ?1 t

t ?1 ? L(0) + F .

(27)

L(?0 ) = ?

T? ?

? ? T? R +

This kind of transformation has been called the “star product” [7]. Since one knows the differential

equation (23) for the propagator, one is able to get the one-point problem with initial conditions of the

matrix Riccati equation for the matrix elements of the scatterers [7], e.g. the reflection:

t

t

t

t

t t t

t t t t t

t t t t t

t ? b1

b2 ?

d R? t t t

d R+

= b 2 + b1R ? + R ? b1 + R ? b 2 R ? ;

= ? b 2 + b1R + ? R + b1 ? R + b 2 R + ; B ? ? t

t ? , (28)

d?

d?

?? ? b 2 ? b1 ??

t

and similarly for the transmission T± .

It is easy to see that (28) is an equivalent to the equations of Ambartsumyan-Chandrasekhar [8],

derived from the principle of invariant embedding. Thus, after eliminating the solution anisotropic part

and the equation discretization we reach a unique analytical solution in the matrix form (26). In this

sense the method of successive orders of scattering is an iterative method for finding the inverse

t

matrix A in (26), and the Monte Carlo is a stochastic method.

5

Eurotherm Conference No. 95: Computational Thermal Radiation in Participating Media IV

IOP Publishing

Journal of Physics: Conference Series 369 (2012) 012021

doi:10.1088/1742-6596/369/1/012021

Implementing the suggested solution (26) as the algorithm includes the calculation of the following

components: zeroes and weights of the quadrature formula for the VRTE discretization (18), source

function (4), eigenvectors and values of the system matrix (24), and products of matrices in (26).

The practical implementation of the indicated algorithm depends mainly on the size of matrices in

(26). Note that (26) is obtained for the regular part, therefore, its size is determined mostly not by the

medium optical parameters and the boundary conditions, but by the elimination of the solution

anisotropic part. In general, M?N?K. However, in case of correct elimination of the anisotropic part

solution, i.e. the smooth regular part is near to the isotropic angular distribution, it is possible that

M

Чтобы не видеть никакую рекламу на сайте, нужно стать

**VIP**-пользователем.

Это можно сделать совершенно бесплатно. Читайте подробности тут.