Mathematical Model of a Surface Radiance Factor A.Y. Basov1, V.P. Budak1 callia12@rambler.ru, budakvp@gmail.com 1 National Research University "Moscow Power Engineering Institute", Moscow, Russia The article is devoted to the creation of a surface radiance factor mathematical model. The basis of the model is the solution of the boundary value problem of the radiative transfer equation (RTE). The surface is considered as a structure consisting of several turbid layers, each of which is characterized by its optical parameters. The top of the structure is randomly rough, uncorrelated, Fresnel. The lower boundary reflects perfectly diffusely. The complexity of solving the RTE boundary value problem for real layers is due to the fact that the suspended particles in each layer are always much longer than the wavelength. This leads to a strong anisotropy of the radi- ance angular distribution according to Mie theory. The solution comes down to a system of equations by the discrete ordinates method that consists of several hundred of differential equations. Subtraction of the anisotropic part from the solution based on an approxi- mate analytical solution of the RTE allows avoiding this problem. The approximation is based on a slight decrease in the anisotropic part of the angular spectrum. The matrix-operator method determines the general solution for a complex multilayer structure. The calculation speed can be increased without compromising the accuracy of the solution with the help of the synthetic iterations method. The method consists of two stages: the first one repeats the described one with a small number of ordinates; on the second one the iteration of it is performed. The model is realised in the Matlab software. Keywords: radiative transfer equation, matrix-operator method, synthetic iterations method. z 1. Introduction    ( ) d  ; x(ˆl, ˆl) is indicatrix of scattering;  is single 0 Reflective properties of surfaces play the key role in light- scattering albedo. The problem (1) is defined in the coordinate ing calculations, and in most cases, calculations are impossible system OXYZ, the axis OZ is directed perpendicularly down, ẑ without knowing them. Some objects reflect both diffusely and is the unit vector along OZ. z = 0 on the upper boundary. specularly. Taking into account such surfaces is essential for Numerical solution implies discretization of the RTE. The correct calculations. Radiance factor is the classical value used scattering integral should be replaced by a finite sum [2]. It is to describe spatial reflective properties. In [14] a model of radi- impossible, since there are singularities in the angular radiance ance factor of a uniform layer with diffuse bottom and Fresnel distribution, which cannot be replaced by a finite series in any upper boundary was described. The present article is aimed to basis. The article [9] offers to single out the direct source radia- find the solution for a multi-layered structure and to optimize tion, which became an indispensable element of all RTE solv- calculations (more fast calculations without lowering in accura- ing methods. However, a feature of all natural objects is the cy). presence of suspended particles, which sizes significantly ex- ceed the wavelength, which, in accordance with Mie Theory, 2. Boundary value problem of the RTE for a leads to high angle scattering anisotropy, and the extraction of uniform layer direct radiation is not enough efficient. In case of presence of A partially coherent wave reflection from a structure with a strong anisotropy in the radiance angular distribution, replacing complex optical characteristics analysis shows [4] that the wave the integral with a finite sum can lead to significant uncertain- is dived into two components, after solution averaging: a coher- ties [15]. ent one determining the optical characteristics and a quasiuni- form admitting the beam description and obeys the RTE. If the 3. RTE discretization irradiating inhomogeneities of the medium are located in the The idea is of the accurate discretization is to analytically Fraunhofer region from each other, which is most often realized distinguish all the singularities and the anisotropic part La ( , ˆl ) , in practice, then the optical characteristics of the medium are determined by summing them over the elementary volume of i.e. to represent the solution in the form [1]: the medium. L(, ˆl )  La (, ˆl )  L(, ˆl ) , (2) Radiance factor is the relation between the radiance coming out from a surface in the specified direction and the radiance of here L(, ˆl ) is the regular part of the solution (RPS), which is a an ideal diffuse surface in the same conditions. Therefore, the smooth function. Such a smooth function can be represented by determination of the surface radiance factor is reduced to a a finite basis of elements. boundary value problem for a flat layer. Let us analyze the al- Considering (2) the boundary value problem (1) for L(, ˆl ) gorithm for the numerical solution of the RTE boundary value is [2]: problem for the case of irradiation of a turbid layer with optical thickness 0 by a plane unidirectional source (PU) in the direc-  L(, ˆl )   L(, ˆl )  x(ˆl , ˆl) L( , ˆl ) dˆl   ( , ˆl ); 4   0  0 0  tion ˆl  1  2 ,0,  ,   (zˆ , ˆl )  cos  : 0 0 0    L (, ˆl )  0; L( , ˆl )   La ( , ˆl ), (3)   0,   0  0 ,  0  L(, ˆl )  ˆ ˆ ˆˆ ˆ    L(, l )  4  L(, l) x( l, l)dl, ˆ where the source function (, l ) is defined by the anisotropic  (1)  L(, ˆl )  (ˆl  ˆl0 ), L(, ˆl )  0; part of the solution:   0,  0  0,  0 dL (, ˆl )  (, ˆl )   a  La (, ˆl )  x(ˆl, ˆl) La (, ˆl)dˆl . (4) where L(,,) is the radiance in the viewing direction d 4    ˆl  1  2 cos , 1  2 sin ,  ,   (zˆ , ˆl ) at optical depth As a numerical method implies finding an approximate val- ue, it cannot exactly correspond to boundary conditions. That is why the second boundary condition is changed. Copyright © 2019 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). Since L(, ˆl ) is a smooth function, it can be represented in tion (8) with respect to flows exiting the layer. The solution is get in the form of scatterers [6]: a finite basis. For example, if the RTE is discretized by the discrete ordinate method (DOM), the representation is as fol-  C (0)   F   R T   C (0)      , (9) C (0 )   F   T R  C (0 )  lows: M L(, i , )   (2  0, m )cos(m)Cm (, i ) , (5) F  R T   u12  e 0 u11  m0 where     hJ,    h   0 , T  F  T R e u 22  u 21  where C ()  C (, i ), C()  C (), C ()  m  is a column 1   u11 e 0 u12  vector of discrete values of the azimuthal expansion coefficients h  .   0 of radiance in a Forier series,  j  0,5( j  1) , j are zeros of   e u 21 u 22  Gaussian quadrature of the order N/2 for zenithal discretization The expression (9) in the form of scatterers gives us the re- of the scattering integral. Index m is further absent due to the lation between the streams emerging from the layer and the lack of need for it. incident ones and is a generalization of the radiance factor. The This allows to replace the scattering integral by a finite column F describes the inner radiation of the layer, the matri- sum. The boundary value problem transforms to a boundary ces R and T represent discrete values of the radiance factors value problem for a matrix inhomogeneous linear differential for reflection and transmission. equation of the first order with constant coefficients: d C() 5. Invariance of the solution   BC()  M 1 (), B  M 1 (1  AW) , (6) d In the simplest case, a multilayer structure can be represent- here wj are the weights of the Gaussian quadrature of the zen- ed by only two layers: ithal order N/2, Pln () are associated polynomials Legendre, C1   F1   R1 T1  C1  Pl0 ()  Pl () are polynomials Legendre,  1    1     1 , C   F   T1+ R1+  C    wi 0   0   K (10) W  0 w , M    , A   (2k  1) Pkm (i ) xk Pkm ( j ), i C2   F2   R 2  T2   C2  4 i  0 i  k 0  2    2     2 , K C   F   T2  R 2   C  x()   (2k  1) xk Pk (cos ) . where the subscript defines whether the upper 1 or lower 2 be- k 0 longs to. Vertical arrows indicate the direction of radiation inci- dence on the layer. 4. Propagator and scatterer The radiation emerging from the first layer is the radiation entering the second layer and vice versa, that is: The solution of the matrix equation (6) can be presented [5] as the sum of the general solution of the homogeneous equation C1  C2  C , C2  C1  C . (11) and the particular solution of the inhomogeneous: Solving the system with respect to the reflected radiation 0  C(0)  P(0, 0 )C(0 )   P(0, ) M 1 (, 0 )d  , (7) C  and transmitted radiation C 2 with respect to the incident 1 0 radiations on the system, we obtain: where the function P(t , )  eB(  t ) is the solution of the homo- 1   C1   F  T1  R 2  F  F  1 2   2      geneous equation describing the relationship of the radiance distributions at two points t and  of the medium without inter- C  T2   F1  R 1 F2  F2    (12) nal sources. Such a function is called propagator.  R  T  R 2  T1 T1  T2   C1  There are problems with the solving of (7), because the   1 1   , propagator contains both positive and negative exponent func-  T2   T1 R 2   T2   R 1 T2   C2    . tions. This physically corresponds to flows propagating top-to- 1 bottom and bottom-top. Thus, the matrix is poorly conditioned, where   1  R 2  R1 and the calculations for fields with >1 become impossible. The expression (12) for two adjacent layers is completely Scale conversion is proposed to avoid this effect [7]: equivalent to the expression for one layer (9) in the form of 0 scatterers, but with effective parameters (coefficient) that can be  SU 1 C(0)  HU 1 C(0 )  J, J  S  et U 1M 1 (t ) dt , (8) obtained from the parameters of each layer separately. This 0 shows the invariance property of the solution of scatterers, where U is the eigenvector matrix of the matrix B ; which allows us to logically go over to the invariant immersion   diag(  ,   ) is the eigenvalues matrix,      ; of V.A. Ambartsymian [13]. On the other hand, invariance allows the calculation of a u u  1 0  e 0  0 layer with an arbitrary vertical inhomogeneity, breaking it into U 1   11 12  , S   , H   . u 21 u 22  0 e    0   0 1 an arbitrary number of homogeneous layers. In this case, two adjacent layers can be replaced by one layer described by ex- The boundary value problem (1) and the corresponding sys- pression (12). This approach in transport theory is called the tem (6) are two-point problems, as the conditions on the bound- matrix operator method [10] (MOM). The advantage of the aries are given, and the solution inside the layer needs to be obtained expression (12) is the allocation of the anisotropic part found. Therefore, the solution (7) based on propagators is not of the solution in arbitrary form (2). complete [11]. Column vectors C  (0), C  ( 0 ) in (8) describe the flows 6. Anisotropic part of the solution falling on the layer and are defined by the boundary conditions. If we talk about the complexity of calculations and, accord- The vectors C  (0), C  ( 0 ) correspond to the radiation flows ingly, the speed of calculations using some software, the key reflected from and transmitted the layer. Let us solve the equa- role for this has the sizes of the matrices included in the matrix solution (9) and (12), i.e. the values of the constants N, M, K, where the index 1 related to the layer is omitted due to the ab- where N is the number of discrete ordinates, M is the number of sence of need. All other azimuthal harmonics m>0 are deter- azimuthal harmonics, K is the number of members of the de- mined by the expression of a single-layer medium (9). composition of the indicatrix into polynomials Legendre. This approach does not work on the boundary with refrac- In the general case of an arbitrary incidence angle, the val- tion, since the directions of the rays vary according to the ues are approximately equal: M ≈ N ≈ K [2]. However, with a Snell’s law: successful choice of the anisotropic part, it can be achieved that n1 sin 1  n2 sin 2 , (19) the angular dependence of the RPS will be close to isotropic. where n1, n2 are refractive indices of the medium, and the corre- Thus, the calculation speed can be significantly increased be- spondence of ordinate directions is violated. A solution to this cause K >> N >> M. problem was proposed in [12]. The isotropic part was distinguished. However, it should be Let us consider in more detail the refraction on the practi- defined too. How can it be defined? The unequivocal answer in cally important case when n1 < n2. The first medium is called the spatial-angular representation is difficult, but the task is atmosphere for certainty (index a), na=1, and the second media greatly simplified for the spectral representation of the angular is the ocean (index o) with no>1. The cosines of the rays with distribution. The more anisotropic the angular distribution is, the axis OZ in both media according to (19) will be related to the more smooth is its spectrum. The radiance distribution can each other by the expression: be represented as polynomial Legendre series:  a  1  no2 (1  o2 ) . (20)  2k  1 La (, )   Z k () Pk () , (13) An important feature of the problem and its further solution k 1 4 is the presence of a region of total internal reflection. It is seen where the assumption is made that the anisotropy in the region of small angles is much stronger than the azimuthal asymmetry from (20) that this region appears for o  t  1  1 no2 im [2],   (ˆl , ˆl0 ) . the ocean medium: the rays do not exit the ocean, but are ideal- The spectrum Zk() of the anisotropic part of the solution ly reflected again into the ocean. Boundary conditions for the slowly monotonically decreases from the index k. This allows total internal reflection region are formulated without problems. Then the scattering integral can be represented as the sum of us entering a continuous function Z(k,). Since the function three integrals taking into account the total internal reflection slowly and monotonously decreases, the following expansion in region: Taylor series is valid: 1 t Z ( k, )  Q ()C (, )d    Q ()C (, )d m m m m Z (k  1, ) Z k ()  . (14) k k k k k 1 1 (21) t 1 Substitution of (14) into an infinite system of differential   Q ()C (, )d    Q ()C (, )d . m k m m k m equations for solving the RTE by the spherical harmonics t t method leads [3] to one equation of mathematical physics that The first and last integrals are related to the region of re- allows an analytical solution fraction, and the second relates to the total internal reflection Z k ()  exp  (1  xk )   0  , (15) region. For the second integral, we perform the transformation: which is called the small-angle modification of the spherical t 1  Q ()C (, )d   Q ()C (, )d ,     , (22) m m m m harmonics method (SHM). k k t In [8], a comparison was performed of the solution (9) with t 1 the separation of the anisotropic part based on the SHM (the which makes it possible to apply a double Gaussian quadrature MDOM program) with the main known programs: MDOM with with Nt nodes and subsequently move in this zone to two the same accuracy in calculating the angular distribution of streams of ordinates Ct , Ct , which are connected by an ideal radiance exceeds other programs by 1-2 orders of magnitude in mirror reflection at the boundary. computational speed. For the first and last integrals in (21), we make the trans- 7. Reflection and refraction at the interface of formation of the integration variable to a by expression (20): two media o  1  (1  a2 ) no2 , d o   a d  a no2  (1   a2 ) . (23) Let us consider a special case of a layer, in which the lower It is obvious that in the transition to discrete ordinates a boundary reflects according to Lambert’s law and has a reflec- complete correspondence is established between atmospheric tance In this case, the boundary condition in (1) on the bot- ordinates С a , С a and ocean ordinates Сo , С o . If we combine tom is the following: the vectors into one in accordance with the rules of the Matlab  L(0 , ˆl ) L(0 , ˆl)dˆl .  t  ocn  o t   (0)    С ; С  , С  С ; С  , then all the rela- software Сocn o (16)  0 tion in MOM will be valid for them. The introduced values also The following matrix expressions can be obtained by substi- allow us to write down the condition at the atmosphere-ocean tuting the azimuthal representation of radiance (5) and discrete boundary: ordinates, and the integral (16) by Gaussian quadrature:  Ca   R Tao   Ca  m  0 : C  (0 )  2 R LC  (0 ); m  0 : C  (0 )  0 , (17)  ocn      ocn  , (24) C  Toa R oo  C  where the matrix of Lambert reflection R L consists of the same 0  0 1   N/2 lines  j w j . where Tao  [T 0], Toa    , R oo    , R, T are Fresnel   T  R 0 In accordance with MOM (12), we obtain the following ex- reflective matrices. pression for the reflected component of zero harmonic:   RF, 1 C (0)  F  2 T 1  2 R LR L  (18) 8. Synthetic iterations method of calculations without loss of quality. The number of discrete ordinates N determines the size of the matrices with which the Fig. 1 shows a comparison of the calculations in MDOM of calculations are performed. That is, a decrease in N would radiance angular distribution reflected by a layer for different speed up the calculations. Thus, the main question is: how not sets of parameters N and M. It is easy to see that the calculation to lose in quality? The synthetic iterations method is the answer of the almost complete distribution of the viewing angle is to this question. much faster than the calculation of individual small sharp The synthetic iterations (SI) method was proposed in nucle- peaks. The difference in calculation time t is more than 150 ar physics [1]. In this case, the iteration splits into two stages. times. What is the reason? At the first stage, an approximate solution is sought that con- The angular distribution of the RPS is actually close to iso- verges well in the average energy metric, and at the second, the tropic, but with some small ripples. A fair question arises: how usual iteration is performed, which significantly increases the many discrete ordinates N are necessary to represent this small convergence in a uniform angle metric. Since the developed ripple? Since the RPS is a smooth function, its expansion into a method for solving MDOM has good convergence in the aver- series of polynomials Legendre has a finite number N: age metric, we should count on its significant increase in con- N 2k  1 m vergence after iteration. Lm (, )   Lk Pk () . (25) A numerical comparison of the reflected radiance in the k 1 2 All polynomials Legendre of the order < N can be expressed first iteration with MDOM is presented in Fig. 2, where t is through an N+1 polynomial order: the computation time. N 1 PN 1 () Pk ()   Pk (i ) , (26) i 1 (  i ) PN 1 (i ) where i are the roots of the polynomial PN+1(). Accordingly, this leads to the expression: N 1 PN 1 () N 2k  1 m Lm (, )    Lk Pk (i ) i 1 (   i ) P N 1 ( ) i k 0 2 or N 1 PN 1 () Lm (, )   Lm (, i ) , (27) i 1 (  i ) PN 1 (i ) which corresponds to the Lagrange interpolation formula for the function L(,). a) full range of viewing angles Fig. 1. Comparison of the reflected radiance angular distribu- tions. Direction to Nadir = 180о. Solid line: N=801, K=299, M=256, Δt=52.3 s. Dashed line: N=101, K=299, M=8, Δt=0.32 s. The latter relations are the analogue of the Nyquist-Shannon theorem on samples for the angular spectrum of the angular distribution over polynomials Legendre. Hence: b) range in the vicinity of gloria, the designations are the same 1. MDOM method provides average convergence; as in Fig. 2a 2. All methods for isolating the anisotropic part are equivalent Fig. 2. Comparison of synthetic iteration (SI) with MDOM for to each other in a uniform metric; the radiance of the reflected radiation of a layer 3. To achieve good convergence in a uniform metric, the sam- pling interval should correspond to the angular size of the It is seen from the figure that the greatest difficulty for cal- smallest part, which should be reproduced on the radiance culating in MDOM is the region near gloria (Fig. 2b). To calcu- distribution. late this region in the MDOM program, N = 801 and M = 256 When implementing a multilayer surface model taking into are required, which corresponds to a sampling step of less than account diffuse reflection from the lower boundary and reflec- 0.5°. To achieve the same accuracy within the framework of tion from the upper boundary according to the Fresnel law in synthetic iteration, only N = 11, M = 4 is necessary, which re- the Matlab software, the question arose of possible acceleration duces the computation time by almost 60 times. Accordingly, the synthetic iteration from MDOM allows us to calculate the [Model of a scattering layer with a diffuse bottom and a angular distribution of radiance with an accuracy in the uniform Fresnel boundary] // GraphiCon 2018. Conference pro- metric of no worse than 1% at a counting time of no more than ceedings. Tomsk, Russia. P. 399-401 (in Russian). 1 second. [15] Krylov V.I. Priblizhennoe vychislenie integralov [Ap- This method allowed to significantly increase the speed of proximate calculation of integrals] // Nauka Publ., Mos- calculations carried out in the Matlab software using the created cow, 1967 (in Russian). model. 9. Conclusion The mathematical model of the luminance factor for a mul- tilayer medium bounded above by the Fresnel and below Lam- bert surfaces is implemented. The calculation is optimized in terms of speed and accuracy of calculation. The obtained de- pendencies and characteristics qualitatively coincide with the expected ones. The model must be filled with parameters of real media, for which it is planned to experimentally test the model. In the future, it is also planned to take into account reflection from a randomly uneven border and polarization. 10. References [1] Adams M.L., Larsen E.W. Fast iterative methods for dis- crete-ordinates particle transport calculations // Progress in Nuclear Energy, 2002. Vol.40, No.1. P.3-159. [2] Budak V.P., Klyuykov D.A., Korkin S.V. Convergence acceleration of radiative transfer equation solution at strongly anisotropic scattering // In Light Scattering Re- views 5. Single Light Scattering and Radiative Transfer / Ed. A.A. Kokhanovsky. Springer Praxis Books, 2010. P.147-204. [3] Budak V.P., Korkin S.V. On the solution of a vectorial radiative transfer equation in an arbitrary three- dimensional turbid medium with anisotropic scattering // J. Quant. Spectrosc. Radial. Transfer., 2008. Vol. 109. P. 220–234. [4] Budak V.P., Veklenko B.A. Boson peak, flickering noise, backscattering processes and radiative transfer in random media // J. Quant. Spectrosc. Radial. Transfer., 2011. Vol. 112. P.864-875. [5] Flatau P.J., Stephens G.L. On the Fundamental Solution of the Radiative Transfer Equation // JGR, 1988. V.93, No.D9. P.11,037-11,050. [6] Fryer G.J., Frazer L.N. Seismic waves in stratified aniso- tropic media // Geophys. J. R. Astr. Soc., 1984. V.78, P.691-698. [7] Karp A.H., Greenstadt J., Fillmore J.A. Radiative transfer through an arbitrarily thick, scattering atmosphere // J. Quant. Spectrosc. Radial. Transfer., 1980. Vol. 24. P. 391-406. [8] Kokhanovsky A.A. et al. Benchmark results in vector at- mospheric radiative transfer // J. Quant. Spectrosc. Radi- al. Transfer., 2010. Vol.111. P.1931-1946. [9] Milne E.A. The reflection effect of the eclipse binaries // Mon. Not. Roy. Astrophys. Soc., 1926. Vol. LXXXVII. P.43-49. [10] Plass G.N., Kattawar G.W., Catchings F.E. Matrix Opera- tor Theory of Radiative Transfer // Appl. Opt., 1973. Vol.12. P.314-326. [11] Redheffer R.M. Inequalities for a matrix Riccati Equation // Journal of Math. and Mechanics, 1959. V.8. P.349–367 [12] Tanaka M., Nakajima T. Effects of oceanic turbidity and index refraction of hydrosols on the flux of solar radiation in the atmosphere-ocean system // J. Quant. Spectrosc. Radial. Transfer., 1977. Vol.18, No1. P.93-111 [13] Ambartsumian V.A. K zadache o diffuznom otrazhenii sveta [To the problem of diffuse reflection of light] // JETF, 1943. Vol. 132, No. 9-10, P.323-334 (in Russian). [14] Basov A.Yu., Budak V.P. Model' rasseivayushchego sloya s diffuznoj podlozhkoj i frenelevskoj granicej