Polarized radiative transfer in an arbitrary multilayer semitransparent medium Xun Ben,1 Hong-Liang Yi,1,2 and He-Ping Tan1,3 1
School of Energy Science and Engineering, Harbin Institute of Technology, 92 West Dazhi Street, Harbin 150001, China 2
e-mail: [email protected] 3
e-mail: [email protected]
Received 27 September 2013; revised 16 January 2014; accepted 20 January 2014; posted 31 January 2014 (Doc. ID 198456); published 27 February 2014
Polarized radiative transfer in a multilayer system is an important problem and has wide applications in various fields. In this work, a Monte Carlo (MC) model is developed to simulate polarized radiative transfer in a semitransparent arbitrary multilayer medium with different refractive indices in each layer. Two kinds of polarization mechanisms are considered: scattering by particles and reflection and refraction at the Fresnel surfaces or interfaces. The MC method has an obvious superiority in that complex mathematical derivations can be avoided in solving the polarization by Fresnel reflection and refraction in an arbitrary multilayer system. We define the vector radiative transfer matrix (VRTM), which describes the polarization characteristics of radiative transfer, and obtain four elements of Stokes vector using the VRTM. The results for the two-layer model by MC method are compared against those for coupled atmosphere–ocean model by the discrete–ordinate method available in the literature, which validates the correctness of the MC multilayer model of polarized radiative transfer. Finally, the results for three-layer, five-layer, and ten-layer models are presented in graphical form. Results show that in the multilayer system, total reflections occurring at the surfaces/interfaces have significant effects on the polarized radiative transfer, which causes abrupt changes or fluctuations like waves in the curves of the Stokes vector. © 2014 Optical Society of America OCIS codes: (030.5620) Radiative transfer; (230.4170) Multilayers; (290.4210) Multiple scattering; (290.7050) Turbid media. http://dx.doi.org/10.1364/AO.53.001427
1. Introduction
Polarized radiative transfer has important applications in atmospheric and ocean optics, remote sensing [1,2], biomedical optics [3,4], particle characterization [5], and many other engineering fields and research disciplines. Over more than 40 years, many numerical methods for vector radiative transfer have been proposed and developed, including the Monte Carlo (MC) method [6], F N method [7], adding– doubling method [8], discrete–ordinate (DO) method [9], successive order scattering (SOS) method [10], 1559-128X/14/071427-15$15.00/0 © 2014 Optical Society of America
Markov chain method [11], and so on. Among these methods, the MC method is very popular for polarized radiative transfer, especially in the 21st century [12–19]. The MC method is time consuming, and in order to overcome this disadvantage, some other numerical methods have been developed for solving the polarized radiative transfer. Garcia and Siewert [7] used the F N method to solve the complete and general polarization problem for a plane–parallel slab, and tabulated the results for the problem of L 13 and L 60. Evans and Stephens [8] built a polarized plane–parallel radiative transfer model for randomly oriented particles of arbitrary shape using a direct method combined with the doubling and adding method. Siewert [9] developed a DO 1 March 2014 / Vol. 53, No. 7 / APPLIED OPTICS
1427
method for polarized radiative transfer in a one-dimensional atmosphere of randomly orientated oblate spheroids. Lenoble et al. [10] presented a SOS code to solve the equation of vector radiative transfer in a plane–parallel earth atmosphere, and considered polarizing effects caused by aerosol scattering and reflectance by water or terrestrial surfaces. Xu et al. [11] developed the Markov chain method to simulate polarized radiative transfer in a pseudospherical atmosphere with solar illumination. Kokhanovsky et al. [20] 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, and presented the benchmark data for one-dimensional vector atmospheric radiative transfer. Most of the work focused on the polarized radiative transfer in scattering media with uniform refractive index. However, vector radiative transfer in a medium with spatial changes in the refractive index is a topic of importance in many areas of investigation. An atmosphere–ocean two-layer model of polarized radiative transfer is a classical and important problem. Kattawar and co-workers [21,22] first examined the problem using the MC method. They calculated the degree and direction of polarization, the ellipticity, and the radiance of the radiation at various levels in the atmosphere–ocean system [21], and obtained the complete four-component Stokes vector at any region in a fully inhomogeneous atmosphere– ocean system with inclusion of a stochastic interface [22]. After that, especially in the past five years, this problem has attracted the attention of many scholars. Zhai and co-workers [23,24] developed the SOS method for the problem. In 2013, Zhai et al. [25] applied the Iteration of Source Matrix method and single scattering method to the SOS method to solving the problem, and obtained the polarized radiation field at arbitrary viewing angles. He et al. [26,27] used the matrix–operator method for vector radiative transfer in the atmosphere–ocean system. Ota et al. [28] developed a solution, based on the DO and matrix–operator methods, to the problem of vector radiative transfer in the coupled atmosphere– ocean system. Sommersten et al. [29] presented an algorithm based on DO method for vector radiative transfer in a coupled system of atmosphere and ocean with different refractive indices, and compared the results by the DO method with those by the MC method. In the past two years, Garcia [30] investigated the vector radiative transfer in a multilayer medium with Fresnel boundaries or interfaces using DO method. His derivation of the reflection and transmission matrices gave rise to some arguments [31], and a consensus on this issue has not been reached [32]. In 2013, Garcia [33] applied the analytical discrete ordinates (ADO) method to calculate the polarized radiative transfer in a multilayer medium with Fresnel boundaries and interfaces, and tabulated the results of four Stokes elements for a nine-layer system. Up to now, there has been 1428
APPLIED OPTICS / Vol. 53, No. 7 / 1 March 2014
little work reported on polarized radiative transfer in a multilayer system with two or more interfaces considering different refractive indices in each layer. Although the MC method is time consuming, it has an obvious superiority over other numerical methods in that complex mathematical derivations can be avoided for solving the polarization problem caused by Fresnel reflection and refraction in an arbitrary multilayer system. In this work, we develop a MC model of vector radiative transfer in an arbitrary multilayer medium with consideration of Fresnel or Lambertian surfaces/interfaces. In different layers, not only can scattering type be different but physical parameters such as refractive index, scattering albedo, extinction coefficient, and so on can also be set to different values. In Section 2, we discuss the basic theory of vector radiative transfer, mainly including the descriptions for scattering, reflection, and transmission of photons. We will focus on how to determine scattering directions using the improved reject method in the MC simulation and the definition of vector radiative transfer matrix (VRTM). In Section 3, since most numerical results available in the literature are for two-layer models of the atmosphere–ocean system, we compare our results with those for atmosphere–ocean systems from the literature to verify our MC method multilayer model. Finally, we investigate polarized radiative transfer in three-layer, five-layer, and ten-layer models, respectively, and analyze the effects of total reflection at the Fresnel surfaces/interfaces on the results of Stokes vector. 2. Descriptions of the Problem
We consider the multilayer medium with an oblique incidence of a beam of parallel light on the top surface, as shown in Fig. 1. The medium includes K layers, and every layer has its own independent physical properties. In Fig. 1, nk is the refractive index of the k layer, and τk is the optical thickness of this layer; the refractive indices outside the top and bottom surfaces of the multilayer medium are S0 nK+1
θ0
nK
Layer K
τΚ
nk
Layer k
τκ
n2
Layer 2
τ2
n1
Layer 1
τ1
n0
τ
Fig. 1. Geometry of the multilayer problem.
denoted by nk1 and n0 , respectively. There are two surfaces for each layer except for layers nk1 and n0 , called top surface and bottom surface. For layer 0, there is only one top surface, and for layer K 1, there is only one bottom surface. Two boundary surfaces and each interface between the adjacent two layers are semitransparent. Two types of the interfaces/boundaries are considered in this model, and each interface or boundary can be Lambertian one or Fresnel one. The Stokes vector S I; Q; U; V0 is usually used to describe the intensity and polarized states of light. The specific form of Stokes vector is given by 0 1 1 El El Er Er I B Q C 1 B El E − Er Er C l B C C SB @ U A η @ El Er Er E A; l iEl Er − Er El V 0
(1)
where η μm ∕ε1∕2 ; μm is the magnetic permeability, and ε is the electric permittivity of the medium; El and Er are the complex amplitudes of the electric field along the two vectors which are parallel and perpendicular to a reference plane; “∗”denotes the complex conjugate. In this work, we use the MC technique to solve the polarized radiative transfer in the arbitrary multilayer medium as shown in Fig. 1. Note that we consider the cold medium in the present work and ignore the thermal emission of the medium. A.
Scattering
The photons are tracked in the multilayer medium, which experiences free propagation, scattering, reflection, transmission, and other events until the photons escape from the boundaries or are absorbed. The deterministic behavior of light is determined by the random behavior of a finite number of photons in the medium, and each photon randomly propagates along the possible path. The distance a photon travels without being absorbed, scattered, reflected, or refracted is determined by L−
ln1 − RL ; κ
photon will be absorbed by the medium and not tracked again; otherwise the photon will be scattered by the particle. The scattering problem for polarized radiative transfer is much more complicated than that for the scalar one. Chandrasekhar [34] first put forward the idea of photon scattering from one location to the next in the meridian and scattering planes. Figure 2 shows the geometry of the scattering problem. The directions of the photon propagation before and after scattering are represented as those pointed to by the lines OA and OB, respectively, on the unit sphere. The incident Stokes vector Si and the resulting scattering Stokes vector Ss determine a plane AOB, defined as the scattering plane. The vector of incident Stokes vector Si and the z axis determine a plane AOC, called the meridian plane. When a scattering event occurs, a new meridian plane is created by the new direction of propagation Ss (scattering Stokes vector) and the z axis, and the Stokes vector must be transformed so that the polarization is properly represented with respect to the new meridian plane BOC. Two rotations of the Stokes vector, in and out of the scattering plane, are required in the MC method. Based on the above analysis, the resulting Stokes vector after scattering is Ss Lπ − i2 MΘL−i1 Si ;
(3)
where the matrix M is the single scattering Mueller matrix defined in Van de Hulst [35], and L is a rotation matrix; i1 denotes the angle of rotation into the scattering plane AOB from the meridian plane AOC, and i2 denotes the angle of radiation into the meridian plane BOC from the scattering plane AOB. The term of Lπ − i2 MΘL−i1 is called the phase matrix. The matrix M has 4 × 4 elements 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 one that has eight nonzero elements [36], of which only four are independent. The form is
(2)
where L is the free distance, κ is the extinction coefficient, and RL is a pseudorandom number with a uniform distribution. The following MC analysis of multiple scattering events of the photons in a semitransparent medium is based on the assumption that the scattering event is independent and does not consider coherence effects. When the photons propagate in each layer of the medium, they may probably interact with the particles in the medium. The photons may be absorbed by the medium or scattered to other directions by the particles. A pseudorandom number is compared with the scattering albedo ω of the medium to determine whether the photon is absorbed or scattered. If the pseudorandom number Rω > ω the
z
Ss
C
Si
i2 B
θ2
Θ
θ1
i1 A
O
y
ϕ2 ϕ1 ϕ2
ϕ1
x
Fig. 2. Geometry of scattering plane and meridian planes. The photon’s direction of propagation before and after scattering is θ1 ; φ1 and θ2 ; φ2 , respectively. 1 March 2014 / Vol. 53, No. 7 / APPLIED OPTICS
1429
0
M1 B M2 B MΘ @ 0 0
M2 M1 0 0
0 0 M3 −M 4
1 0 0 C C; M4 A M3
(4a)
where Θ is the scattering angle between the incident and scattering directions. The simplified Mueller matrix Eq. (4a) can be used for Mie scattering and pure Rayleigh scattering as well. For Rayleigh scattering, the matrix M can also be written by 0
cos2 Θ 1 cos2 Θ − 1
0
2 2 0 3Δ B B cos Θ − 1 cos Θ 1 MΘ B 4 @ 0 0 2 cos Θ
0 0 01 0 0 01
0
1
0
C C C A
0 0 2Δ0 cos Θ
B0 0 0 0C C 1 − ΔB @ 0 0 0 0 A;
(4b)
0 0 0 0 in which Δ 1 − δ∕1 δ∕2 and Δ0 1 − 2δ∕ 1 − δ, where δ is the depolarization factor [37] which is equal to zero for pure Rayleigh scattering. Equation (4b) can be summarized as 0 1 0 M1 M2 0 B M2 M3 0 0 C C: MΘ B (4c) @ 0 0 M4 0 A 0 0 0 M5 The matrix L represents a rotation in the clockwise direction with respect to an observer looking into the direction of the photon propagation [38], and can be written as 0 1 1 0 0 0 B 0 cos 2ϕ sin 2ϕ 0 C C Lϕ B (5) @ 0 − sin 2ϕ cos 2ϕ 0 A: 0 0 0 1 The physical meaning of Eq. (3) is easy to understand. First of all, the incident Stokes vector Si is rotated at an i1 angle from the meridian plane AOC to the scattering plane AOB counterclockwise. Then, the incident Stokes vector Si interacts with the particle and is scattered in a new direction. Finally, the scattering Stokes vector is rotated at a π − i2 angle from the scattering plane AOB to the meridian plane BOC clockwise. Eventually, the scattering Stokes vector Ss is obtained. For a given incident direction θ1 ; φ1 as shown in Fig. 2, if the angles i1 and Θ are determined, then the angles i2, θ2 , and φ1 − φ2 can be calculated by the spherical laws of sines and cosines [39]. In order to obtain these angles, the cumulative distribution method [19] was used. Whitney [40] presented a rejection method to obtain i1 and Θ by using a brute-force numerical method in which all angles are calculated and compared with each other. The 1430
APPLIED OPTICS / Vol. 53, No. 7 / 1 March 2014
brute-force numerical method is time consuming to find i1 and Θ. In this work, we develop an optimized rejection method used for obtaining i1 and Θ, which obviously saves the computing time. The detailed steps are given as follows. Before we rotate the scattering Stokes vector from the scattering plane to the meridian plane BOC, Eq. (3) can be rewritten as S0s MΘL−i1 Si :
(6a)
Then the I Stokes element can be solved from Eq. (4a) (both for Mie and Rayleigh scattering), giving I 0s I i M 1 Θ Qi M 2 Θ cos 2i1 − U i M 2 Θ sin 2i1 ; (6b) where Θ ∈ 0; π and i1 ∈ 0; 2π. By using a trigonometric function formula, Eq. (6b) can be rewritten as I 0s I i M 1 Θ
q Q2i U 2i M 2 Θ cos2i1 ϕ; (6c)
where ϕ is an angle we are not interested in. Regarding i1 as independent variable and considering the maximum absolute value of the cosine function, we can obtain the maximum of I 0s;Θ, I 0s;Θ;max I i M 1 Θ
q Q2i U 2i jM 2 Θj;
(6d)
where Θ ∈ 0; π. Then we regard Θ as an independent variable and compute the maximum of I 0s with the numerical method over Θ. Ultimately, the resulting expression of the maximum of I 0s is written as I 0s;max maxfI i M 1 Θ
q Q2i U 2i jM 2 Θjg;
(6e)
where Θ ∈ 0; π. Now, based on Eq. (6e), we can use the rejection method to obtain scattering angles, while in the brute-force rejection method [40], it is based on the maximum of Eq. (6b) to determine the scattering angles. Compared with the time taken to find the maximum of Eq. (6d) over Θ ∈ 0; π, much more time is needed to get the maximum of Eq. (6b) over Θ ∈ 0; π and i1 ∈ 0; 2π. The detailed procedures of the rejection method for determining the scattering angles are presented as follows. (1) i1 and Θ are sampled from a uniform distribution Θ ∈ 0; π;
i1 ∈ 0; 2π:
(7)
(2) i1 and Θ are substituted into Eq. (6b) to obtain I 0s i1 ; Θ. Then the value of I 0s i1 ; Θ is compared with that of ξI 0s;max, where ξ is a uniform random number between 0 and 1. If I 0s i1 ; Θ ≥ ξI 0s;max, then i1 and Θ are accepted as the new scattering angles; otherwise, we reject i1 and Θ as the new scattering angles and go back to step 1 to resample them again.
(3) Once getting the new scattering angles i1 and Θ, we are able to determine the new propagation direction and Stokes vectors in the observer’s frame (the x–y–z frame in Fig. 2). The angles i2, θ2 , φ1 − φ2 , and φ2 can be derived from the spherical laws of sines and cosines. In the Stokes vector, the first component denotes the energy carried by photons, and the other three components store the states of polarization. So in the MC simulation, in order to ensure the conservation of energy, the normalization must be applied to the Stokes vector in the propagation of the photon; that is to say, the I-Stokes element of the photon is always equal to 1.0 as it propagates through. The normalization procedure should be carried out for every state-change event such as emitting, scattering, reflection, and transmission. B.
in which n nh ∕ns (nh is the refractive index with the higher value and ns is the one with the lower value). When going in the reverse direction, from a higher refractive index to a lower one, the ρ is given by ρn 1 −
1 1 − Fn: n2
(10b)
We compare a uniform random number Rref with the ρ to determine whether the reflection takes place or not at the interface/boundary. If Rref < ρ, the photon is reflected back to the incident layer; otherwise, it is refracted into another adjacent layer. The reflection direction is determined by the equations p 1 − Rθ ; θr arccos
φr 2πRφ ;
(11)
Reflection and Refraction
When the photon hits the semitransparent boundary surface or the interface between the two adjacent layers, the reflection or refraction will take place. According to the reflection characteristics, the medium interfaces or boundary mainly include Fresnel and Lambertian surfaces. The reflection Stokes vector and refraction Stokes vector are, respectively, written as Sref RSi ;
(8a)
Srefr TSi ;
(8b)
(12a)
where R and T are reflection matrix and refraction matrix, respectively. Note that Sref and Srefr also need to be normalized. For the Lambertian surface, the reflection matrix and refraction matrix are given as 0 1 0 1 ρ 0 0 0 1−ρ 0 0 0 B0 0 0 0C B 0 0 0 0C C C RB TB @ 0 0 0 0 A; @ 0 0 0 0 A; (9) 0 0 0 0 0 0 0 0 where ρ is the reflection coefficient of the interface/ surface which is independent to the incident beam and just related to the properties of the surface. For the semitransparent interface/surface, ρ is determined by the refractive indices of the adjacent two layers. Siegel [41] presented the expressions of ρ for semitransparent Lambertian surface. When the photon passes into a layer with higher refractive index, ρ is given by ρn ≡ Fn −
n2 n2 − 12
where Rθ and Rφ are uniform pseudorandom numbers. Both the incident angle and reflection angle are measured with respect to the normal direction of the reflection plane. The Fresnel reflection and refraction matrices for a flat dielectric interface have been presented in a few papers [22,24,42]. The reflection matrix is written as 0 2 1 0 0 rl r2r r2l − r2r C 1 B r2l − r2r r2l r2r 0 0 C Rθi B g 2 Imfr r g A; @ 0 0 2 Refr r 2 r l l r 0 0 2 Imfrr rl g 2 Refrr rl g
3n 1n − 1 n−1 ln 2 2 3 n1 6n 1 n 1
2n3 n2 2n − 1 8n4 n4 1 lnn 0.5; 2 2 4 n 1n − 1 n 1n4 − 12
in which
q n2t − n2i sin2 θi q ; rr ni cosθi n2t − n2i sin2 θi q n2t cosθi − ni n2t − n2i sin2 θi q ; rl n2t cosθi ni n2t − n2i sin2 θi ni cosθi −
(12b)
where r and l denote the components perpendicular and parallel to the plane of incidence, respectively; Refg and Imfg denote the real and imaginary parts, respectively; and ni and nt are the refractive indices of the media for the incident side and the refraction side, respectively. The refraction matrix is expressed as cosθt nt 3 Tθi 2 cosθi ni 0 t2 t2r t2l − t2r B l B t2 − t2 t2 t2 r r B l ×B l B 0 0 @ 0 0
0
0
0
0
1
C C C C; 2 Reftr tl g 2 Imftl tr g C A 2 Imftr tl g 2 Reftr tl g
(10a)
(12c) 1 March 2014 / Vol. 53, No. 7 / APPLIED OPTICS
1431
in which tr
2ni cosθi q ; ni cosθi n2t − n2i sin2 θi
tl
2ni nt cosθi q : n2t cosθi ni n2t − n2i sin2 θi
(12d)
The reflection and refraction events are controlled by the laws of reflection and retraction. The angles after reflection and refraction are given by the equations θr θi ;
nt sin θt ni sin θi ;
φr φt φi : (13)
We also need to compare the reflection coefficient ρ with a uniform pseudorandom number Rref to determine whether the reflection or refraction take place when the photon hits the Fresnel interface/boundary. For the Fresnel interface/boundary, the reflection coefficient is written as ρ C.
r2l r2r : 2
(14)
Photon Statistics and Post-processing
In order to get the real Stokes vector passing through some plane along some direction, we need to record the number of the photons passing through this plane along the θ; φ direction. In this work, we present the concept of VRTM to solve this problem. For a one-dimensional problem, the VRTM denoted by Di;θ;φ is defined as Di;θ;φ
Nt X Sp p1
N
;
(15)
in which N denotes the total number of emitting photons at the top surface of the multilayer system, as shown in Fig. 1; N t denotes the total number of the photons which pass through plane i from the θ; φ direction in the solid angle ΔΩ; Sp is the normalized Stokes vector of photon p. The physical meaning of VRTM is clear and it represents the transfer of intensity and polarized states of photons from the emission position to the receiving position. Finally, we obtain the real Stokes vector, written as Si;θ;φ
I 0 j cos θ0 j ·D ; j cos θjΔΩ i;θ;φ
(16)
in which I 0 is the intensity of incident beam and cos θ0 is the incident-direction cosine, as shown in Fig. 1. 3. Results and Discussion
In Section 2, we developed an arbitrary multilayer model for polarized radiative transfer. In Section 3, we first need to validate the correctness of the model. 1432
APPLIED OPTICS / Vol. 53, No. 7 / 1 March 2014
The MC model for polarized radiative transfer in the arbitrary multilayer system built in this work can be validated against the two-layer model for the atmosphere–ocean system. Then we will give the cases for the polarized radiative transfer in the three-layer, five-layer, and ten-layer system considering different optical boundary/interface conditions, refractive indices, optical thicknesses, and scattering albedos in different layers. Note that the Stokes vector is composed of direct light and diffuse light [20], and the results presented in this work are for the diffuse Stokes vector. According to Ref. [20], the direct Stokes vector is defined as Sdir S0 δcos θ − cos θ0 δφ − φ0 exp−τ∕ cos θ; (17a) where δ is the Dirac delta function; θ0 ; φ0 is the direction of the beam propagation without experiencing scattering. The diffuse Stokes vector is then defined as Sdif S − Sdir :
(17b)
It should be noted that the three layer, five-layer, and ten-layer cases are all fictitious problems, and do not have any physical significance. They are mainly used to test the computational implementation of the developed MC model for polarized radiative transfer in a multilayer (K > 2) system with complex properties. Case 1: A two-layer model of the atmosphere– ocean system In this case, we present the results compared against those obtained by the coupled atmosphere– water vector discrete ordinate radiative transfer (CAW-VDISORT) code in Ref. [29] to verify the MC polarized radiative transfer model with consideration of Fresnel interface developed in this work. The simulation is performed for polarized radiative transfer in a plane–parallel medium composed of Rayleigh scattering atmosphere layer and water body. A relative refractive index of water to the atmosphere is taken as 1.338. The optical depth is taken to be 0.15 for the atmosphere and 1.0 for the water. Both the atmosphere layer and water body are assumed to be nonabsorbing (ω 1.0), and the bottom of the water is assumed to be completely absorbing. A two-layer model is selected for the comparison. In the two-layer model, the top surface of the atmosphere and the bottom surface of the water are supposed to be Fresnel surfaces with n0 n1 and n3 n2 . The downward radiation field at the top of the atmosphere is taken to be a collimated and polarized beam of irradiance S0 1; 0; 0.5; 0.50 with an incident direction of θ0 ; φ0 60°; 0°. The number of the photons used in the MC simulation is 1.0 × 1010 and the angular discretization is taken as N θ × N φ 180 × 40 for this test. Figure 3 plots the curves for the four elements of the Stokes vector just below the water surface. Results for three different azimuth angles 0°, 90°,
1.4
×10
−1
4.0
×10
−2
ΜC ϕ = 0° ΜC ϕ = 90° ΜC ϕ = 180° Ref. [29] ϕ = 0° Ref. [29] ϕ = 90° Ref. [29] ϕ = 180°
3.0
1.2
2.0
1.0 ΜC ϕ = 0° ΜC ϕ = 90° ΜC ϕ = 180° Ref. [29] ϕ = 0° Ref. [29] ϕ = 90° Ref. [29] ϕ = 180°
0.6 0.4 0.2 -1.0 5.0
×10
-0.5 −2
0.0 cosθ
0.5
Q
I
1.0
0.8
0.0 -1.0 -2.0
1.0
-3.0 -1.0 5.0
4.0
4.0
3.0
3.0
2.0 0.0
0.5
1.0
0.5
1.0
−2
-1.0 -2.0
1.0 V
U
×10
0.0 cosθ
2.0 ΜC ϕ = 0° ΜC ϕ = 90° ΜC ϕ = 180° Ref. [29] ϕ = 0° Ref. [29] ϕ = 90° Ref. [29] ϕ = 180°
1.0
-1.0 -2.0 -3.0
-4.0
-4.0 -0.5
0.0 cosθ
0.5
ΜC ϕ = 0° ΜC ϕ = 90° ΜC ϕ = 180° Ref. [29] ϕ = 0° Ref. [29] ϕ = 90° Ref. [29] ϕ = 180°
0.0
-3.0 -5.0 -1.0
-0.5
1.0
-5.0 -1.0
-0.5
0.0 cosθ
Fig. 3. Stokes vector elements just below the atmosphere–water interface for a collimated and polarized incident beam.
and 180° are given for each element. As shown in Fig. 3, an excellent agreement between the results obtained by the MC model of this work and those by CAW-VDISORT [29] can be achieved. Case 2: A three-layer model with Fresnel surfaces/ interfaces In this case, we present a three-layer model with Fresnel surfaces/interfaces. The model is just for testing the computational implementation of our solution to polarized radiative transfer in a multilayer (K > 2) semitransparent medium with phase matrices more general than the pure Rayleigh scattering phase matrix considered in the first case. Note that the phase matrices used in layers 1 and 3 are presented in Ref. [20], including the cloud scattering and aerosol scattering. The Rayleigh scattering is considered in layer 2. The corresponding incident beam for these scattering matrices is at a wavelength of 0.412 μm. The refractive indices of the aerosol scattering, pure Rayleigh scattering, and cloud scattering layers are 1.385, 1, and 1.339, respectively. The refractive index background environment surrounding the three-layer system is equal to 1. Scattering albedos for layer 1, layer 2, and layer 3 are 0.3, 0.5, and 0.2, respectively; optical thickness for layer 1, layer 2, and layer 3 are 1, 0.3, and 0.5, respectively. The downward radiation field at the top of the multilayer medium (outside the top surface) is taken to be a collimated and polarized beam of irradiance S0 1; 0; 0.5; 0.50 with an incident direction of θ0 ; φ0 60°; 0°. The number of the photons used
in the model is 2.0 × 109 and the angular discretization is taken as N θ × N φ 90 × 20 for the simulation. Figures 4–6 plot the curves for the four Stokes vector elements at the interfaces of every layer. Results for two different azimuth angles 90° and 180° are given for each element. Resulting from the total reflection occurring at the interfaces of layers 1 and 3, Figs. 4 and 6 demonstrate the obvious peaks and valleys of abrupt change in the Stokes vector. At the two interfaces of layer 1, the peaks and valleys appear at the polar angle with a cosine of μ −0.67, and at the two interfaces of layer 3, they appear at μ 0.69. From Fig. 5, we can also find the peaks and valleys of the Stokes vector in the direction approximately parallel to the two interfaces of layer 2. Case 3: A three-layer model with Fresnel and Lambertian surface/interface In this case, we consider a three-layer model with Fresnel and Lambertian surface. The phase matrix for layer 1 and layer 3 is for 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 form of the phase matrix is as in Eq. (4a), and the Legendre coefficients for the four unique elements of the phase matrix have been obtained by Evens and Stephens [8]. Pure Rayleigh scattering is considered in layer 2, and the form of the phase matrix is given by Eq. (4b). In the three-layer system, the interface between layers 0 and 1 is a Lambertian one, and the other surfaces and interfaces are 1 March 2014 / Vol. 53, No. 7 / APPLIED OPTICS
1433
2.5
×10
−2
8.0
1.0
1-layer-top, ϕ = 90° 1-layer-bottom, ϕ = 90° 1-layer-top, ϕ = 180° 1-layer-bottom, ϕ = 180°
4.0 2.0 0.0
Q
I
1.5
−3
6.0
1-layer-top, ϕ = 90° 1-layer-bottom, ϕ = 90° 1-layer-top, ϕ = 180° 1-layer-bottom, ϕ = 180°
2.0
×10
-2.0
0.5
-4.0 0.0
-6.0
-0.5 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ 1.5
×10
-8.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
−2
0.4 1-layer-top, ϕ = 90° 1-layer-bottom, ϕ = 90° 1-layer-top, ϕ = 180° 1-layer-bottom, ϕ = 180°
1.0
−2
0.2 0.0 -0.2
1-layer-top, ϕ = 90° 1-layer-bottom, ϕ = 90° 1-layer-top, ϕ = 180° 1-layer-bottom, ϕ = 180°
V
U
0.5
×10
-0.4
0.0
-0.6 -0.5
-0.8
-1.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
-1.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
Fig. 4. Stokes vector elements at the top and bottom surfaces of layer 1 for the three-layer model.
3.5
×10
−2
1.0 2-layer-top, ϕ = 90° 2-layer-bottom, ϕ = 90° 2-layer-top, ϕ = 180° 2-layer-bottom, ϕ = 180°
3.0 2.5
×10
−2
2-layer-top, ϕ = 90° 2-layer-bottom, ϕ = 90° 2-layer-top, ϕ = 180° 2-layer-bottom, ϕ = 180°
0.8 0.6 0.4
I
Q
2.0 1.5
0.2
1.0
0.0
0.5
-0.2
0.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
-0.4 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
2.0 1.5 1.0
×10
−2
0.4 2-layer-top, ϕ = 90° 2-layer-bottom, ϕ = 90° 2-layer-top, ϕ = 180° 2-layer-bottom, ϕ = 180°
0.0 -0.2 -0.4 V
U -0.5
-0.6 -0.8 -1.0
-1.0 -1.5 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
−2
0.2
0.5 0.0
×10
-1.2
2-layer-top, ϕ = 90° 2-layer-bottom, ϕ = 90° 2-layer-top, ϕ = 180° 2-layer-bottom, ϕ = 180°
-1.4 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
Fig. 5. Stokes vector elements at the top and bottom surfaces of layer 2 for the three-layer model. 1434
APPLIED OPTICS / Vol. 53, No. 7 / 1 March 2014
3.5
×10
−2
1.4
2.5
0.8 0.6
1.5
0.4
1.0
0.2
0.5
0.0
0.0
-0.2
-0.5 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
1.0
U
0.5
×10
-0.4 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
−2
0.4
0.0 -0.5 -1.0 -1.5 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
×10
−2
0.2
3-layer-top, ϕ = 90° 3-layer-bottom, ϕ = 90° 3-layer-top, ϕ = 180° 3-layer-bottom, ϕ = 180°
0.0 -0.2 -0.4 V
1.5
3-layer-top, ϕ = 90° 3-layer-bottom, ϕ = 90° 3-layer-top, ϕ = 180° 3-layer-bottom, ϕ = 180°
1.0
Q
I
2.0
−2
1.2
3-layer-up, ϕ = 90° 3-layer-down, ϕ = 90° 3-layer-up, ϕ = 180° 3-layer-down, ϕ = 180°
3.0
×10
-0.6 -0.8 -1.0
3-layer-top, ϕ = 90° 3-layer-bottom, ϕ = 90° 3-layer-top, ϕ = 180° 3-layer-bottom, ϕ = 180°
-1.2 -1.4 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
Fig. 6. Stokes vector elements at the top and bottom surfaces of layer 3 for the three-layer model.
Fresnel ones. The refractive index of the background environment surrounding the multilayer system is equal to 1. The refractive indices of the pure Rayleigh scattering and Mie scattering layers are 1.338 and 1.44, respectively. Scattering albedos for layers 1, 2, and 3 are 1, 0.5, and 0.8, respectively; optical thickness for layers 1, 2, and 3 are 0.4, 0.3, and 0.3, respectively. The downward radiation field at the top of the three-layer system is taken to be a collimated and unpolarized beam of irradiance S0 1; 0; 0; 00 with an incident direction of θ0 ; φ0 60°; 0°. The number of the photons used for this case is 7.0 × 109 and the angular discretization is taken as N θ × N φ 90 × 20 for the simulation. Figures 7–9 show the curves for the four components of the Stokes vector at the top and bottom surfaces of every layer, and results for the azimuth angle φ 90° are presented. In Fig. 7, Q, U, and V components at the bottom surface of layer 1 for μ > 0 are also equal to zero due to the diffuse reflection on the Lambertian surface eliminating the polarization. The refractive index of layer 3 is larger than that of layer 2 and the background environment, so the total reflection will take place on the two Fresnel surfaces of layer 3 when the incident angle is greater than the critical angle. At the two surfaces of layer 3, the abrupt changes in the Stokes vector are observed at μ −0.72 and μ 0.37, as shown in Fig. 9. For layer 1, the total reflection takes place on the top surface, and the abrupt changes in the Stokes vector are found at μ −0.37 and μ −0.72; resulting from the
effect of total reflection occurring at the top surface of layer 1, similar changes in the Stokes vector can also observed at the bottom surface of layer 1, as shown in Fig. 7. Compared with the two adjacent layers, the refractive index of layer 2 is smaller, but compared with the background environments, the refractive index of layer 2 is bigger. So the Stokes vector at the surfaces of layer 2 is still influenced by the total reflection occurring at the surfaces of other layers. From Fig. 8, we can also see the abrupt changes in the Stokes vector at the two surfaces of layer 2 at μ −0.6644. Case 4: Five-layer and ten-layer models with Fresnel surfaces/interfaces In this case, we consider a five-layer model and a ten-layer model with Fresnel surfaces/interfaces. For the multilayer system, every layer has the same physical property parameters except for the refractive index. The total optical thickness of the multilayer system is 2.0, and the scattering albedo ω is equal to 0.5 for the multilayer system. The scattering type of the multilayer system is pure Rayleigh scattering. The refractive indices of the multilayer system are considered as a linear step distribution along the vertical direction. Two refractive index distributions are presented, and from the bottom layer to the top layer the distributions of refractive indices are two arithmetic progressions from 2.0 to 1.2 and from 1.2 to 2.0, respectively. The downward radiation field at the top of the multilayer system is taken to be a collimated and polarized beam of irradiance 1 March 2014 / Vol. 53, No. 7 / APPLIED OPTICS
1435
9.0 ×10
−2
3.0
×10
−3
2.0
8.0
1-layer-top,ϕ = 90° 1-layer-bottom,ϕ = 90°
1.0
7.0 Q
I
0.0
6.0
-1.0
5.0
-2.0
1-layer-top,ϕ = 90° 1-layer-bottom,ϕ = 90°
4.0
-3.0
3.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ 2.0
×10
-4.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
−2
4.0 3.0
1.0
2.0
−3
1-layer-top,ϕ = 90° 1-layer-bottom,ϕ = 90°
V
U
1.5
×10
0.5
1.0
0.0
0.0
1-layer-top,ϕ = 90° 1-layer-bottom,ϕ = 90°
-0.5 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
-1.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
Fig. 7. Stokes vector elements at the top and bottom surfaces of layer 1 for the three-layer model.
6.5
×10
−2
2.0
6.0
0.0
5.0
-1.0
4.5
Q
I
−3
1.0
5.5
-2.0
4.0
-3.0
3.5 2-layer-top,ϕ = 90° 2-layer-bottom,ϕ = 90°
3.0
×10
2-layer-top,ϕ = 90° 2-layer-bottom,ϕ = 90°
-4.0
2.5 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ 2.0
×10
-5.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
−2
6.0 5.0
1.5
×10
−3
2-layer-top,ϕ = 90° 2-layer-bottom,ϕ = 90°
4.0 3.0 V
U
1.0
2.0 0.5
0.0
1.0 2-layer-top,ϕ = 90° 2-layer-bottom,ϕ = 90°
-1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
0.0 -1.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
Fig. 8. Stokes vector elements at the top and bottom surfaces of layer 2 for the three-layer model. 1436
APPLIED OPTICS / Vol. 53, No. 7 / 1 March 2014
8.0
×10
−2
3.0
−3
2.0
7.0
3-layer-top,ϕ = 90° 3-layer-bottom,ϕ = 90°
1.0
6.0
0.0
5.0
-1.0
4.0
-2.0
Q
I
×10
-3.0
3.0
3-layer-top,ϕ = 90° 3-layer-bottom,ϕ = 90°
2.0
-4.0 -5.0
1.0
-6.0
0.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
-7.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
3.0
×10
−2
1.5
×10
−2
2.5 1.0
2.0
U
V
1.5 0.5
1.0 0.5 0.0
3-layer-top,ϕ = 90° 3-layer-bottom,ϕ = 90°
-0.5 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
0.0
3-layer-top,ϕ = 90° 3-layer-bottom,ϕ = 90°
-0.5 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
Fig. 9. Stokes vector elements at the top and bottom surfaces of layer 3 for the three-layer model.
S0 1; 0.5; 0.5; 0.50 with an incident direction of θ0 ; φ0 30°; 0°. The number of the photons used for these two cases is 5.0 × 109 and the angular discretization is taken as N θ × N φ 90 × 20 for the simulations. The results only for φ 90° are plotted. Figures 10 and 11 illustrate the curves for the four elements I, Q, U, and V of the Stokes vector at the bottom surface and the top surface of the five-layer model, respectively. It is well known that, as light beams enter into the medium with smaller refractive index from the medium with bigger refractive index, total reflection will happen at the interface only if the incident angle is not smaller than the critical angle of total reflection. It should also be noted that the scattering can change the direction of the light propagation. Based on these two points, we can give a qualitative analysis of the effects of total reflection occurring at the interfaces on the Stokes vector at the top and bottom surfaces of the multilayer system. In the multilayer system with a distribution of refractive indices increasing along the incident direction of external beams, Stokes vector at the top surface of the fifth layer is mainly affected by the total reflection at the top surface corresponding to the critical-angle cosine of cos θ ≈ −0.553, which gives rise to a big abrupt change in curves of Stokes vector just at the direction with cos θ ≈ −0.553, as shown in Fig. 11. In fact, the Stokes vector at the top surface can also be affected by the total reflection from the bottom surface, although the effect is so weak that it can be neglected. At the bottom surface,
the Stokes vector is affected by the total reflection from all six interfaces in the multilayer system. With respect to the bottom surface of the first layer with n1 2, corresponding to the critical angles for total reflection at the top surfaces of layers 1 to 4, the bottom surface of layer 1, and the top surface of layer 5, the cosines of the incident angles can be obtained, equal to 0.436, 0.6, 0.714, 0.8, and 0.866, respectively, where abrupt changes happen in the curves of the Stokes vector at the bottom surface of the multilayer system, and consequently, in the ranges of [0.436, 0.866] and −0.866; −0.436, fluctuations like waves can be seen, as shown in Fig. 10. In the multilayer system with a distribution of refractive indices decreasing along the incident direction of external beams, a big abrupt change in the Stokes vector at the bottom surface can be found at the direction with cos θ ≈ 0.553, as shown in Fig. 10, and abrupt changes in the Stokes vector at the top surface can be observed at the directions with 0.436, 0.6, 0.714, 0.8, and 0.866, respectively, as shown in Fig. 11. Figures 12 and 13 plot the curves of the Stokes vector at the top and bottom surfaces of the ten-layer system for the two distributions of refractive indices. From Figs. 10 to 13, we can see that the curves of the Stokes vector for the ten-layer model are similar to those for the five-layer model, and compared with the curves of the Stokes vector for the five-layer model, the fluctuations like waves in the curves of the Stokes vector for the ten-layer model are of 1 March 2014 / Vol. 53, No. 7 / APPLIED OPTICS
1437
−2
−2
4.5 ×10
0.0 ×10
4.0 3.5
nbottom=2.0, ntop=1.2, ϕ=90°
3.0
nbottom=1.2, ntop=2.0, ϕ=90°
-0.5 -1.0
2.5
Q
I
z/L=0
2.0
nbottom=2.0, ntop=1.2, ϕ=90°
-1.5
nbottom=1.2, ntop=2.0, ϕ=90°
1.5 1.0
-2.0
z/L=0
0.5 0.0 -1.0
-0.5 −2
×10
0.0 cosθ
0.5
-2.5 -1.0
1.0
-0.5 −2
1.4 ×10
0.5
0.5
1.0
1.2 1.0
0.0
0.8 -0.5
nbottom=2.0, ntop=1.2, ϕ=90°
V
U
0.0 cosθ
nbottom=1.2, ntop=2.0, ϕ=90°
nbottom=2.0, ntop=1.2, ϕ=90° nbottom=1.2, ntop=2.0, ϕ=90°
0.6 0.4
-1.0
z/L=0
0.2
z/L=0
0.0
-1.5 -1.0
-0.5
0.0 cosθ
0.5
-0.2 -1.0
1.0
-0.5
0.0 cosθ
0.5
1.0
Fig. 10. Stokes vector elements at the bottom interface for the five-layer model.
−2
−2
6.0 ×10
0.0 ×10
5.0
nbottom=2.0, ntop=1.2, ϕ=90°
-1.0
nbottom=1.2, ntop=2.0, ϕ=90°
4.0 z/L=1.0
3.0
nbottom=2.0, ntop=1.2, ϕ=90°
Q
I
-2.0 -3.0
nbottom=1.2, ntop=2.0, ϕ=90°
2.0
z/L=1.0
-4.0
1.0 0.0 -1.0
-0.5
0.0 cosθ
0.5
-5.0 -1.0
1.0
−2
-0.5
0.0 cosθ
0.5
1.0
−2
4.0 ×10
3.0 ×10
3.0
2.0
nbottom=2.0, ntop=1.2, ϕ=90° nbottom=1.2, ntop=2.0, ϕ=90°
1.0
z/L=1.0
V
U
2.0 nbottom=2.0, ntop=1.2, ϕ=90°
0.0
nbottom=1.2, ntop=2.0, ϕ=90°
1.0
-1.0
z/L=1.0
0.0 -1.0
-0.5
0.0 cosθ
0.5
1.0
-2.0 -1.0
-0.5
0.0 cosθ
Fig. 11. Stokes vector elements at the top interface for the five-layer model. 1438
APPLIED OPTICS / Vol. 53, No. 7 / 1 March 2014
0.5
1.0
−2
−2
5.0 ×10
0.0 ×10
4.0
-0.5 nbottom=2.0, ntop=1.2, ϕ=90° nbottom=1.2, ntop=2.0, ϕ=90°
-1.0 Q
I
3.0
z/L=0
2.0
nbottom=2.0, ntop=1.2, ϕ=90°
-1.5
nbottom=1.2, ntop=2.0, ϕ=90°
1.0
z/L=0
-2.0
0.0 -1.0
-0.5 −2
0.0 cosθ
0.5
-2.5 -1.0
1.0
-0.5 −2
1.0 ×10
0.0 cosθ
0.5
1.0
1.6 ×10 1.4
0.5
1.2 1.0
0.0 V
U
nbottom=2.0, ntop=1.2, ϕ=90° nbottom=1.2, ntop=2.0, ϕ=90°
-0.5
0.8
nbottom=2.0, ntop=1.2, ϕ=90°
0.6
nbottom=1.2, ntop=2.0, ϕ=90° z/L=0
0.4
z/L=0
0.2
-1.0
0.0 -1.5 -1.0
-0.5
0.0 cosθ
0.5
-0.2 -1.0
1.0
-0.5
0.0 cosθ
0.5
1.0
Fig. 12. Stokes vector elements at the bottom interface for the ten-layer model.
−2
−2
6.0 ×10
0.0 ×10
5.0
nbottom=2.0, ntop=1.2, ϕ=90° nbottom=1.2, ntop=2.0, ϕ=90°
-1.0
z/L=1.0
4.0 Q
I
-2.0 3.0 nbottom=2.0, ntop=1.2, ϕ=90°
2.0
-3.0
nbottom=1.2, ntop=2.0, ϕ=90°
-4.0
z/L=1.0
1.0 0.0 -1.0
-0.5
0.0
0.5
-5.0 -1.0
1.0
cosθ
−2
×10
-0.5
0.0
0.5
1.0
cosθ
−2
3.0 ×10
4.0 nbottom=2.0, ntop=1.2, ϕ=90°
2.0
nbottom=1.2, ntop=2.0, ϕ=90°
3.0
z/L=1.0
1.0 V
U
2.0 nbottom=2.0, ntop=1.2, ϕ=90°
1.0 0.0 -1.0 -1.0
0.0
nbottom=1.2, ntop=2.0, ϕ=90°
-1.0
z/L=1.0
-0.5
0.0 cosθ
0.5
1.0
-2.0 -1.0
-0.5
0.0 cosθ
0.5
1.0
Fig. 13. Stokes vector elements at the top interface for the ten-layer model. 1 March 2014 / Vol. 53, No. 7 / APPLIED OPTICS
1439
higher frequency and smaller amplitudes. Such calculations took 6 h 2 min and 7 h 34 min of CPU time for the five-layer and ten-layer cases, respectively. The program was running at 3.0 GHz on a single core of an Intel Core i5 processor. The computation time increases significantly with the number of layers increasing, since reflection and transmission events increase many times. The ten-layer model is closer to a continuum, while it needs more computation time to obtain converged results. It is expected that the approximate solution to vector radiative transfer will be implemented in a gradient-index medium by dividing the medium into a sufficient number of sublayers. Based on the computation time spent by the ten-layer case, more or unacceptable time may be needed for the gradient-index problem. Optimizing the code to reduce the computation time will be our next work. 4. Conclusions
In this work, we developed a MC model for solving arbitrary multilayer vector radiative transfer problems with Fresnel and Lambertian surfaces/ interfaces. We considered two kinds of polarization mechanisms, scattering by particles and reflection/ refraction by Fresnel surfaces/interfaces. We presented the detailed steps of determining scattering angles by using the modified rejection method. We put forward the concept of VRTM for describing the polarized radiative transfer. Using the VRTM solved by the MC method, we obtained four components of the Stokes vector. By comparing against the results of the two-layer model for the atmosphere–ocean system in the literature, we verified our multilayer MC method model with high accuracy and flexible adaptability. By comparison, we can see that all of the four components of the Stokes vector obtained by MC simulation are in good agreement with the results by DOM reported in the published literature. Next, we presented the results for two three-layer models with different optical boundary conditions and scattering phase matrices. Fresnel and Lambertian surfaces/interfaces were considered, and Rayleigh scattering and different types of Mie scattering from the literature were chosen for our models. Finally, two numerical cases of five-layer and ten-layer models were also presented at the top and bottom surfaces. By analyzing the results, we find the optical boundary/interface conditions have a significant effect on the polarized radiative transfer. In the multilayer system, due to inhomogeneous distribution of refractive index, the complex total reflection occurs at the surfaces/interfaces when the incident angle is equal to or greater than the critical angle, which generates abrupt changes in Stokes vector. Resulting from the total reflection at the surfaces/interfaces, fluctuations like waves can be observed in the curves of the Stokes vector at the surfaces/interfaces. With the increasing of the number of layers and the critical angles for total reflection, the frequency of fluctuations increases 1440
APPLIED OPTICS / Vol. 53, No. 7 / 1 March 2014
and the amplitudes of fluctuations decrease. The novelty of our work is that the arbitrary multilayer model developed can be applied to simulating vector radiative transfer in the inhomogeneous media which may be continuum or noncontinuum. For instance, it may be used to implement the approximate solution to vector radiative transfer in a gradientindex medium. In this work, we have made some attempts with the five-layer and ten-layer models to try to solve this problem. However, it is clear that we do not obtain smooth and stable solutions to polarized radiative transfer in a gradient-index medium with such a small number of layers. This work is supported by the National Natural Science Foundation of China (Grant Nos. 51176040 and 51276050). References 1. L. Tsang and J. A. Kong, “Radiative transfer theory for active remote sensing of half-space random media,” Radio Sci. 13, 763–773 (1978). 2. L. Tsang, J. A. Kong, and R. T. Shin, Theory of Microwave Remote Sensing (Wiley, 1985). 3. L. V. Wang and H. Wu, Biomedical Optics: Principles and Imaging (Wiley, 2007). 4. N. Ghosh, M. F. G. Wood, and I. A. Vitkin, Handbook of Photonics for Biomedical Science, V. V. Tuchin, ed. (Taylor and Francis, 2010), Chap. 9, pp. 253–282. 5. M. P. Mengüç and S. Manickavasagam, “Characterization of size and structure of agglomerates and inhomogeneous particles via polarized light,” Int. J. Eng. Sci. 36, 1569–1593 (1998). 6. G. W. Kattawar and G. N. Plass, “Radiance and polarization of multiple scattered light from haze and clouds,” Appl. Opt. 7, 1519–1527 (1968). 7. R. D. M. Garcia and C. E. Siewert, “The FN method for radiative transfer models that include polarization effects,” J. Quant. Spectrosc. Radiat. Transfer 41, 117–145 (1989). 8. K. F. Evans and G. L. Stephens, “A new polarized atmospheric radiative transfer model,” J. Quant. Spectrosc. Radiat. Transfer 46, 413–423 (1991). 9. C. E. Siewert, “A discrete-ordinates solution for radiativetransfer models that include polarization effects,” J. Quant. Spectrosc. Radiat. Transfer 64, 227–254 (2000). 10. J. Lenoble, M. Herman, J. L. Deuzé, B. Lafrance, R. Santer, and D. Tanré, “A successive order of scattering code for solving the vector equation of transfer in the earth’s atmosphere with aerosols,” J. Quant. Spectrosc. Radiat. Transfer 107, 479–507 (2007). 11. F. Xu, R. A. West, and A. B. Davis, “A hybrid method for modeling polarized radiative transfer in a spherical-shell planetary atmosphere,” J. Quant. Spectrosc. Radiat. Transfer 117, 59–70 (2013). 12. H. Ishimoto and K. Masuda, “A Monte Carlo approach for the calculation of polarized light: application to an incident narrow beam,” J. Quant. Spectrosc. Radiat. Transfer 72, 467–483 (2002). 13. X. D. Wang and L. V. Wang, “Propagation of polarized light in birefringent turbid media: a Monte Carlo study,” J. Biomed. Opt. 7, 279–290 (2002). 14. M. Xu, “Electric field Monte Carlo simulation of polarized light propagation in turbid media,” Opt. Express 12, 6530–6539 (2004). 15. R. Vaillon, B. T. Wong, and M. P. Mengüç, “Polarized radiative transfer in a particle-laden semi-transparent medium via a vector Monte Carlo method,” J. Quant. Spectrosc. Radiat. Transfer 84, 383–394 (2004). 16. C. Davis, C. Emde, and R. Harwood, “A 3-D polarized reversed Monte Carlo radiative transfer model for millimeter and submillimeter passive remote sensing in cloudy atmospheres,”
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
IEEE Trans. Geosci. Remote Sens. 43, 1096–1101 (2005). J. C. Ramella-Roman, S. A. Prahl, and S. L. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part I,” Opt. Express 13, 4420–4438 (2005). J. N. Swamy, C. Crofcheck, and M. P. Mengüç, “A Monte Carlo ray tracing study of polarized light propagation in liquid foams,” J. Quant. Spectrosc. Radiat. Transfer 104, 277–287 (2007). C. Cornet, L. C. Labonnote, and F. Szczap, “Three-dimensional polarized Monte Carlo atmospheric radiative transfer model (3DMCPOL): 3D effects on polarized visible reflectances of a cirrus cloud,” J. Quant. Spectrosc. Radiat. Transfer 111, 174–186 (2010). A. A. Kokhanovsky, V. P. Budak, C. Cornet, M. Z. Duan, C. Emde, I. L. Katsev, D. A. Klyukov, S. V. Korkin, L. C. Labonnote, B. Mayer, Q. L. Min, T. Nakajima, Y. Ota, A. S. Prikhach, V. V. Rozanov, T. Yokota, and E. P. Zege, “Benchmark results in vector atmospheric radiative transfer,” J. Quant. Spectrosc. Radiat. Transfer 111, 1931–1946 (2010). G. W. Kattawar, G. N. Plass, and J. A. Guinn, “Monte Carlo calculations of the polarization of radiation in the Earth’s atmosphere-ocean system,” J. Phys. Oceanogr. 3, 353–372 (1973). G. W. Kattawar and C. N. Adams, “Stokes vector calculations of the submarine light field in an atmosphere–ocean with scattering according to a Rayleigh phase matrix: effect of interface refractive index on radiance and polarization,” Limnol. Oceanogr. 34, 1453–1472 (1989). P. W. Zhai, Y. X. Hu, C. R. Trepte, and P. L. Lucker, “A vector radiative transfer model for coupled atmosphere and ocean systems based on successive order of scattering method,” Opt. Express 17, 2057–2079 (2009). P. W. Zhai, Y. X. Hu, J. Cdhary, C. R. Trepte, P. L. Lucker, and D. B. Josset, “A vector radiative transfer model for coupled atmosphere and ocean systems with a rough interface,” J. Quant. Spectrosc. Radiat. Transfer 111, 1025–1040 (2010). P. W. Zhai, Y. X. Hu, D. B. Josset, C. R. Trepte, P. L. Lucker, and B. Lin, “Advanced angular interpolation in the vector radiative transfer for coupled atmosphere and ocean systems,” J. Quant. Spectrosc. Radiat. Transfer 115, 19–27 (2013). X. Q. He, D. L. Pan, Y. Bai, Q. K. Zhu, and F. Gong, “Vector radiative transfer numerical model of coupled oceanatmosphere system using matrix-operator method,” Sci. China Ser. D 50, 442–452 (2007). X. Q. He, Y. Bai, Q. Zhu, and F. Gong, “A vector radiative transfer model of coupled ocean-atmosphere system using
28.
29.
30.
31. 32. 33.
34. 35. 36. 37. 38.
39. 40. 41. 42.
matrix-operator method for rough sea-surface,” J. Quant. Spectrosc. Radiat. Transfer 111, 1426–1448 (2010). Y. Ota, A. Higurashi, T. Nakajima, and T. Yokota, “Matrix formulations of radiative transfer including the polarization effect in a coupled atmosphere-ocean system,” J. Quant. Spectrosc. Radiat. Transfer 111, 878–894 (2010). E. R. Sommersten, J. K. Lotsberg, K. Stamnes, and J. J. Stamnes, “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. Transfer 111, 616–633 (2010). R. D. M. Garcia, “Fresnel boundary and interface conditions for polarized radiative transfer in a multilayer medium,” J. Quant. Spectrosc. Radiat. Transfer 113, 306–317 (2012). P. W. Zhai, G. W. Kattawar, and Y. X. Hu, “Comment on transmission matrix for dielectric interface,” J. Quant. Spectrosc. Radiat. Transfer 113, 1981–1984 (2012). R. D. M. Garcia, “Response to ‘Comment on transmission matrix for a dielectric interface’,” J. Quant. Spectrosc. Radiat. Transfer 113, 2251–2254 (2012). R. D. M. Garcia, “Radiative transfer with polarization in a multi-layer medium subject to Fresnel boundary and interface conditions,” J. Quant. Spectrosc. Radiat. Transfer 115, 28–45 (2013). S. Chandrasekhar, Radiative Transfer (Oxford University, 1950). H. C. Van de Hulst, Light Scattering by Small Particles (Dover, 1981). M. I. Mishchenko, Scattering, Absorption, and Emission of Light by Small Particles (NASA, 2002). A. T. Young, “Revised depolarization corrections for atmospheric extinction,” Appl. Opt. 19, 1333 (1980). H. H. Tynes, G. W. Kattawar, E. P. Zege, I. L. Katsev, A. S. Prikhach, and L. I. Chaikovskaya, “Monte Carlo and multicomponent approximation methods for vector radiative transfer by use of effective Mueller matrix calculations,” Appl. Opt. 40, 400–412 (2001). R. M. Green, Spherical Astronomy (Cambridge University, 1985). B. A. Whitney, “Monte Carlo radiative transfer,” arXiv: 1104.4990 (2011). R. Siegel and C. M. Spuckler, “Refractive index effects on radiation in an absorbing, emitting, and scattering laminated layer,” J. Heat Transfer 115, 194–200 (1993). K. Masuda and T. Takashima, “Computational accuracy of radiation emerging from the ocean surface in the model atmosphere–ocean system,” Pap. Met. Geophys. 37, 1–13 (1986).
1 March 2014 / Vol. 53, No. 7 / APPLIED OPTICS
1441
School of Energy Science and Engineering, Harbin Institute of Technology, 92 West Dazhi Street, Harbin 150001, China 2
e-mail: [email protected] 3
e-mail: [email protected]
Received 27 September 2013; revised 16 January 2014; accepted 20 January 2014; posted 31 January 2014 (Doc. ID 198456); published 27 February 2014
Polarized radiative transfer in a multilayer system is an important problem and has wide applications in various fields. In this work, a Monte Carlo (MC) model is developed to simulate polarized radiative transfer in a semitransparent arbitrary multilayer medium with different refractive indices in each layer. Two kinds of polarization mechanisms are considered: scattering by particles and reflection and refraction at the Fresnel surfaces or interfaces. The MC method has an obvious superiority in that complex mathematical derivations can be avoided in solving the polarization by Fresnel reflection and refraction in an arbitrary multilayer system. We define the vector radiative transfer matrix (VRTM), which describes the polarization characteristics of radiative transfer, and obtain four elements of Stokes vector using the VRTM. The results for the two-layer model by MC method are compared against those for coupled atmosphere–ocean model by the discrete–ordinate method available in the literature, which validates the correctness of the MC multilayer model of polarized radiative transfer. Finally, the results for three-layer, five-layer, and ten-layer models are presented in graphical form. Results show that in the multilayer system, total reflections occurring at the surfaces/interfaces have significant effects on the polarized radiative transfer, which causes abrupt changes or fluctuations like waves in the curves of the Stokes vector. © 2014 Optical Society of America OCIS codes: (030.5620) Radiative transfer; (230.4170) Multilayers; (290.4210) Multiple scattering; (290.7050) Turbid media. http://dx.doi.org/10.1364/AO.53.001427
1. Introduction
Polarized radiative transfer has important applications in atmospheric and ocean optics, remote sensing [1,2], biomedical optics [3,4], particle characterization [5], and many other engineering fields and research disciplines. Over more than 40 years, many numerical methods for vector radiative transfer have been proposed and developed, including the Monte Carlo (MC) method [6], F N method [7], adding– doubling method [8], discrete–ordinate (DO) method [9], successive order scattering (SOS) method [10], 1559-128X/14/071427-15$15.00/0 © 2014 Optical Society of America
Markov chain method [11], and so on. Among these methods, the MC method is very popular for polarized radiative transfer, especially in the 21st century [12–19]. The MC method is time consuming, and in order to overcome this disadvantage, some other numerical methods have been developed for solving the polarized radiative transfer. Garcia and Siewert [7] used the F N method to solve the complete and general polarization problem for a plane–parallel slab, and tabulated the results for the problem of L 13 and L 60. Evans and Stephens [8] built a polarized plane–parallel radiative transfer model for randomly oriented particles of arbitrary shape using a direct method combined with the doubling and adding method. Siewert [9] developed a DO 1 March 2014 / Vol. 53, No. 7 / APPLIED OPTICS
1427
method for polarized radiative transfer in a one-dimensional atmosphere of randomly orientated oblate spheroids. Lenoble et al. [10] presented a SOS code to solve the equation of vector radiative transfer in a plane–parallel earth atmosphere, and considered polarizing effects caused by aerosol scattering and reflectance by water or terrestrial surfaces. Xu et al. [11] developed the Markov chain method to simulate polarized radiative transfer in a pseudospherical atmosphere with solar illumination. Kokhanovsky et al. [20] 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, and presented the benchmark data for one-dimensional vector atmospheric radiative transfer. Most of the work focused on the polarized radiative transfer in scattering media with uniform refractive index. However, vector radiative transfer in a medium with spatial changes in the refractive index is a topic of importance in many areas of investigation. An atmosphere–ocean two-layer model of polarized radiative transfer is a classical and important problem. Kattawar and co-workers [21,22] first examined the problem using the MC method. They calculated the degree and direction of polarization, the ellipticity, and the radiance of the radiation at various levels in the atmosphere–ocean system [21], and obtained the complete four-component Stokes vector at any region in a fully inhomogeneous atmosphere– ocean system with inclusion of a stochastic interface [22]. After that, especially in the past five years, this problem has attracted the attention of many scholars. Zhai and co-workers [23,24] developed the SOS method for the problem. In 2013, Zhai et al. [25] applied the Iteration of Source Matrix method and single scattering method to the SOS method to solving the problem, and obtained the polarized radiation field at arbitrary viewing angles. He et al. [26,27] used the matrix–operator method for vector radiative transfer in the atmosphere–ocean system. Ota et al. [28] developed a solution, based on the DO and matrix–operator methods, to the problem of vector radiative transfer in the coupled atmosphere– ocean system. Sommersten et al. [29] presented an algorithm based on DO method for vector radiative transfer in a coupled system of atmosphere and ocean with different refractive indices, and compared the results by the DO method with those by the MC method. In the past two years, Garcia [30] investigated the vector radiative transfer in a multilayer medium with Fresnel boundaries or interfaces using DO method. His derivation of the reflection and transmission matrices gave rise to some arguments [31], and a consensus on this issue has not been reached [32]. In 2013, Garcia [33] applied the analytical discrete ordinates (ADO) method to calculate the polarized radiative transfer in a multilayer medium with Fresnel boundaries and interfaces, and tabulated the results of four Stokes elements for a nine-layer system. Up to now, there has been 1428
APPLIED OPTICS / Vol. 53, No. 7 / 1 March 2014
little work reported on polarized radiative transfer in a multilayer system with two or more interfaces considering different refractive indices in each layer. Although the MC method is time consuming, it has an obvious superiority over other numerical methods in that complex mathematical derivations can be avoided for solving the polarization problem caused by Fresnel reflection and refraction in an arbitrary multilayer system. In this work, we develop a MC model of vector radiative transfer in an arbitrary multilayer medium with consideration of Fresnel or Lambertian surfaces/interfaces. In different layers, not only can scattering type be different but physical parameters such as refractive index, scattering albedo, extinction coefficient, and so on can also be set to different values. In Section 2, we discuss the basic theory of vector radiative transfer, mainly including the descriptions for scattering, reflection, and transmission of photons. We will focus on how to determine scattering directions using the improved reject method in the MC simulation and the definition of vector radiative transfer matrix (VRTM). In Section 3, since most numerical results available in the literature are for two-layer models of the atmosphere–ocean system, we compare our results with those for atmosphere–ocean systems from the literature to verify our MC method multilayer model. Finally, we investigate polarized radiative transfer in three-layer, five-layer, and ten-layer models, respectively, and analyze the effects of total reflection at the Fresnel surfaces/interfaces on the results of Stokes vector. 2. Descriptions of the Problem
We consider the multilayer medium with an oblique incidence of a beam of parallel light on the top surface, as shown in Fig. 1. The medium includes K layers, and every layer has its own independent physical properties. In Fig. 1, nk is the refractive index of the k layer, and τk is the optical thickness of this layer; the refractive indices outside the top and bottom surfaces of the multilayer medium are S0 nK+1
θ0
nK
Layer K
τΚ
nk
Layer k
τκ
n2
Layer 2
τ2
n1
Layer 1
τ1
n0
τ
Fig. 1. Geometry of the multilayer problem.
denoted by nk1 and n0 , respectively. There are two surfaces for each layer except for layers nk1 and n0 , called top surface and bottom surface. For layer 0, there is only one top surface, and for layer K 1, there is only one bottom surface. Two boundary surfaces and each interface between the adjacent two layers are semitransparent. Two types of the interfaces/boundaries are considered in this model, and each interface or boundary can be Lambertian one or Fresnel one. The Stokes vector S I; Q; U; V0 is usually used to describe the intensity and polarized states of light. The specific form of Stokes vector is given by 0 1 1 El El Er Er I B Q C 1 B El E − Er Er C l B C C SB @ U A η @ El Er Er E A; l iEl Er − Er El V 0
(1)
where η μm ∕ε1∕2 ; μm is the magnetic permeability, and ε is the electric permittivity of the medium; El and Er are the complex amplitudes of the electric field along the two vectors which are parallel and perpendicular to a reference plane; “∗”denotes the complex conjugate. In this work, we use the MC technique to solve the polarized radiative transfer in the arbitrary multilayer medium as shown in Fig. 1. Note that we consider the cold medium in the present work and ignore the thermal emission of the medium. A.
Scattering
The photons are tracked in the multilayer medium, which experiences free propagation, scattering, reflection, transmission, and other events until the photons escape from the boundaries or are absorbed. The deterministic behavior of light is determined by the random behavior of a finite number of photons in the medium, and each photon randomly propagates along the possible path. The distance a photon travels without being absorbed, scattered, reflected, or refracted is determined by L−
ln1 − RL ; κ
photon will be absorbed by the medium and not tracked again; otherwise the photon will be scattered by the particle. The scattering problem for polarized radiative transfer is much more complicated than that for the scalar one. Chandrasekhar [34] first put forward the idea of photon scattering from one location to the next in the meridian and scattering planes. Figure 2 shows the geometry of the scattering problem. The directions of the photon propagation before and after scattering are represented as those pointed to by the lines OA and OB, respectively, on the unit sphere. The incident Stokes vector Si and the resulting scattering Stokes vector Ss determine a plane AOB, defined as the scattering plane. The vector of incident Stokes vector Si and the z axis determine a plane AOC, called the meridian plane. When a scattering event occurs, a new meridian plane is created by the new direction of propagation Ss (scattering Stokes vector) and the z axis, and the Stokes vector must be transformed so that the polarization is properly represented with respect to the new meridian plane BOC. Two rotations of the Stokes vector, in and out of the scattering plane, are required in the MC method. Based on the above analysis, the resulting Stokes vector after scattering is Ss Lπ − i2 MΘL−i1 Si ;
(3)
where the matrix M is the single scattering Mueller matrix defined in Van de Hulst [35], and L is a rotation matrix; i1 denotes the angle of rotation into the scattering plane AOB from the meridian plane AOC, and i2 denotes the angle of radiation into the meridian plane BOC from the scattering plane AOB. The term of Lπ − i2 MΘL−i1 is called the phase matrix. The matrix M has 4 × 4 elements 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 one that has eight nonzero elements [36], of which only four are independent. The form is
(2)
where L is the free distance, κ is the extinction coefficient, and RL is a pseudorandom number with a uniform distribution. The following MC analysis of multiple scattering events of the photons in a semitransparent medium is based on the assumption that the scattering event is independent and does not consider coherence effects. When the photons propagate in each layer of the medium, they may probably interact with the particles in the medium. The photons may be absorbed by the medium or scattered to other directions by the particles. A pseudorandom number is compared with the scattering albedo ω of the medium to determine whether the photon is absorbed or scattered. If the pseudorandom number Rω > ω the
z
Ss
C
Si
i2 B
θ2
Θ
θ1
i1 A
O
y
ϕ2 ϕ1 ϕ2
ϕ1
x
Fig. 2. Geometry of scattering plane and meridian planes. The photon’s direction of propagation before and after scattering is θ1 ; φ1 and θ2 ; φ2 , respectively. 1 March 2014 / Vol. 53, No. 7 / APPLIED OPTICS
1429
0
M1 B M2 B MΘ @ 0 0
M2 M1 0 0
0 0 M3 −M 4
1 0 0 C C; M4 A M3
(4a)
where Θ is the scattering angle between the incident and scattering directions. The simplified Mueller matrix Eq. (4a) can be used for Mie scattering and pure Rayleigh scattering as well. For Rayleigh scattering, the matrix M can also be written by 0
cos2 Θ 1 cos2 Θ − 1
0
2 2 0 3Δ B B cos Θ − 1 cos Θ 1 MΘ B 4 @ 0 0 2 cos Θ
0 0 01 0 0 01
0
1
0
C C C A
0 0 2Δ0 cos Θ
B0 0 0 0C C 1 − ΔB @ 0 0 0 0 A;
(4b)
0 0 0 0 in which Δ 1 − δ∕1 δ∕2 and Δ0 1 − 2δ∕ 1 − δ, where δ is the depolarization factor [37] which is equal to zero for pure Rayleigh scattering. Equation (4b) can be summarized as 0 1 0 M1 M2 0 B M2 M3 0 0 C C: MΘ B (4c) @ 0 0 M4 0 A 0 0 0 M5 The matrix L represents a rotation in the clockwise direction with respect to an observer looking into the direction of the photon propagation [38], and can be written as 0 1 1 0 0 0 B 0 cos 2ϕ sin 2ϕ 0 C C Lϕ B (5) @ 0 − sin 2ϕ cos 2ϕ 0 A: 0 0 0 1 The physical meaning of Eq. (3) is easy to understand. First of all, the incident Stokes vector Si is rotated at an i1 angle from the meridian plane AOC to the scattering plane AOB counterclockwise. Then, the incident Stokes vector Si interacts with the particle and is scattered in a new direction. Finally, the scattering Stokes vector is rotated at a π − i2 angle from the scattering plane AOB to the meridian plane BOC clockwise. Eventually, the scattering Stokes vector Ss is obtained. For a given incident direction θ1 ; φ1 as shown in Fig. 2, if the angles i1 and Θ are determined, then the angles i2, θ2 , and φ1 − φ2 can be calculated by the spherical laws of sines and cosines [39]. In order to obtain these angles, the cumulative distribution method [19] was used. Whitney [40] presented a rejection method to obtain i1 and Θ by using a brute-force numerical method in which all angles are calculated and compared with each other. The 1430
APPLIED OPTICS / Vol. 53, No. 7 / 1 March 2014
brute-force numerical method is time consuming to find i1 and Θ. In this work, we develop an optimized rejection method used for obtaining i1 and Θ, which obviously saves the computing time. The detailed steps are given as follows. Before we rotate the scattering Stokes vector from the scattering plane to the meridian plane BOC, Eq. (3) can be rewritten as S0s MΘL−i1 Si :
(6a)
Then the I Stokes element can be solved from Eq. (4a) (both for Mie and Rayleigh scattering), giving I 0s I i M 1 Θ Qi M 2 Θ cos 2i1 − U i M 2 Θ sin 2i1 ; (6b) where Θ ∈ 0; π and i1 ∈ 0; 2π. By using a trigonometric function formula, Eq. (6b) can be rewritten as I 0s I i M 1 Θ
q Q2i U 2i M 2 Θ cos2i1 ϕ; (6c)
where ϕ is an angle we are not interested in. Regarding i1 as independent variable and considering the maximum absolute value of the cosine function, we can obtain the maximum of I 0s;Θ, I 0s;Θ;max I i M 1 Θ
q Q2i U 2i jM 2 Θj;
(6d)
where Θ ∈ 0; π. Then we regard Θ as an independent variable and compute the maximum of I 0s with the numerical method over Θ. Ultimately, the resulting expression of the maximum of I 0s is written as I 0s;max maxfI i M 1 Θ
q Q2i U 2i jM 2 Θjg;
(6e)
where Θ ∈ 0; π. Now, based on Eq. (6e), we can use the rejection method to obtain scattering angles, while in the brute-force rejection method [40], it is based on the maximum of Eq. (6b) to determine the scattering angles. Compared with the time taken to find the maximum of Eq. (6d) over Θ ∈ 0; π, much more time is needed to get the maximum of Eq. (6b) over Θ ∈ 0; π and i1 ∈ 0; 2π. The detailed procedures of the rejection method for determining the scattering angles are presented as follows. (1) i1 and Θ are sampled from a uniform distribution Θ ∈ 0; π;
i1 ∈ 0; 2π:
(7)
(2) i1 and Θ are substituted into Eq. (6b) to obtain I 0s i1 ; Θ. Then the value of I 0s i1 ; Θ is compared with that of ξI 0s;max, where ξ is a uniform random number between 0 and 1. If I 0s i1 ; Θ ≥ ξI 0s;max, then i1 and Θ are accepted as the new scattering angles; otherwise, we reject i1 and Θ as the new scattering angles and go back to step 1 to resample them again.
(3) Once getting the new scattering angles i1 and Θ, we are able to determine the new propagation direction and Stokes vectors in the observer’s frame (the x–y–z frame in Fig. 2). The angles i2, θ2 , φ1 − φ2 , and φ2 can be derived from the spherical laws of sines and cosines. In the Stokes vector, the first component denotes the energy carried by photons, and the other three components store the states of polarization. So in the MC simulation, in order to ensure the conservation of energy, the normalization must be applied to the Stokes vector in the propagation of the photon; that is to say, the I-Stokes element of the photon is always equal to 1.0 as it propagates through. The normalization procedure should be carried out for every state-change event such as emitting, scattering, reflection, and transmission. B.
in which n nh ∕ns (nh is the refractive index with the higher value and ns is the one with the lower value). When going in the reverse direction, from a higher refractive index to a lower one, the ρ is given by ρn 1 −
1 1 − Fn: n2
(10b)
We compare a uniform random number Rref with the ρ to determine whether the reflection takes place or not at the interface/boundary. If Rref < ρ, the photon is reflected back to the incident layer; otherwise, it is refracted into another adjacent layer. The reflection direction is determined by the equations p 1 − Rθ ; θr arccos
φr 2πRφ ;
(11)
Reflection and Refraction
When the photon hits the semitransparent boundary surface or the interface between the two adjacent layers, the reflection or refraction will take place. According to the reflection characteristics, the medium interfaces or boundary mainly include Fresnel and Lambertian surfaces. The reflection Stokes vector and refraction Stokes vector are, respectively, written as Sref RSi ;
(8a)
Srefr TSi ;
(8b)
(12a)
where R and T are reflection matrix and refraction matrix, respectively. Note that Sref and Srefr also need to be normalized. For the Lambertian surface, the reflection matrix and refraction matrix are given as 0 1 0 1 ρ 0 0 0 1−ρ 0 0 0 B0 0 0 0C B 0 0 0 0C C C RB TB @ 0 0 0 0 A; @ 0 0 0 0 A; (9) 0 0 0 0 0 0 0 0 where ρ is the reflection coefficient of the interface/ surface which is independent to the incident beam and just related to the properties of the surface. For the semitransparent interface/surface, ρ is determined by the refractive indices of the adjacent two layers. Siegel [41] presented the expressions of ρ for semitransparent Lambertian surface. When the photon passes into a layer with higher refractive index, ρ is given by ρn ≡ Fn −
n2 n2 − 12
where Rθ and Rφ are uniform pseudorandom numbers. Both the incident angle and reflection angle are measured with respect to the normal direction of the reflection plane. The Fresnel reflection and refraction matrices for a flat dielectric interface have been presented in a few papers [22,24,42]. The reflection matrix is written as 0 2 1 0 0 rl r2r r2l − r2r C 1 B r2l − r2r r2l r2r 0 0 C Rθi B g 2 Imfr r g A; @ 0 0 2 Refr r 2 r l l r 0 0 2 Imfrr rl g 2 Refrr rl g
3n 1n − 1 n−1 ln 2 2 3 n1 6n 1 n 1
2n3 n2 2n − 1 8n4 n4 1 lnn 0.5; 2 2 4 n 1n − 1 n 1n4 − 12
in which
q n2t − n2i sin2 θi q ; rr ni cosθi n2t − n2i sin2 θi q n2t cosθi − ni n2t − n2i sin2 θi q ; rl n2t cosθi ni n2t − n2i sin2 θi ni cosθi −
(12b)
where r and l denote the components perpendicular and parallel to the plane of incidence, respectively; Refg and Imfg denote the real and imaginary parts, respectively; and ni and nt are the refractive indices of the media for the incident side and the refraction side, respectively. The refraction matrix is expressed as cosθt nt 3 Tθi 2 cosθi ni 0 t2 t2r t2l − t2r B l B t2 − t2 t2 t2 r r B l ×B l B 0 0 @ 0 0
0
0
0
0
1
C C C C; 2 Reftr tl g 2 Imftl tr g C A 2 Imftr tl g 2 Reftr tl g
(10a)
(12c) 1 March 2014 / Vol. 53, No. 7 / APPLIED OPTICS
1431
in which tr
2ni cosθi q ; ni cosθi n2t − n2i sin2 θi
tl
2ni nt cosθi q : n2t cosθi ni n2t − n2i sin2 θi
(12d)
The reflection and refraction events are controlled by the laws of reflection and retraction. The angles after reflection and refraction are given by the equations θr θi ;
nt sin θt ni sin θi ;
φr φt φi : (13)
We also need to compare the reflection coefficient ρ with a uniform pseudorandom number Rref to determine whether the reflection or refraction take place when the photon hits the Fresnel interface/boundary. For the Fresnel interface/boundary, the reflection coefficient is written as ρ C.
r2l r2r : 2
(14)
Photon Statistics and Post-processing
In order to get the real Stokes vector passing through some plane along some direction, we need to record the number of the photons passing through this plane along the θ; φ direction. In this work, we present the concept of VRTM to solve this problem. For a one-dimensional problem, the VRTM denoted by Di;θ;φ is defined as Di;θ;φ
Nt X Sp p1
N
;
(15)
in which N denotes the total number of emitting photons at the top surface of the multilayer system, as shown in Fig. 1; N t denotes the total number of the photons which pass through plane i from the θ; φ direction in the solid angle ΔΩ; Sp is the normalized Stokes vector of photon p. The physical meaning of VRTM is clear and it represents the transfer of intensity and polarized states of photons from the emission position to the receiving position. Finally, we obtain the real Stokes vector, written as Si;θ;φ
I 0 j cos θ0 j ·D ; j cos θjΔΩ i;θ;φ
(16)
in which I 0 is the intensity of incident beam and cos θ0 is the incident-direction cosine, as shown in Fig. 1. 3. Results and Discussion
In Section 2, we developed an arbitrary multilayer model for polarized radiative transfer. In Section 3, we first need to validate the correctness of the model. 1432
APPLIED OPTICS / Vol. 53, No. 7 / 1 March 2014
The MC model for polarized radiative transfer in the arbitrary multilayer system built in this work can be validated against the two-layer model for the atmosphere–ocean system. Then we will give the cases for the polarized radiative transfer in the three-layer, five-layer, and ten-layer system considering different optical boundary/interface conditions, refractive indices, optical thicknesses, and scattering albedos in different layers. Note that the Stokes vector is composed of direct light and diffuse light [20], and the results presented in this work are for the diffuse Stokes vector. According to Ref. [20], the direct Stokes vector is defined as Sdir S0 δcos θ − cos θ0 δφ − φ0 exp−τ∕ cos θ; (17a) where δ is the Dirac delta function; θ0 ; φ0 is the direction of the beam propagation without experiencing scattering. The diffuse Stokes vector is then defined as Sdif S − Sdir :
(17b)
It should be noted that the three layer, five-layer, and ten-layer cases are all fictitious problems, and do not have any physical significance. They are mainly used to test the computational implementation of the developed MC model for polarized radiative transfer in a multilayer (K > 2) system with complex properties. Case 1: A two-layer model of the atmosphere– ocean system In this case, we present the results compared against those obtained by the coupled atmosphere– water vector discrete ordinate radiative transfer (CAW-VDISORT) code in Ref. [29] to verify the MC polarized radiative transfer model with consideration of Fresnel interface developed in this work. The simulation is performed for polarized radiative transfer in a plane–parallel medium composed of Rayleigh scattering atmosphere layer and water body. A relative refractive index of water to the atmosphere is taken as 1.338. The optical depth is taken to be 0.15 for the atmosphere and 1.0 for the water. Both the atmosphere layer and water body are assumed to be nonabsorbing (ω 1.0), and the bottom of the water is assumed to be completely absorbing. A two-layer model is selected for the comparison. In the two-layer model, the top surface of the atmosphere and the bottom surface of the water are supposed to be Fresnel surfaces with n0 n1 and n3 n2 . The downward radiation field at the top of the atmosphere is taken to be a collimated and polarized beam of irradiance S0 1; 0; 0.5; 0.50 with an incident direction of θ0 ; φ0 60°; 0°. The number of the photons used in the MC simulation is 1.0 × 1010 and the angular discretization is taken as N θ × N φ 180 × 40 for this test. Figure 3 plots the curves for the four elements of the Stokes vector just below the water surface. Results for three different azimuth angles 0°, 90°,
1.4
×10
−1
4.0
×10
−2
ΜC ϕ = 0° ΜC ϕ = 90° ΜC ϕ = 180° Ref. [29] ϕ = 0° Ref. [29] ϕ = 90° Ref. [29] ϕ = 180°
3.0
1.2
2.0
1.0 ΜC ϕ = 0° ΜC ϕ = 90° ΜC ϕ = 180° Ref. [29] ϕ = 0° Ref. [29] ϕ = 90° Ref. [29] ϕ = 180°
0.6 0.4 0.2 -1.0 5.0
×10
-0.5 −2
0.0 cosθ
0.5
Q
I
1.0
0.8
0.0 -1.0 -2.0
1.0
-3.0 -1.0 5.0
4.0
4.0
3.0
3.0
2.0 0.0
0.5
1.0
0.5
1.0
−2
-1.0 -2.0
1.0 V
U
×10
0.0 cosθ
2.0 ΜC ϕ = 0° ΜC ϕ = 90° ΜC ϕ = 180° Ref. [29] ϕ = 0° Ref. [29] ϕ = 90° Ref. [29] ϕ = 180°
1.0
-1.0 -2.0 -3.0
-4.0
-4.0 -0.5
0.0 cosθ
0.5
ΜC ϕ = 0° ΜC ϕ = 90° ΜC ϕ = 180° Ref. [29] ϕ = 0° Ref. [29] ϕ = 90° Ref. [29] ϕ = 180°
0.0
-3.0 -5.0 -1.0
-0.5
1.0
-5.0 -1.0
-0.5
0.0 cosθ
Fig. 3. Stokes vector elements just below the atmosphere–water interface for a collimated and polarized incident beam.
and 180° are given for each element. As shown in Fig. 3, an excellent agreement between the results obtained by the MC model of this work and those by CAW-VDISORT [29] can be achieved. Case 2: A three-layer model with Fresnel surfaces/ interfaces In this case, we present a three-layer model with Fresnel surfaces/interfaces. The model is just for testing the computational implementation of our solution to polarized radiative transfer in a multilayer (K > 2) semitransparent medium with phase matrices more general than the pure Rayleigh scattering phase matrix considered in the first case. Note that the phase matrices used in layers 1 and 3 are presented in Ref. [20], including the cloud scattering and aerosol scattering. The Rayleigh scattering is considered in layer 2. The corresponding incident beam for these scattering matrices is at a wavelength of 0.412 μm. The refractive indices of the aerosol scattering, pure Rayleigh scattering, and cloud scattering layers are 1.385, 1, and 1.339, respectively. The refractive index background environment surrounding the three-layer system is equal to 1. Scattering albedos for layer 1, layer 2, and layer 3 are 0.3, 0.5, and 0.2, respectively; optical thickness for layer 1, layer 2, and layer 3 are 1, 0.3, and 0.5, respectively. The downward radiation field at the top of the multilayer medium (outside the top surface) is taken to be a collimated and polarized beam of irradiance S0 1; 0; 0.5; 0.50 with an incident direction of θ0 ; φ0 60°; 0°. The number of the photons used
in the model is 2.0 × 109 and the angular discretization is taken as N θ × N φ 90 × 20 for the simulation. Figures 4–6 plot the curves for the four Stokes vector elements at the interfaces of every layer. Results for two different azimuth angles 90° and 180° are given for each element. Resulting from the total reflection occurring at the interfaces of layers 1 and 3, Figs. 4 and 6 demonstrate the obvious peaks and valleys of abrupt change in the Stokes vector. At the two interfaces of layer 1, the peaks and valleys appear at the polar angle with a cosine of μ −0.67, and at the two interfaces of layer 3, they appear at μ 0.69. From Fig. 5, we can also find the peaks and valleys of the Stokes vector in the direction approximately parallel to the two interfaces of layer 2. Case 3: A three-layer model with Fresnel and Lambertian surface/interface In this case, we consider a three-layer model with Fresnel and Lambertian surface. The phase matrix for layer 1 and layer 3 is for 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 form of the phase matrix is as in Eq. (4a), and the Legendre coefficients for the four unique elements of the phase matrix have been obtained by Evens and Stephens [8]. Pure Rayleigh scattering is considered in layer 2, and the form of the phase matrix is given by Eq. (4b). In the three-layer system, the interface between layers 0 and 1 is a Lambertian one, and the other surfaces and interfaces are 1 March 2014 / Vol. 53, No. 7 / APPLIED OPTICS
1433
2.5
×10
−2
8.0
1.0
1-layer-top, ϕ = 90° 1-layer-bottom, ϕ = 90° 1-layer-top, ϕ = 180° 1-layer-bottom, ϕ = 180°
4.0 2.0 0.0
Q
I
1.5
−3
6.0
1-layer-top, ϕ = 90° 1-layer-bottom, ϕ = 90° 1-layer-top, ϕ = 180° 1-layer-bottom, ϕ = 180°
2.0
×10
-2.0
0.5
-4.0 0.0
-6.0
-0.5 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ 1.5
×10
-8.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
−2
0.4 1-layer-top, ϕ = 90° 1-layer-bottom, ϕ = 90° 1-layer-top, ϕ = 180° 1-layer-bottom, ϕ = 180°
1.0
−2
0.2 0.0 -0.2
1-layer-top, ϕ = 90° 1-layer-bottom, ϕ = 90° 1-layer-top, ϕ = 180° 1-layer-bottom, ϕ = 180°
V
U
0.5
×10
-0.4
0.0
-0.6 -0.5
-0.8
-1.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
-1.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
Fig. 4. Stokes vector elements at the top and bottom surfaces of layer 1 for the three-layer model.
3.5
×10
−2
1.0 2-layer-top, ϕ = 90° 2-layer-bottom, ϕ = 90° 2-layer-top, ϕ = 180° 2-layer-bottom, ϕ = 180°
3.0 2.5
×10
−2
2-layer-top, ϕ = 90° 2-layer-bottom, ϕ = 90° 2-layer-top, ϕ = 180° 2-layer-bottom, ϕ = 180°
0.8 0.6 0.4
I
Q
2.0 1.5
0.2
1.0
0.0
0.5
-0.2
0.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
-0.4 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
2.0 1.5 1.0
×10
−2
0.4 2-layer-top, ϕ = 90° 2-layer-bottom, ϕ = 90° 2-layer-top, ϕ = 180° 2-layer-bottom, ϕ = 180°
0.0 -0.2 -0.4 V
U -0.5
-0.6 -0.8 -1.0
-1.0 -1.5 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
−2
0.2
0.5 0.0
×10
-1.2
2-layer-top, ϕ = 90° 2-layer-bottom, ϕ = 90° 2-layer-top, ϕ = 180° 2-layer-bottom, ϕ = 180°
-1.4 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
Fig. 5. Stokes vector elements at the top and bottom surfaces of layer 2 for the three-layer model. 1434
APPLIED OPTICS / Vol. 53, No. 7 / 1 March 2014
3.5
×10
−2
1.4
2.5
0.8 0.6
1.5
0.4
1.0
0.2
0.5
0.0
0.0
-0.2
-0.5 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
1.0
U
0.5
×10
-0.4 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
−2
0.4
0.0 -0.5 -1.0 -1.5 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
×10
−2
0.2
3-layer-top, ϕ = 90° 3-layer-bottom, ϕ = 90° 3-layer-top, ϕ = 180° 3-layer-bottom, ϕ = 180°
0.0 -0.2 -0.4 V
1.5
3-layer-top, ϕ = 90° 3-layer-bottom, ϕ = 90° 3-layer-top, ϕ = 180° 3-layer-bottom, ϕ = 180°
1.0
Q
I
2.0
−2
1.2
3-layer-up, ϕ = 90° 3-layer-down, ϕ = 90° 3-layer-up, ϕ = 180° 3-layer-down, ϕ = 180°
3.0
×10
-0.6 -0.8 -1.0
3-layer-top, ϕ = 90° 3-layer-bottom, ϕ = 90° 3-layer-top, ϕ = 180° 3-layer-bottom, ϕ = 180°
-1.2 -1.4 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
Fig. 6. Stokes vector elements at the top and bottom surfaces of layer 3 for the three-layer model.
Fresnel ones. The refractive index of the background environment surrounding the multilayer system is equal to 1. The refractive indices of the pure Rayleigh scattering and Mie scattering layers are 1.338 and 1.44, respectively. Scattering albedos for layers 1, 2, and 3 are 1, 0.5, and 0.8, respectively; optical thickness for layers 1, 2, and 3 are 0.4, 0.3, and 0.3, respectively. The downward radiation field at the top of the three-layer system is taken to be a collimated and unpolarized beam of irradiance S0 1; 0; 0; 00 with an incident direction of θ0 ; φ0 60°; 0°. The number of the photons used for this case is 7.0 × 109 and the angular discretization is taken as N θ × N φ 90 × 20 for the simulation. Figures 7–9 show the curves for the four components of the Stokes vector at the top and bottom surfaces of every layer, and results for the azimuth angle φ 90° are presented. In Fig. 7, Q, U, and V components at the bottom surface of layer 1 for μ > 0 are also equal to zero due to the diffuse reflection on the Lambertian surface eliminating the polarization. The refractive index of layer 3 is larger than that of layer 2 and the background environment, so the total reflection will take place on the two Fresnel surfaces of layer 3 when the incident angle is greater than the critical angle. At the two surfaces of layer 3, the abrupt changes in the Stokes vector are observed at μ −0.72 and μ 0.37, as shown in Fig. 9. For layer 1, the total reflection takes place on the top surface, and the abrupt changes in the Stokes vector are found at μ −0.37 and μ −0.72; resulting from the
effect of total reflection occurring at the top surface of layer 1, similar changes in the Stokes vector can also observed at the bottom surface of layer 1, as shown in Fig. 7. Compared with the two adjacent layers, the refractive index of layer 2 is smaller, but compared with the background environments, the refractive index of layer 2 is bigger. So the Stokes vector at the surfaces of layer 2 is still influenced by the total reflection occurring at the surfaces of other layers. From Fig. 8, we can also see the abrupt changes in the Stokes vector at the two surfaces of layer 2 at μ −0.6644. Case 4: Five-layer and ten-layer models with Fresnel surfaces/interfaces In this case, we consider a five-layer model and a ten-layer model with Fresnel surfaces/interfaces. For the multilayer system, every layer has the same physical property parameters except for the refractive index. The total optical thickness of the multilayer system is 2.0, and the scattering albedo ω is equal to 0.5 for the multilayer system. The scattering type of the multilayer system is pure Rayleigh scattering. The refractive indices of the multilayer system are considered as a linear step distribution along the vertical direction. Two refractive index distributions are presented, and from the bottom layer to the top layer the distributions of refractive indices are two arithmetic progressions from 2.0 to 1.2 and from 1.2 to 2.0, respectively. The downward radiation field at the top of the multilayer system is taken to be a collimated and polarized beam of irradiance 1 March 2014 / Vol. 53, No. 7 / APPLIED OPTICS
1435
9.0 ×10
−2
3.0
×10
−3
2.0
8.0
1-layer-top,ϕ = 90° 1-layer-bottom,ϕ = 90°
1.0
7.0 Q
I
0.0
6.0
-1.0
5.0
-2.0
1-layer-top,ϕ = 90° 1-layer-bottom,ϕ = 90°
4.0
-3.0
3.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ 2.0
×10
-4.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
−2
4.0 3.0
1.0
2.0
−3
1-layer-top,ϕ = 90° 1-layer-bottom,ϕ = 90°
V
U
1.5
×10
0.5
1.0
0.0
0.0
1-layer-top,ϕ = 90° 1-layer-bottom,ϕ = 90°
-0.5 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
-1.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
Fig. 7. Stokes vector elements at the top and bottom surfaces of layer 1 for the three-layer model.
6.5
×10
−2
2.0
6.0
0.0
5.0
-1.0
4.5
Q
I
−3
1.0
5.5
-2.0
4.0
-3.0
3.5 2-layer-top,ϕ = 90° 2-layer-bottom,ϕ = 90°
3.0
×10
2-layer-top,ϕ = 90° 2-layer-bottom,ϕ = 90°
-4.0
2.5 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ 2.0
×10
-5.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
−2
6.0 5.0
1.5
×10
−3
2-layer-top,ϕ = 90° 2-layer-bottom,ϕ = 90°
4.0 3.0 V
U
1.0
2.0 0.5
0.0
1.0 2-layer-top,ϕ = 90° 2-layer-bottom,ϕ = 90°
-1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
0.0 -1.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
Fig. 8. Stokes vector elements at the top and bottom surfaces of layer 2 for the three-layer model. 1436
APPLIED OPTICS / Vol. 53, No. 7 / 1 March 2014
8.0
×10
−2
3.0
−3
2.0
7.0
3-layer-top,ϕ = 90° 3-layer-bottom,ϕ = 90°
1.0
6.0
0.0
5.0
-1.0
4.0
-2.0
Q
I
×10
-3.0
3.0
3-layer-top,ϕ = 90° 3-layer-bottom,ϕ = 90°
2.0
-4.0 -5.0
1.0
-6.0
0.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
-7.0 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
3.0
×10
−2
1.5
×10
−2
2.5 1.0
2.0
U
V
1.5 0.5
1.0 0.5 0.0
3-layer-top,ϕ = 90° 3-layer-bottom,ϕ = 90°
-0.5 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
0.0
3-layer-top,ϕ = 90° 3-layer-bottom,ϕ = 90°
-0.5 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 cosθ
Fig. 9. Stokes vector elements at the top and bottom surfaces of layer 3 for the three-layer model.
S0 1; 0.5; 0.5; 0.50 with an incident direction of θ0 ; φ0 30°; 0°. The number of the photons used for these two cases is 5.0 × 109 and the angular discretization is taken as N θ × N φ 90 × 20 for the simulations. The results only for φ 90° are plotted. Figures 10 and 11 illustrate the curves for the four elements I, Q, U, and V of the Stokes vector at the bottom surface and the top surface of the five-layer model, respectively. It is well known that, as light beams enter into the medium with smaller refractive index from the medium with bigger refractive index, total reflection will happen at the interface only if the incident angle is not smaller than the critical angle of total reflection. It should also be noted that the scattering can change the direction of the light propagation. Based on these two points, we can give a qualitative analysis of the effects of total reflection occurring at the interfaces on the Stokes vector at the top and bottom surfaces of the multilayer system. In the multilayer system with a distribution of refractive indices increasing along the incident direction of external beams, Stokes vector at the top surface of the fifth layer is mainly affected by the total reflection at the top surface corresponding to the critical-angle cosine of cos θ ≈ −0.553, which gives rise to a big abrupt change in curves of Stokes vector just at the direction with cos θ ≈ −0.553, as shown in Fig. 11. In fact, the Stokes vector at the top surface can also be affected by the total reflection from the bottom surface, although the effect is so weak that it can be neglected. At the bottom surface,
the Stokes vector is affected by the total reflection from all six interfaces in the multilayer system. With respect to the bottom surface of the first layer with n1 2, corresponding to the critical angles for total reflection at the top surfaces of layers 1 to 4, the bottom surface of layer 1, and the top surface of layer 5, the cosines of the incident angles can be obtained, equal to 0.436, 0.6, 0.714, 0.8, and 0.866, respectively, where abrupt changes happen in the curves of the Stokes vector at the bottom surface of the multilayer system, and consequently, in the ranges of [0.436, 0.866] and −0.866; −0.436, fluctuations like waves can be seen, as shown in Fig. 10. In the multilayer system with a distribution of refractive indices decreasing along the incident direction of external beams, a big abrupt change in the Stokes vector at the bottom surface can be found at the direction with cos θ ≈ 0.553, as shown in Fig. 10, and abrupt changes in the Stokes vector at the top surface can be observed at the directions with 0.436, 0.6, 0.714, 0.8, and 0.866, respectively, as shown in Fig. 11. Figures 12 and 13 plot the curves of the Stokes vector at the top and bottom surfaces of the ten-layer system for the two distributions of refractive indices. From Figs. 10 to 13, we can see that the curves of the Stokes vector for the ten-layer model are similar to those for the five-layer model, and compared with the curves of the Stokes vector for the five-layer model, the fluctuations like waves in the curves of the Stokes vector for the ten-layer model are of 1 March 2014 / Vol. 53, No. 7 / APPLIED OPTICS
1437
−2
−2
4.5 ×10
0.0 ×10
4.0 3.5
nbottom=2.0, ntop=1.2, ϕ=90°
3.0
nbottom=1.2, ntop=2.0, ϕ=90°
-0.5 -1.0
2.5
Q
I
z/L=0
2.0
nbottom=2.0, ntop=1.2, ϕ=90°
-1.5
nbottom=1.2, ntop=2.0, ϕ=90°
1.5 1.0
-2.0
z/L=0
0.5 0.0 -1.0
-0.5 −2
×10
0.0 cosθ
0.5
-2.5 -1.0
1.0
-0.5 −2
1.4 ×10
0.5
0.5
1.0
1.2 1.0
0.0
0.8 -0.5
nbottom=2.0, ntop=1.2, ϕ=90°
V
U
0.0 cosθ
nbottom=1.2, ntop=2.0, ϕ=90°
nbottom=2.0, ntop=1.2, ϕ=90° nbottom=1.2, ntop=2.0, ϕ=90°
0.6 0.4
-1.0
z/L=0
0.2
z/L=0
0.0
-1.5 -1.0
-0.5
0.0 cosθ
0.5
-0.2 -1.0
1.0
-0.5
0.0 cosθ
0.5
1.0
Fig. 10. Stokes vector elements at the bottom interface for the five-layer model.
−2
−2
6.0 ×10
0.0 ×10
5.0
nbottom=2.0, ntop=1.2, ϕ=90°
-1.0
nbottom=1.2, ntop=2.0, ϕ=90°
4.0 z/L=1.0
3.0
nbottom=2.0, ntop=1.2, ϕ=90°
Q
I
-2.0 -3.0
nbottom=1.2, ntop=2.0, ϕ=90°
2.0
z/L=1.0
-4.0
1.0 0.0 -1.0
-0.5
0.0 cosθ
0.5
-5.0 -1.0
1.0
−2
-0.5
0.0 cosθ
0.5
1.0
−2
4.0 ×10
3.0 ×10
3.0
2.0
nbottom=2.0, ntop=1.2, ϕ=90° nbottom=1.2, ntop=2.0, ϕ=90°
1.0
z/L=1.0
V
U
2.0 nbottom=2.0, ntop=1.2, ϕ=90°
0.0
nbottom=1.2, ntop=2.0, ϕ=90°
1.0
-1.0
z/L=1.0
0.0 -1.0
-0.5
0.0 cosθ
0.5
1.0
-2.0 -1.0
-0.5
0.0 cosθ
Fig. 11. Stokes vector elements at the top interface for the five-layer model. 1438
APPLIED OPTICS / Vol. 53, No. 7 / 1 March 2014
0.5
1.0
−2
−2
5.0 ×10
0.0 ×10
4.0
-0.5 nbottom=2.0, ntop=1.2, ϕ=90° nbottom=1.2, ntop=2.0, ϕ=90°
-1.0 Q
I
3.0
z/L=0
2.0
nbottom=2.0, ntop=1.2, ϕ=90°
-1.5
nbottom=1.2, ntop=2.0, ϕ=90°
1.0
z/L=0
-2.0
0.0 -1.0
-0.5 −2
0.0 cosθ
0.5
-2.5 -1.0
1.0
-0.5 −2
1.0 ×10
0.0 cosθ
0.5
1.0
1.6 ×10 1.4
0.5
1.2 1.0
0.0 V
U
nbottom=2.0, ntop=1.2, ϕ=90° nbottom=1.2, ntop=2.0, ϕ=90°
-0.5
0.8
nbottom=2.0, ntop=1.2, ϕ=90°
0.6
nbottom=1.2, ntop=2.0, ϕ=90° z/L=0
0.4
z/L=0
0.2
-1.0
0.0 -1.5 -1.0
-0.5
0.0 cosθ
0.5
-0.2 -1.0
1.0
-0.5
0.0 cosθ
0.5
1.0
Fig. 12. Stokes vector elements at the bottom interface for the ten-layer model.
−2
−2
6.0 ×10
0.0 ×10
5.0
nbottom=2.0, ntop=1.2, ϕ=90° nbottom=1.2, ntop=2.0, ϕ=90°
-1.0
z/L=1.0
4.0 Q
I
-2.0 3.0 nbottom=2.0, ntop=1.2, ϕ=90°
2.0
-3.0
nbottom=1.2, ntop=2.0, ϕ=90°
-4.0
z/L=1.0
1.0 0.0 -1.0
-0.5
0.0
0.5
-5.0 -1.0
1.0
cosθ
−2
×10
-0.5
0.0
0.5
1.0
cosθ
−2
3.0 ×10
4.0 nbottom=2.0, ntop=1.2, ϕ=90°
2.0
nbottom=1.2, ntop=2.0, ϕ=90°
3.0
z/L=1.0
1.0 V
U
2.0 nbottom=2.0, ntop=1.2, ϕ=90°
1.0 0.0 -1.0 -1.0
0.0
nbottom=1.2, ntop=2.0, ϕ=90°
-1.0
z/L=1.0
-0.5
0.0 cosθ
0.5
1.0
-2.0 -1.0
-0.5
0.0 cosθ
0.5
1.0
Fig. 13. Stokes vector elements at the top interface for the ten-layer model. 1 March 2014 / Vol. 53, No. 7 / APPLIED OPTICS
1439
higher frequency and smaller amplitudes. Such calculations took 6 h 2 min and 7 h 34 min of CPU time for the five-layer and ten-layer cases, respectively. The program was running at 3.0 GHz on a single core of an Intel Core i5 processor. The computation time increases significantly with the number of layers increasing, since reflection and transmission events increase many times. The ten-layer model is closer to a continuum, while it needs more computation time to obtain converged results. It is expected that the approximate solution to vector radiative transfer will be implemented in a gradient-index medium by dividing the medium into a sufficient number of sublayers. Based on the computation time spent by the ten-layer case, more or unacceptable time may be needed for the gradient-index problem. Optimizing the code to reduce the computation time will be our next work. 4. Conclusions
In this work, we developed a MC model for solving arbitrary multilayer vector radiative transfer problems with Fresnel and Lambertian surfaces/ interfaces. We considered two kinds of polarization mechanisms, scattering by particles and reflection/ refraction by Fresnel surfaces/interfaces. We presented the detailed steps of determining scattering angles by using the modified rejection method. We put forward the concept of VRTM for describing the polarized radiative transfer. Using the VRTM solved by the MC method, we obtained four components of the Stokes vector. By comparing against the results of the two-layer model for the atmosphere–ocean system in the literature, we verified our multilayer MC method model with high accuracy and flexible adaptability. By comparison, we can see that all of the four components of the Stokes vector obtained by MC simulation are in good agreement with the results by DOM reported in the published literature. Next, we presented the results for two three-layer models with different optical boundary conditions and scattering phase matrices. Fresnel and Lambertian surfaces/interfaces were considered, and Rayleigh scattering and different types of Mie scattering from the literature were chosen for our models. Finally, two numerical cases of five-layer and ten-layer models were also presented at the top and bottom surfaces. By analyzing the results, we find the optical boundary/interface conditions have a significant effect on the polarized radiative transfer. In the multilayer system, due to inhomogeneous distribution of refractive index, the complex total reflection occurs at the surfaces/interfaces when the incident angle is equal to or greater than the critical angle, which generates abrupt changes in Stokes vector. Resulting from the total reflection at the surfaces/interfaces, fluctuations like waves can be observed in the curves of the Stokes vector at the surfaces/interfaces. With the increasing of the number of layers and the critical angles for total reflection, the frequency of fluctuations increases 1440
APPLIED OPTICS / Vol. 53, No. 7 / 1 March 2014
and the amplitudes of fluctuations decrease. The novelty of our work is that the arbitrary multilayer model developed can be applied to simulating vector radiative transfer in the inhomogeneous media which may be continuum or noncontinuum. For instance, it may be used to implement the approximate solution to vector radiative transfer in a gradientindex medium. In this work, we have made some attempts with the five-layer and ten-layer models to try to solve this problem. However, it is clear that we do not obtain smooth and stable solutions to polarized radiative transfer in a gradient-index medium with such a small number of layers. This work is supported by the National Natural Science Foundation of China (Grant Nos. 51176040 and 51276050). References 1. L. Tsang and J. A. Kong, “Radiative transfer theory for active remote sensing of half-space random media,” Radio Sci. 13, 763–773 (1978). 2. L. Tsang, J. A. Kong, and R. T. Shin, Theory of Microwave Remote Sensing (Wiley, 1985). 3. L. V. Wang and H. Wu, Biomedical Optics: Principles and Imaging (Wiley, 2007). 4. N. Ghosh, M. F. G. Wood, and I. A. Vitkin, Handbook of Photonics for Biomedical Science, V. V. Tuchin, ed. (Taylor and Francis, 2010), Chap. 9, pp. 253–282. 5. M. P. Mengüç and S. Manickavasagam, “Characterization of size and structure of agglomerates and inhomogeneous particles via polarized light,” Int. J. Eng. Sci. 36, 1569–1593 (1998). 6. G. W. Kattawar and G. N. Plass, “Radiance and polarization of multiple scattered light from haze and clouds,” Appl. Opt. 7, 1519–1527 (1968). 7. R. D. M. Garcia and C. E. Siewert, “The FN method for radiative transfer models that include polarization effects,” J. Quant. Spectrosc. Radiat. Transfer 41, 117–145 (1989). 8. K. F. Evans and G. L. Stephens, “A new polarized atmospheric radiative transfer model,” J. Quant. Spectrosc. Radiat. Transfer 46, 413–423 (1991). 9. C. E. Siewert, “A discrete-ordinates solution for radiativetransfer models that include polarization effects,” J. Quant. Spectrosc. Radiat. Transfer 64, 227–254 (2000). 10. J. Lenoble, M. Herman, J. L. Deuzé, B. Lafrance, R. Santer, and D. Tanré, “A successive order of scattering code for solving the vector equation of transfer in the earth’s atmosphere with aerosols,” J. Quant. Spectrosc. Radiat. Transfer 107, 479–507 (2007). 11. F. Xu, R. A. West, and A. B. Davis, “A hybrid method for modeling polarized radiative transfer in a spherical-shell planetary atmosphere,” J. Quant. Spectrosc. Radiat. Transfer 117, 59–70 (2013). 12. H. Ishimoto and K. Masuda, “A Monte Carlo approach for the calculation of polarized light: application to an incident narrow beam,” J. Quant. Spectrosc. Radiat. Transfer 72, 467–483 (2002). 13. X. D. Wang and L. V. Wang, “Propagation of polarized light in birefringent turbid media: a Monte Carlo study,” J. Biomed. Opt. 7, 279–290 (2002). 14. M. Xu, “Electric field Monte Carlo simulation of polarized light propagation in turbid media,” Opt. Express 12, 6530–6539 (2004). 15. R. Vaillon, B. T. Wong, and M. P. Mengüç, “Polarized radiative transfer in a particle-laden semi-transparent medium via a vector Monte Carlo method,” J. Quant. Spectrosc. Radiat. Transfer 84, 383–394 (2004). 16. C. Davis, C. Emde, and R. Harwood, “A 3-D polarized reversed Monte Carlo radiative transfer model for millimeter and submillimeter passive remote sensing in cloudy atmospheres,”
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
IEEE Trans. Geosci. Remote Sens. 43, 1096–1101 (2005). J. C. Ramella-Roman, S. A. Prahl, and S. L. Jacques, “Three Monte Carlo programs of polarized light transport into scattering media: part I,” Opt. Express 13, 4420–4438 (2005). J. N. Swamy, C. Crofcheck, and M. P. Mengüç, “A Monte Carlo ray tracing study of polarized light propagation in liquid foams,” J. Quant. Spectrosc. Radiat. Transfer 104, 277–287 (2007). C. Cornet, L. C. Labonnote, and F. Szczap, “Three-dimensional polarized Monte Carlo atmospheric radiative transfer model (3DMCPOL): 3D effects on polarized visible reflectances of a cirrus cloud,” J. Quant. Spectrosc. Radiat. Transfer 111, 174–186 (2010). A. A. Kokhanovsky, V. P. Budak, C. Cornet, M. Z. Duan, C. Emde, I. L. Katsev, D. A. Klyukov, S. V. Korkin, L. C. Labonnote, B. Mayer, Q. L. Min, T. Nakajima, Y. Ota, A. S. Prikhach, V. V. Rozanov, T. Yokota, and E. P. Zege, “Benchmark results in vector atmospheric radiative transfer,” J. Quant. Spectrosc. Radiat. Transfer 111, 1931–1946 (2010). G. W. Kattawar, G. N. Plass, and J. A. Guinn, “Monte Carlo calculations of the polarization of radiation in the Earth’s atmosphere-ocean system,” J. Phys. Oceanogr. 3, 353–372 (1973). G. W. Kattawar and C. N. Adams, “Stokes vector calculations of the submarine light field in an atmosphere–ocean with scattering according to a Rayleigh phase matrix: effect of interface refractive index on radiance and polarization,” Limnol. Oceanogr. 34, 1453–1472 (1989). P. W. Zhai, Y. X. Hu, C. R. Trepte, and P. L. Lucker, “A vector radiative transfer model for coupled atmosphere and ocean systems based on successive order of scattering method,” Opt. Express 17, 2057–2079 (2009). P. W. Zhai, Y. X. Hu, J. Cdhary, C. R. Trepte, P. L. Lucker, and D. B. Josset, “A vector radiative transfer model for coupled atmosphere and ocean systems with a rough interface,” J. Quant. Spectrosc. Radiat. Transfer 111, 1025–1040 (2010). P. W. Zhai, Y. X. Hu, D. B. Josset, C. R. Trepte, P. L. Lucker, and B. Lin, “Advanced angular interpolation in the vector radiative transfer for coupled atmosphere and ocean systems,” J. Quant. Spectrosc. Radiat. Transfer 115, 19–27 (2013). X. Q. He, D. L. Pan, Y. Bai, Q. K. Zhu, and F. Gong, “Vector radiative transfer numerical model of coupled oceanatmosphere system using matrix-operator method,” Sci. China Ser. D 50, 442–452 (2007). X. Q. He, Y. Bai, Q. Zhu, and F. Gong, “A vector radiative transfer model of coupled ocean-atmosphere system using
28.
29.
30.
31. 32. 33.
34. 35. 36. 37. 38.
39. 40. 41. 42.
matrix-operator method for rough sea-surface,” J. Quant. Spectrosc. Radiat. Transfer 111, 1426–1448 (2010). Y. Ota, A. Higurashi, T. Nakajima, and T. Yokota, “Matrix formulations of radiative transfer including the polarization effect in a coupled atmosphere-ocean system,” J. Quant. Spectrosc. Radiat. Transfer 111, 878–894 (2010). E. R. Sommersten, J. K. Lotsberg, K. Stamnes, and J. J. Stamnes, “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. Transfer 111, 616–633 (2010). R. D. M. Garcia, “Fresnel boundary and interface conditions for polarized radiative transfer in a multilayer medium,” J. Quant. Spectrosc. Radiat. Transfer 113, 306–317 (2012). P. W. Zhai, G. W. Kattawar, and Y. X. Hu, “Comment on transmission matrix for dielectric interface,” J. Quant. Spectrosc. Radiat. Transfer 113, 1981–1984 (2012). R. D. M. Garcia, “Response to ‘Comment on transmission matrix for a dielectric interface’,” J. Quant. Spectrosc. Radiat. Transfer 113, 2251–2254 (2012). R. D. M. Garcia, “Radiative transfer with polarization in a multi-layer medium subject to Fresnel boundary and interface conditions,” J. Quant. Spectrosc. Radiat. Transfer 115, 28–45 (2013). S. Chandrasekhar, Radiative Transfer (Oxford University, 1950). H. C. Van de Hulst, Light Scattering by Small Particles (Dover, 1981). M. I. Mishchenko, Scattering, Absorption, and Emission of Light by Small Particles (NASA, 2002). A. T. Young, “Revised depolarization corrections for atmospheric extinction,” Appl. Opt. 19, 1333 (1980). H. H. Tynes, G. W. Kattawar, E. P. Zege, I. L. Katsev, A. S. Prikhach, and L. I. Chaikovskaya, “Monte Carlo and multicomponent approximation methods for vector radiative transfer by use of effective Mueller matrix calculations,” Appl. Opt. 40, 400–412 (2001). R. M. Green, Spherical Astronomy (Cambridge University, 1985). B. A. Whitney, “Monte Carlo radiative transfer,” arXiv: 1104.4990 (2011). R. Siegel and C. M. Spuckler, “Refractive index effects on radiation in an absorbing, emitting, and scattering laminated layer,” J. Heat Transfer 115, 194–200 (1993). K. Masuda and T. Takashima, “Computational accuracy of radiation emerging from the ocean surface in the model atmosphere–ocean system,” Pap. Met. Geophys. 37, 1–13 (1986).
1 March 2014 / Vol. 53, No. 7 / APPLIED OPTICS
1441
