Numerical test of an inverse polarized radiative transfer algorithm

Journal of Quantitative Spectroscopy & Radiative Transfer 78 (2003) 235 – 241

www.elsevier.com/locate/jqsrt

Numerical test of an inverse polarized radiative transfer algorithm T. Viika;∗ , N.J. McCormickb a

b

Tartu Observatory, Toravere, Tartumaa 61602, Estonia University of Washington, Mechanical Engineering, Seattle, WA 98195-2600, USA Received 6 June 2002; accepted 3 September 2002

Abstract A procedure is tested with which to determine the single-scattering albedo from polarization measurements of the angle-dependent intensity at two locations within, or on the boundaries of, a homogeneous 6nite or in6nite atmosphere that scatters radiation according to the Rayleigh law with true absorption. ? 2003 Elsevier Science Ltd. All rights reserved. Keywords: Radiative transfer; Polarization; Rayleigh scattering

1. Introduction Polarized radiative transfer for Rayleigh scattering is especially simple because the analysis can be broken up into a matrix equation for the analysis of the azimuthally averaged light components I‘ (; ) and Ir (; ) in the two-vector I(; ) plus scalar equations for U (; ) and V (; ). The classic Stokes parameters are given by I (; ) = I‘ (; ) + Ir (; ) and Q(; ) = I‘ (; ) − Ir (; ). The equation of transfer for I(; ) that we consider here is [1,2]
1 @   I(; ) + I(; ) = Q() Q( )I(;  ) d ; 0 6  6 0 ; (1) @ 2 −1 where  is measured from the positive axis of the optical depth . The matrix Q() is de6ned as [2–4]  2 2  c + 3 (1 − c) (2c)1=2 (1 − 2 ) 3 Q() = : (2) 1 2(c + 2)1=2 (c + 2) 0 3

∗

Corresponding author. Tel.: +372-7-410-154; fax: +372-7-410-205. E-mail addresses: [email protected] (T. Viik), [email protected] (N.J. McCormick).

0022-4073/03/$ - see front matter ? 2003 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 2 - 4 0 7 3 ( 0 2 ) 0 0 2 1 5 - 7

236 T. Viik, N.J. McCormick / Journal of Quantitative Spectroscopy & Radiative Transfer 78 (2003) 235 – 241

The general scattering law (=2)Q()QT (), for QT () the transpose of Q(), is for the singlescattering albedo  ∈ [0; 1] and parameter c ∈ [0; 1]. Two special cases of this law are the purely Rayleigh-scattering model [1,5] for c = 1 and the conservative model for a mixture of Rayleigh- and isotropic-scattering laws [1] for  = 1 and c ∈ [0; 1]. The objective of this paper is to numerically test two explicit inverse radiative transfer equations for determining  and c from measurements of I(; ); −1 6  6 1, at depths 1 and 2 for 0 6 1 ¡ 2 6 0 6 ∞. Two problems will be considered—the planetary problem with boundary conditions I(0; ) = U( − 0 ); I(0 ; ) = 0;

0 6  6 1;

−1 6  6 0;

(3)

for unit matrix U, and the Milne problem with a source deep in the interior, where I(0; ) = 0;

0661

(4)

and the intensity deep in the interior increases no more rapidly than exp(). The inverse equations for this problem were derived with the same procedures applied to derive the inverse equations for general anisotropic scattering [6,7] and are [8]:  = [P0T (1 )P0 (1 ) − P0T (2 )P0 (2 )]−1 [S0 (1 ) − S0 (2 )]

(5)

[P1T (1 )(U − R)−1 P1 (1 ) − P1T (2 )(U − R)−1 P1 (2 )] = S1 (1 ) − S1 (2 ):

(6)

and Here 1 ¿ 0 and 2 6 0 are the locations of measurements and the vector Pn () is
1 n QT ()I(; ) d; n = 0; 1 Pn () =

(7)

and the scalar Sn () is
1 2n IT (; −)I(; ) d; Sn () = 4

(8)

−1

0

and

 (U − R)−1 = −1

n = 0; 1

c + 2 − 12c=5

(2c)1=2 (1 − 7c=10)

(2c)1=2 (1 − 7c=10)

c + 2 − (2 + 7c2 =10)

 (9)

for  = (1 − )(c + 2)(1 − 7c=10). Eq. (5) is a linear equation, independent of c, from which  can be immediately computed, whereas Eq. (6) is a highly nonlinear equation that must be evaluated for c once  has been computed. 2. Numerical tests The numerical tests were performed by solving the forward radiative transfer problem for I(; ) at Gauss nodes with the procedure outlined earlier [9–12]. Then those results were used for the mean values about which random errors were introduced by Monte Carlo sampling of 10,000

T. Viik, N.J. McCormick / Journal of Quantitative Spectroscopy & Radiative Transfer 78 (2003) 235 – 241 237 1.0


obs

0.8

0 = 0.10 c = 0.10  = 0.10 1 = 0.0 2 = 1.0

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4 0.4

0.5

0 = 1.00 c = 0.10  = 0.10 1 = 0.0 2 = 1.0

0.9


obs

0.9

1.0

0.6

0.7

0.8

0.9

1.0


calc

Fig. 1. obs obtained from Eq. (5) versus calc of forward problem for c = 0:1 with error factor  = 0:1 and planetary illumination 0 = 0:1.

0.4 0.4

0.5

0.6

0.7
calc

0.8

0.9

1.0

Fig. 2. Same as Fig. 1 except 0 = 1:0.

“measurements”. The numerically simulated errors were obtained by using the RAN1 program [13] to generate pseudorandom numbers uniformly distributed on (0,1). Then the Box–Muller method [13] was used to obtain a normal distribution of the errors, (x) =  exp(−x2 =22 );

(10)

where  is the dispersion of random individual errors and  is the amplitude error factor (i.e.,  = 0:01 corresponds to a magnitude of 1%). The errors were added to the intensities computed for the forward problem. Eq. (5) was 6rst solved for the value of , which could be done even if the error amplitude was a∼0:25– 0.3. Then the value of  was used in Eq. (6) to 6rst determine the number of possible values of c using BRAK [13] and then the values were obtained using BRENT [13]. In instances where there were two or three roots we had to use each c-value in our model to determine which minimized the 1 values of −1 [IT (2 ; )I(2 ; ) − IT (1 ; )I(1 ; )] d, where I(; ) is the diMerence between measured values and those computed with the selected value of c. For the numerical tests, the values of 1 = 0 and 2 = 0 = 1 were used for the planetary problem, but similar results were obtained for 1 ¿ 0 and values 2 ¡ 1 unless the two locations were very close to each other. Values of 1 = 0 and 2 = 0 = ∞ were selected for the Milne problem. The following 6gures show values observed from running the inversion algorithm, obs , versus the corresponding values used in the forward problem calculations, calc , or similarly cobs versus ccalc . Each line of mean values obtained from the inversion algorithm, with I(; ) corrupted by random noise, is bounded on either side by the values of the mean plus or minus two standard deviations. Figs. 1 and 2 are for the planetary atmosphere with c = 0:1, a 10% error factor (i.e.,  = 0:1), and 0 = 0:1 and 1.0, respectively. The results are very insensitive to the direction of the incident

238 T. Viik, N.J. McCormick / Journal of Quantitative Spectroscopy & Radiative Transfer 78 (2003) 235 – 241 1.0

1.0


obs

0.8

0.9

0.7

0.6

0.5

0.5

0.5

0.6

0.7
calc

0.8

0.9

0.4 0.4

1.0

Fig. 3. Same as Fig. 1 except c = 1:0 and 0 = 0:5.

0.8 0.7 0.6

0.6

0.7
calc

0.8

0.9

1.0

1.0 0 = 0.10
= 0.40  = 0.01 1 = 0.0 2 = 1.0

0.5

0.9 0.8 0.7

0 = 0.50
= 0.40  = 0.01 1 = 0.0 2 = 1.0

0.6 cobs

0.9

0.5

Fig. 4. Same as Fig. 1 except for Milne problem.

1.0

cobs

0.7

0.6

0.4 0.4

c = 1.00  = 0.10 1 = 0.0 2 = 1.0

0.8
obc

0.9

0 = 0.50 c = 1.00  = 0.10 1 = 0.0 2 = 1.0

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ccalc

0.0 0.0 0.1 0.2 0.3

Fig. 5. cobs obtained from Eq. (6) versus ccalc of forward problem for =0:4 with error factor =0:01 and planetary illumination 0 = 0:1.

0.4 0.5 0.6 0.7 0.8 0.9 1.0 ccalc

Fig. 6. Same as Fig. 5 0 = 0:5.

illumination. For a planetary atmosphere with c = 1 and 0 = 0:5, shown in Fig. 3, and for the corresponding Milne problem in Fig. 4, it is seen that the dispersion of the results is less for the Milne problem. The algorithm for determining the single scattering albedo was found to be very stable, although it is surprising that the dispersion bounds of the output are so broad. Figs. 5–7 show the results for the planetary atmosphere with  = 0:4;  = 0:01, and 0 = 0:1; 0:5, and 1.0, respectively. Even for only 1% errors, the results in Fig. 7 are quite poor for c larger than ∼0:7. To see if the poor results of Fig. 7 were due to the value of  selected, we also tested

T. Viik, N.J. McCormick / Journal of Quantitative Spectroscopy & Radiative Transfer 78 (2003) 235 – 241 239 1.0 0.9 0.8 0.7

1.0 0 = 1.00
= 0.40  = 0.01 1 = 0.0 2 = 1.0

0.9 0.8 0.7 0.6 cobs

cobs

0.6 0.5

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0.0 0.0 0.1 0.2 0.3

0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ccalc

0.4 0.5 0.6 0.7 0.8 0.9 1.0 ccalc

Fig. 7. Same as Fig. 5 except 0 = 1:0.

Fig. 8. Same as Fig. 5 except  = 0:7 and 0 = 0:5.

1.0

1.0

0.8 0.7

0.9

0 = 1.00
= 0.70  = 0.01 1 = 0.0 2 = 1.0

0.8 0.7 0.6

0.5

cobs

0.6 cobs

0.5

0.4

0.9

0 = 0.50
= 0.70  = 0.01 1 = 0.0 2 = 1.0

0 = 0.10
= 0.99  = 0.01 1 = 0.0 2 = 1.0

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ccalc

0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ccalc

Fig. 9. Same as Fig. 5 except  = 0:7 and 0 = 1:0.

Fig. 10. Same as Fig. 5 except  = 0:99.

the algorithm for  = 0:7 with 0 = 0:5 and 1.0, as shown in Figs. 8 and 9, and for  = 0:99 in Figs. 10–12. It is surprising to see that for  = 0:99 the results for 0 = 0:5 are worse than those for 0 = 0:1 and 1.0. We conclude that the value of c for the planetary atmosphere problem is diOcult to determine for values above ∼0:7. Fig. 13 shows, however, that the value of c for the Milne problem can be reasonably extracted up to c∼0:9.

240 T. Viik, N.J. McCormick / Journal of Quantitative Spectroscopy & Radiative Transfer 78 (2003) 235 – 241 1.0 0.9 0.8 0.7

1.0 0 = 0.50
= 0.99  = 0.01 1 = 0.0 2 = 1.0

0.9 0.8 0.7 0.6 cobs

cobs

0.6 0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0.0 0.0 0.1

0 = 1.00
= 0.99  = 0.01 1 = 0.0 2 = 1.0

0.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ccalc

0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ccalc

Fig. 11. Same as Fig. 5 except  = 0:99 and 0 = 0:5.

Fig. 12. Same as Fig. 5 except  = 0:99 and 0 = 1:0.

1.0

0.9


obs

0.8

c = 1.00  = 0.10 1 = 0.0 2 = 1.0

0.7

0.6

0.5

0.4 0.4

0.5

0.6

0.7

0.8

0.9

1.0


calc

Fig. 13. cobs obtained from Eq. (6) versus ccalc of forward problem for  =0:9 with error factor =0:01 for Milne problem.

The numerical instabilities encountered in solving Eq. (6) for determining c are analogous to those encountered by Siewert and Maiorino [14] who considered a 6nite atmosphere with c = 1 with a Lambertian rePecting surface. They too found that  could be more accurately determined than their second unknown, the Lambertian rePection coeOcient. It is tempting to ask whether it might be possible to use the Stokes parameters U (; ) or V (; ) to develop an inverse equation to replace Eq. (6) that does not work well. But for a

T. Viik, N.J. McCormick / Journal of Quantitative Spectroscopy & Radiative Transfer 78 (2003) 235 – 241 241

planetary atmosphere illuminated by natural light then V (; ) ≡ 0, so we are left with considering whether information can be obtained from U (; ). Unfortunately, however, U (; ) also vanishes unless U also depends on the azimuthal angle  through an azimuthally nonsymmetric boundary condition. Acknowledgements This work has been supported by the Estonian Ministry of Education within Project No. TO 0060059S98 and by the Estonian Science Foundation under Grant No. 4701. References [1] Chandrasekhar S. Radiative transfer. Oxford: Oxford University Press, 1950; New York: Dover, 1960. [2] Siewert CE, Maiorino JR. The complete solution for the scattering of polarized light in a Rayleigh and isotropically scattering atmosphere. Astrophys Space Sci 1980;72:189–201. [3] Bond GR, Siewert CE. On the nonconservative equation of transfer for a combination of Rayleigh and isotropic scattering. Astrophys J 1971;164:97–110. [4] Siewert CE, Burniston EE. On existence and uniqueness theorems concerning the H -matrix of radiative transfer. Astrophys J 1972;174:629–41. [5] Schnatz TW, Siewert CE. On the transfer of polarized light in Rayleigh-scattering half spaces with true absorption. Mon Not R Astron Soc 1971;152:491–508. [6] Siewert CE. Solutions to an inverse problem in radiative transfer with polarization—I. JQSRT 1983;30:523– 6. [7] McCormick NJ, Sanchez R. Solutions to an inverse problem in radiative transfer with polarization—II. JQSRT 1983;30:527–35. [8] McCormick NJ. Determination of the single-scattering albedo of a dense Rayleigh-scattering atmosphere with true absorption. Astrophys Space Sci 1980;71:235–8. [9] Viik T. Rayleigh scattering in planetary atmospheres. Earth Moon Planets 1989;46:261–84. [10] Viik T. The Pomraning phase function: how good is it? JQSRT 2000;68:319–26. [11] Viik T. Rayleigh–Cabannes scattering in planetary atmospheres II. The Milne problem in conservative atmospheres. Earth Moon Planets 1990;49:163–75. [12] Viik T. Radiation 6eld in Rayleigh–Cabannes scattering atmospheres. The nonconservative Milne problem. JQSRT 2000;66:581–90. [13] Press HW, Flannery BP, Teukolsky SA, Vetterling TW. Numerical recipes. Cambridge: Cambridge University Press, 1986. [14] Siewert CE, Maiorino JR. The inverse problem for a 6nite Rayleigh-scattering atmosphere. J Appl Math Phys (ZAMP) 1980;31:767–70.