Matrix riccati equation solution of the 1d radiative transfer equation
JJournal of Computational and Theoretical Transport MATRIX RICCATI EQUATION SOLUTION OF THE 1D RADIATIVE TRANSFER EQUATION
Barry D. Ganapol
Department of Aerospace and Mechanical Engineering University of Arizona Tucson, AZ, 85721 [email protected]
Japan Patel
Department of Mechanical and Aerospace Engineering, The Ohio State University Columbus, OH, 43210 [email protected]
ABSTRACT
In recent years, the first author has developed three successful numerical methods to solve the 1D radiative transport equation yielding highly precise benchmarks. The second author has shown a keen interest in novel solution methodologies and an ability for their implementation. Here, we combine talents to generate yet another high precision solution, the Matrix Riccati Equation Method (MREM). MREM features the solution to two of the four matrix Riccati ODEs that arise from the interaction principle of particle transport. Through interaction coefficients, the interaction principle describes how particles reflect from- and transmit through- a single slab. On combination with Taylor series and doubling, a high quality numerical benchmark, to nearly seven places, is established.
Key Words : Interaction principle, Doubling, Cloud C1 and HAZE L phase functions, Benchmarks.
1. INTRODUCTION
While the 1D radiative transfer equation is probably the most solved transport equation with numerous numerical and analytical solutions available, few solutions can claim to be an extreme benchmark. In context, an extreme benchmark means high confidence in six to seven places (seven to eight significant figures) of the numerical solution for angular intensity. In addition, 1D means exclusively a homogeneous plane parallel medium with a highly anisotropic scattering phase function represented either by a Legendre expansion or point-to-point discrete directions. The radiative transfer equation in question will simulate azimuthally integrated unpolarized radiation only. Thus, when referenced, the transport equation is for monochromatic- 1D- plane parallel- azithmuthally integrated- anistropically scattering- neutrons or unpolarized photons in a homogeneous medium. In addition, we limit our investigation to a perpendicular beam source impinging on the left (near) surface of a conservative ( ω = 1) slab medium in vacuum. Our intent is to develop a radiative transfer benchmark to extreme precision for anisotropically scattering slabs of small to modest optical thicknesses. ournal of Computational and Theoretical Transport Our approach commonly falls under the category of the matrix-exponential method, which includes the discrete ordinates method, method of matrix operators and the matrix Riccati equations. These methods all share discretization of particle directions. The solutions include eigenvalue/eigenvector expansion, spatial finite difference and iterative sweeping in particle directions as well as doubling and adding via the interaction principle. For a through summary of the available literature, we refer the reader to the excellent work of Efremenko, et. al. [1].
The interaction principle is the basic transport algorithm that defines Discrete Space Theory [2], which gives the solution to the transport equation with discretization only in direction µ . The spatial variable τ remains continuous. The principle also goes under the name of invariant imbedding and is one of the first analytical methods developed to solve the radiative transfer equation [3]. The interaction principle is an alternative form of the discrete ordinate or S N equations avoiding the ad-hoc diamond difference approximation commonly found with spatial differencing. When combined with doubling, the interaction principle leads to a fast and precise numerical method. Recognizing that S N directional discretization is consistent, i.e ., as the number of discretizations become infinite, the exact solution results; offering an opportunity for convergence acceleration. Previously, convergence acceleration and doubling has been considered in a limited instance [4]. Here, convergence acceleration, through the Wynn-epsilon algorithm, enhances numerical efficiency of the interaction principle. The matrix Riccati equations arise from the solution of the neutral particle transport equation in discrete ordinates form ( ) ( ) ( ) Nm m m m m mm
I N f I N µ τ µ α µ µ τ µ τ τ ττ ′ ′ ′′= ∂ + = < < ∂ ∑ (1a) for directions m µ ,
1, 2, ..., 2 m N = . Thus, µ m replaces µ in the angular intensity ( ) τ µ , I , where µ m takes on both positive and negative values , 1,..., , 0. µ µµ µ + − = > ≡ m mN m m m N (1b) An N th order double Gauss quadrature approximates the collision integral over the half-range intervals [-1,0], [0,1] as ournal of Computational and Theoretical Transport ( ) ( ) ( ) ( ) ( )
11 2 1 , , , , ; , ; , µ µ µ τ µα µ µ τ µ τ µ − ′ ′ ′′= ′ ′ ′ == + ∫ ∑
N m m m Gm d f If I N E N , (1c) where ( ) , ; τ µ G E N is the quadrature error [5] ( ) ( )
210 4 ,, ; 2
NG N
MeE N N τ µτ µ π . (1d) For an integrand with at least 4 N continuously bounded derivatives, the Gauss quadrature converges to its integral and therefore one can show the exact solution, ( ) ( ) ( ) , , ; , ; I I N e N τ µ τ µ τ µ= + , with error ( ) , ; e N τ µ , is the limit ( ) ( ) , lim , ; N I I N τ µ τ µ →∞ = . (1e) The limit plays a prominent role in the application of convergence acceleration applied below, but before doing so, we must find the approximate intensity ( ) , ; I N τ µ . To form Eqs(1) as a vector equation, one defines the vectors and matrices (suppressing N ) ( ) { } ( ) { } ( ) { } ( ) { } ( ) { } { } .... , ; , 1,..., , NN N N
Tm m m mm
I I Idiag f j m Ndiag τ τ τ τµ α µ µα + + −− −+±± ′+− ++ ≡ ≡ ≡ = ± ± = ≡ IM C CC W WC CW (2a-d) where ( ) ( ) , ; m m
I N I τ µ τ ± ± ≡ (3a) ournal of Computational and Theoretical Transport to give ( ) ( ) ( ) ( ) . N d d τ τ ττ ± − ±± ± − ± = − ± I M I C W I M C WI (3b) I N is the size N matrix identity, not to be confused with intensity. Therefore, by stacking the vectors in directions ± µ m ( ) ( )( ) ττ τ −+ ≡ II I , we have ( ) ( ) d d τ ττ =
I AI , (4a) with . − ≡ − A α ββ α (4b) Since from scattering invariance, +− −+ = C C , there results ( ) . N − ++− +− ≡ −≡ I C WC W α Μβ Μ (4c,d) Equations (1) and (4a) are subject to the incoming intensity at τ = τ and no source incoming intensity at τ = τ ( ) ( ) , τ τ + − ≡ ≡ I g I (1f) to complete the problem statement. g is the angular intensity for a perpendicular incoming beam on the near surface τ specified later. Since we consider a conservative medium, there is no re- radiation making the medium source less. ournal of Computational and Theoretical Transport At this point in our development, there are several paths forward to achieve a benchmarking goal. For your interest, APPENDIX A includes brief descriptions of three successful approaches. All methods referenced generate benchmarks to seven places to both the HAZE L 82 and Cloud C1 300 term scattering kernels [6]. We now give the details of the Matrix Riccati Equation Method (MREM) of solution featuring doubling, and offering an entirely different formalism than previous solutions.
Consider a slab layer of thickness τ τ− as shown in Fig. 1. Rather than seek an analytical exponential eigenvalue/eigenfunction based solution to Eq(4a), equally appropriate since A is independent of τ , we derive the solution ( ) , ττ P t o the corresponding fundamental matrix equation ( ) ( ) ( ) , , ; , N d d ττ ττ τ ττ = = P AP P I (5a) by letting ( ) ( ) ( ) ; τ ττ τ= I P I . (5b) ( ) , ττ P is also called the transition matrix and remains valid even when the A matrix Fig. 1. Single Slab. depends upon τ . Therefore, Eqs(5) are the starting point for consideration of variable scattering across the slab (a future effort). Recalling that the superscript of the intensity indicates direction of particle motion, we write the stacked vectors at the boundaries shown in Fig. 1 as τ τ ournal of Computational and Theoretical Transport ( ) ( )( ) ( ) ( )( )
00 0 . ττ τττ τ −+−+ = = II III I (6a,b) The interaction principle says a homogeneous slab will reflect and transmit particles entering either face equally. Thus, if R (τ−τ ) and T (τ−τ ) are the reflectance and transmittance matrix coefficients of the homogeneous medium depending on the medium thickness and scattering and absorbing properties, we can the write from physical reasoning ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) , τ τ τ τ τ τ ττ τ τ τ τ τ τ − + −+ + − = − + −= − + − I R I T II T I R I (7a) or in vector form ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) , τ τ τ τ ττ τ τ τ τ ττ −− ++ − −= − − T R II R T II (7b) which is the interaction principle. Thus, the exiting intensity at surface τ results from reflection from the incoming at τ and that transmitted from the far surface at τ . Likewise, the intensity exiting the surface at τ results from that reflected from the incoming at τ and that transmitted from the near surface at τ . The interaction principle separates the slab influence on the particles within the slab into reflection and transmission. The RHS is the source of particles into the slab and the LHS is the consequence of particle interaction within the slab, or output response. Hence, we call the matrix connecting the output to the input the response matrix ( ) ( ) ( )( ) ( ) τ τ τ ττ τ τ τ τ τ − −− = − − T RQ R T . (8a) to give the interaction principle for a homogeneous slab as ( ) ( ) ( ) ( ) ( ) . ττ τ τ ττ −− ++ = − II Q II (8b) By construction ournal of Computational and Theoretical Transport ( )( )
0 0 N . == RT I (8c) The response matrix is block symmetric and depends upon slab thickness and material scattering and absorption properties. Henceforth, to simplify without loss of generality, we let τ = 0; and therefore in the context used here, τ is the slab thickness. In addition, when obvious, we suppress the interaction coefficient dependence on slab thickness. The interaction principle is an exact formulation of particle transport, as much as the radiative transfer equation is, once we know the reflectance and transmittance, which now becomes our task. Re-arranging the vector equation Eq(8b) gives ( ) ( ) τ = NN T I RI IR I T −− . (9a) or ( ) ( ) τ − = NN T I RI IR I T −− . (9b) From Schur’s complement N N − − − T T=R I RT I − , and therefore ( ) ( ) τ − −− − = T T RI IRT T RT R −− . (10) Then, from Eq(5b) ( ) τ − −− − = T T RP RT T RT R −− . (11) ournal of Computational and Theoretical Transport Finally, the following four Matrix Riccati differential equations result for T and R directly from Eq(5a): ( ) , d d dd d dd d d dd d d d τ τ ττ τ τ τ − − −− − −− −− − −− − − = − + + T T RR T T T RRT T RT RR T TT R RT R − α β −β α −− which give by equating corresponding block matrix components (after some algebra) dd τ R = T T β dd τ T = R T T β −α (12a-d) + dd τ R = R R R R β− α−α β dd τ T = T R T β − α . In addition, there are two obvious identities + T T = R R R R β β− α−α β (12e,f)
R T T = T R T β −α β − α . Next, we need to determine analytical expressions for interaction coefficients R and T .
2. TAYLOR SERIES REPRESENTATIONS FOR R AND T At this point, there are numerous ways to numerically evaluate the interaction coefficients. Possibilities include the multiple collision technique [7], where one partitions intensity into particle collisions whose sum is the intensity at τ and µ . This approach becomes analytically unwieldly with increasing number of collisions [1]. Another method forms integral equations for Eqs(12) and proceeds iteratively [8], which is essentially the Taylor series (TS) formulation. Also, integration by Runge-Kutta methods have been applied to this non-linear set of ODEs [9]. Here, we will follow the TS representation in two distinct forms, direct and indirect. The direct form comes from Eqs(12c,d), where the interaction coefficients have a TS representation in τ even though equations are non-linear. The indirect form comes from a product transform rendering Eqs(12) linear in each factor of the product. A TS then applies to each factor with the resulting product also a TS. ournal of Computational and Theoretical Transport Even though the Matrix Riccati equations are quadratically nonlinear, analytical solutions, possibly not in closed form, are possible through Taylor series.
We now focus on solving the matrix Riccati equation Eq(12c) for the reflectance R since it is uncoupled. The most convenient solution is the TS formulation ( ) nnn τ τ ∞= = ∑ R R , (13a) where ( ) ( ) ( ) nn n τ τ τ = ≡ R R . Then, applying n -1 derivatives to Eq(12c) and dividing by n !, results in ( ) ( ) ( ) ( ) ( )
1+ ! τ δ τ τ τ ττ −− − − nn n n n n dn n d
R = R R R R β − α−α β . (13b) The last term comes from Leibnitz’s differentiation rule for scalars and vectors alike ( ) ( ) ( ) ( ) τ τ τ ττ − − − −− = = ∑ n n n l ln l d f g f gn d (14) giving the recurrence for τ = + nn n n n n l ll n δ −− − − − −= = − ∑ R R R R R β α−α β . (15)
Similarly, if ( ) τ τ ∞= = ∑ nnn T T , (16a) ournal of Computational and Theoretical Transport by applying n -1 derivatives to Eq(12d), dividing by n ! and by application of Leibnitz’s differentiation, ( ) ( ) ( ) ( ) ( ) ( )( )
1! ! 1 ! n nnn dn n d n n τ ττ ττ −−− − T T= T R β − α , or nn n l l nl n − − − −= = − ∑ T T R T β α , (16b) where ( ) ( ) nn n τ τ = ≡ T T . While the direct TS is certainly numerically adequate, it is of analytical interest to derive an alternative.
One can linearize Eqs(12) through the following transformations: ( ) ( ) ( ) τ τ τ − ≡ R Q P (B1) − T = WP (B14) (Do not to confuse Q with the response matrix). As one shows in APPENDIX B, W is unity and the corresponding TS coefficients are nn n n l ll −= − ∑ R = Q P R . (B12)
10 0 nn N n n l ll δ − −= = − ∑ T I P T . (B19) Numerically, the direct and indirect representations give identical results, but linearization is important for the theoretical connection between formulations.
3. NUMERICAL EVALUATION
We are now in position to establish numerical implementation. Built around Eq(8b) with τ = 0 and thickness τ L expressed as ournal of Computational and Theoretical Transport ( )( ) ( ) ( )( ) τττ − −+ + = LLL
I IQI I , (17) our strategy is to first determine the interaction coefficients from their TS representations of Eqs(13a) and (16a). We follow this by application of the method of doubling to find the overall slab response matrix to give the outgoing intensity from the incoming intensity specified at slab boundaries according to Eq(17). In this way, one finds the exiting angular intensity distributions. The interior angular intensity distribution at τ comes from the method of adding responses of two complimentary slabs of thickness τ and τ L - τ . In doing so, we find the intensity distribution at any point within the slab without requiring the unknown exiting intensity distributions as a pre-calculation. We believe this procedure to be knew. As will be discussed, each numerical procedure must include additional features to make MREM worthy of a benchmark. The method of doubling begins with the interaction principle for a slab of thickness τ l -1 in block matrix form l ll ll l ll − − −− −− −−+ + +− −− = = Q QI I IQ Q QI I I . (18) (Note: We suppress the argument in the response matrix.) Figure 2 shows the slab, called slab l -1 with intensities at boundaries at and τ l -1 . Equation (18) takes into account block symmetry of the reflectance and transmittance coefficients , − − − − − − = = = = l l l l l l Q Q R Q Q T . Fig. 2. Single slab of thickness τ l -1 . Note that a single slab can be a combination of many slabs whose equivalent response matrix comes from the following procedure.
Consider a slab of thickness τ l - τ l -1 with different properties and interaction principle τ l -1 ournal of Computational and Theoretical Transport − − −−+ + +− − = = l ll l ll l ll l l T RI I IR TI I I R joined to the first slab as shown in Fig. 3. Writing out the interaction principle , , ; − − +− − −+ − + − − +− − − − − − + − +− = += + = + = + l l ll l l l l l l l l l l l l l I Q I Q II Q I Q I I T I R I
I R I T I (19) and eliminating the intensities at the interface gives the combined response l Q for the composite slab [0, τ l ] and the exiting intensities ( ) τ − −+ + = l lll I IQI I , (20a)
Fig. 3. Two Slabs. from block recurrence ( ) ,1 ,2 1,1 1,2,3 ,4 τ − − ≡ = + l l l ll ll ll l l Q Q Q QQ UQ Q T R , (20b) with
N ll l Nll l −− −+ − ≡ − ≡
I Qw R IQw T (20c,d,e) ( ) , τ τ τ −− − + − l l l l ll U = w w . τ l -1 τ l ournal of Computational and Theoretical Transport Finally, from the interaction principle in terms of slab l interaction coefficients, ( ) τ − ≡ l ll l l l T RQ R T . Thus, the response for joined slabs of any thicknesses and properties results from recurrence for each of the four square matrix partitions of the response Q l . While not initially apparent, the composite response is also block symmetric. In addition, in deriving Eq(20b), an intermediate step gives the intensity at the interface [10] ( ) , τ τ τ + −− −− +− = − l ll l lll I IUI I (21) to be used below. We are ready to double. Doubling arises from the definition of P ( h ) for a single slab (designated 1) from Eq(5b) with τ = 0 of thickness h ( ) ( ) ( ) = h h I P I (22) relating incoming to outgoing intensities. ( ) h P acts as a propagator, propagating particles back and forth through the slab. For a composite of two identical slabs, each with interaction coefficients R ( h )and T ( h ), propagation represents a joint probably of particle transfer between sequential surfaces of the composite slabs, with the exiting intensities ( ) ( ) ( ) = h h h I P I (23a) or from Eq(22) ( ) ( ) ( ) ( ) = h h h I P P I . (23b) Since ( ) ( ) ( ) = h h I P I , Eq(23b) specifies ( ) ( )
22 1 ≡ h h P P . (24) From Eq(11), for the single slab ournal of Computational and Theoretical Transport ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) − −− − = h h hh h h h hh h T T RP R T T R T R −− , (25a) and from Eqs (11), (23c) and (24) ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( )( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( )
2 .2 2 22 2 2 2 2 2 − −− −− −− − ≡ = h h hh h h h h hh h hh h h h h hh
T T RP R T T R T RT T RR T T R T R −−−− (25b) Equating the four matrix partitions across the equal sign results in the following well-known propagated interaction coefficients for a two identical slab composite [11]:
122 1 1 1 122 1 1 1 1 1 , − − = − = + − N N
T T I R TR R T R I R T (26a) where ( )( )( )( )
22 . ≡≡≡≡ h hh h
T TT TR RR R (26b) The identical relation should also be true from Eq(20b). With the determination of the interaction coefficients, the corresponding composite response for two identical slabs follows as ≡
T RQ R T . (26c) We next double the two-slab response to construct a four-slab response and then an eight-slab response etc., and continue until we construct the response matrix for the entire slab of thickness τ . We achieve this by initially dividing the slab into 2 n slabs each of thickness n n h τ≡ . (27a) ournal of Computational and Theoretical Transport Doubling to thickness 2 l h n , for l =1,…, n , the interaction coefficients from Eq(26a) are
122 2 2 2 122 2 2 2 2 2 − − −− − − − − − − = − = + − l l l ll l l l l l
N N
T T I R TR R T R I R T (27b,c) with coverage when l = n n n n nnn nn − − −+ + + = T RI I I= QR TI I I , (28a) or ( )( ) ( ) ( )( ) ( ) ( )( ) ( ) ( )( ) τ τ τ τττ ττ − − −+ + + =
T RI I I= QR TI I I . (28b) The significant numerical advantage of doubling is immediately apparent since rather than marching discretely from the near to the far surface in 2 n steps, one doubles to the far surface in n steps greatly reducing the propagation error. Note that after Eq(25b) ( ) τ P is not required in this formulation, though in other formulations it is [10]. There are several choices to find the interior intensity. The first and least advantageous is from the interaction principle itself. From Eq(5b) ( ) ( ) ( ) τ τ= I P I when applied recursively over l doublings to τ = l l h to give ( ) ( ) ( ) τ = l l l h I P I , (29a) where l P is constructed from ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) − −− − ≡ j j jj j j j j j j j j jl l lj l j j j j j jl l l l l l h h hh h h h h h h T T RP R T T R T R −− (29b) ournal of Computational and Theoretical Transport for j = 1,… l . Thus, the interior intensity distribution is ,1 ,2 0,3 ,4 0 − −+ + = l ll l ll P PI IP PI I . (29c)
The disadvantage here is that the exiting intensity leaving the near face − I must be determined first. In addition, the matrix inversion of the ill-conditioned transmission coefficient is required. More straightforwardly, from Eq(21) becomes, as referenced to Fig. 4 ( )( ) ( ) , τ ττ ττ + −− + = − LL I IUI I . (30) Now only the prescribed incoming intensities are required without any matrix inversion. To apply Eq(30), one needs to find ( ) τ Q , ( ) τ τ− L R by doubling and ( ) , ττ τ− L U from Eqs(20e). Fig. 4. For interior intensity at τ. The fundamental numerical algorithm depends upon the evaluation of the interaction coefficients of reflectance and transmission by TS ( ) nnn h τ ∞= = ∑ F F . (31) F is either of the quantities R or T . The coefficients F n are determined recursively for any n as shown above and in APPENDIX B. Each TS is for small slab of width h to begin doubling as discussed below and nothing more. For this reason, we expect rapid convergence. The series sum continues until the last value of all components of the matrix F n are within a prescribed relative error to the total sum. τ τ L ournal of Computational and Theoretical Transport Several relatively unfamiliar analytical features aid in the implementation of MREM, including faux quadrature interpolation, exact beam source and the grazing angle.
As will be shown, a key element of the benchmark from MREM will be faux (false) quadrature interpolation. Interpolation based on the matrix Riccati equations themselves, proceeds as follows. One can add additional (faux) quadrature abscissae to the quadrature list provided they have zero weight. Thus, the solution algorithm gives interpolated intensities at the faux abscissae and their addition does not modify the solution in any way. Faux quadrature can highlight a particular direction range or insure a particular set of directions is in each grid approximation for error assessment.
Generally, beam source implementation in discrete ordinates methods is thought to be limited by the quadrature set chosen. However, a beam along a specific direction is easily implemented using faux quadrature when one includes the beam direction in the quadrature list and sets the beam intensity to S / ε in that direction and zero for all other directions. The beam quadrature weight is then set to ε so that integration over all directions normalizes the beam source to S . Thus, the source is included in its actual direction. In this way, we have avoided the need to do a separate determination of the uncollided contribution. The intensity approximation is insensitive to the value of ε as long as it is small ( ε < 10 -16 ). µ = 0 The chosen quadrature cannot accommodate µ = 0 either directly or through faux quadrature interpolation. Thus, µ = 0 requires special treatment. In particular, from Eq(1a) with µ m = 0, we have (with argument N suppressed) ( ) ( ) ( ) , 0 0, , , 0 . τ α µ τ µ τ τ ′ ′ ′′= = < < ∑ N m m mm
I f I (32a) or in terms of the directional intensities at τ ( ) ( ) ( ) ( ) ( ) , 0 0, 0, , 0 . τ α µ τ µ τ τ τ + −′ ′ ′ ′ ′+′= = + < < ∑ N m m m N m mm
I f I f I (32b)
Since ( ) I τ ± comes from Eq(30) one can evaluate this expression. In the limit as τ approaches the two surfaces ournal of Computational and Theoretical Transport ( ) ( ) ( ) ( ) ( )( ) ( ) ( )
0, 0 0, 0 0, 00, 0, , α µ µτ α µ τ −′ ′ ′ ′ ′+′= +′ ′ ′′= = + = ∑∑ N m m m N m mmN m m mm
I f g f II f I (32c,d) where the required exiting intensities ( ) ( ) τ − + I I come from Eq(28b).
4. NUMERICAL IMPLEMENTATION
This section considers the additional numerical features to make the MREM a high precision method for up to seven digits of the true solution.
The Wynn-epsilon (
W-e ) convergence algorithm [12] plays a significant role in achieving high precision. The concept of convergence acceleration is to ensure the limit of a sequence of solutions, for example based on number of discrete directions N , converges to the true solution. Convergence however, does not need be sequential with increasing number of discretizations N . If we know the asymptotic behavior of the sequence, then a surrogate sequence can be established to possibly converge to the limit more quickly. The Wynn-epsilon algorithm is one such algorithm based on a non-linear recurrence to find the asymptotic trend in a sequence for a large variety of sequences, however not all. Built on directional discretization and TS, MREM offers several opportunities to apply convergence acceleration beginning with the TS. The TS representation of Eq(31) written as the limit of partial sums is ( ) lim m nnm n τ τ →∞ = = ∑ F F . (33) The partial sums therefore form a sequence to which one can apply the following non-linear
W-e recurrence by component ( )( )( ) ( ) ( ) ( )
10 , , 11 11 1 , 0,..., 2 1 ; 0,..., , ε εε ε ε ε − −+ ++ − = ≡ ≡ = + − = − = mm m i jm m m mk k k k m L f k K m L (34) where , , m i j f is the sequence of components of the matrix sequence ˆ m F of partial sums ournal of Computational and Theoretical Transport ( ) ˆ τ τ = ≡ ∑ m nm nn F F . (35) for m = 0,… . Each ( ) mk ε + represents a generally more quickly convergent approximation to the limit. These quantities form the W-e
Tableau shown in Fig. 5, with every other term on the last diagonal taken to be a more quickly convergent sequence than the original. As an example, Table I shows the advantage of convergence acceleration for the TS
Table I: TS Relative Errors for fixed h and N. m R ( Ori ) Ratio T ( Ori ) Ratio
2 9.48E-01 1.00E+00 9.48E-01 1.00E+00 3 1.39E-05 1.00E+00 1.36E-05 1.00E+00 4 2.91E-07 1.02E+00 2.81E-07 1.02E+00 5 5.03E-09 3.97E+00 4.74E-09 3.94E+00 6 7.37E-11 4.89E+00 6.64E-11 4.85E+00 7 9.49E-13 5.29E+00 7.84E-13 5.12E+00 evaluation of R and T at the slab surface for the Cloud C1 phase function considered in the next section. The table shows the relative errors for a given h used in doubling and quadrature order N and their ratio for the original ( Ori ) and
W-e sequence (not given). In our context, converged means the difference between two TS approximations is less than a desired amount. One observes the ratio (
Ratio ) is generally larger than unity indicating superiority of
W-e convergence.
Since the interaction coefficients are at the center of the calculation, their efficient determination will significantly reduce computational time. For this reason, the number of terms in the TS is limited (~15). If not converged, we then take the next h , as explained in the next section. A second example of convergence acceleration, though less dramatic, comes from the choice of h to determine the slab interaction coefficients leading to intensities m I ± from Fig. 5.
W-e
Tableau. ournal of Computational and Theoretical Transport doubling. As shown in Ref. [4], the determination of h to start doubling is nontrivial. We will exploit some known trends that enable precise determination. In particular, if h is too large, the reflectance and transmittance become unphysically negative. To capture this, we choose an initial h and check for positivity. If the check fails, we replace h by h /2 and check again until positive. Positive interaction coefficients marks the beginning of the sequence choice l h for slab thickness τ ; ,...,2 l l h l l l τ≡ = . (36a) It is generally found that after initial positivity, the coefficients remain positive for smaller h . Thus, in this way we develop a sequence of solutions in h for a fixed quadrature order. One also need be aware of round off error, which destroys the calculation if h l becomes too small. Thus, there is an optimum h . Care must be taken to treat this delicate balance. One accomplishes this by assessing both the original ( Ori ) and
W-e accelerated relative errors between consecutive h -iterations ( ) ( ) ( )( ) ( ) ( ) , 1, 1 ; ; / ;; ; / ; . τ τ ττ τ τ ± ± ± ±−± ± ± ±− ≡ −≡ − Ori l m l m l m l OriWe l m l m l m l We
Rel I h I h I hRel I h I h I h (36b) where h l is now included in the argument. The converged solution is the one satisfying , , , ± ±≤ ≤ Ori l We ll l l min Rel Rel . (36c) Since this condition may not meet the desired relative error, there is a minimum relative error achievable with the MREM found to be 10 -8 . Table II shows the corruption of relative error by round off error with increasing l for the Table II: Relative error sequence for choice of h l R ( Ori ) T ( Ori ) R ( W-e ) T ( W-e )
23 1.09E+00 1.09E+00 1.09E+00 1.09E+00 24 3.57E-09 1.33E-08 3.57E-09 1.33E-08 25 1.07E-08 4.01E-08 5.35E-09 2.00E-08 26 3.19E-08 1.20E-07 7.37E-11 7.93E-11 27 6.51E-08 2.46E-07 1.45E-08 5.57E-08 28 7.12E-08 2.66E-07 6.52E-07 2.77E-06 29 6.50E-08 2.41E-07 1.59E-06 6.08E-06 30 8.32E-08 3.11E-07 9.75E-07 3.40E-06
Cloud C1 phase function. The errors shown are averaged over 22 faux interpolation edits. The table values indicate that
W-e for l = 26 give the minimal relative error shown shaded. Note that l = 23 is the first solution with all positive interaction coefficients. ournal of Computational and Theoretical Transport Recall the solution to the radiative transfer equation is the limit ( ) ( ) lim , ; m mN
I I N τ τ µ ± →∞ ≡ ± , (37) which is precisely where
W-e would most appropriately apply. For a series of discrete N , ( ) , ; m I N τ µ± forms a sequence of solutions in N converging to the true solution as indicated. If for each τ and µ m , one introduces the sequence into the W-e algorithm [Eq(34)] the approximate solution will most likely be accelerated towards its true solution. Since convergence acceleration is an experimental mathematical procedure and experiments are not always successful,
W-e will not always better convergence that the original sequence. For the Cloud C1 phase function, convergence of the quadrature sequence based on relative error between consecutive approximations is equally effective with or without convergence acceleration. Therefore the inclusion of
W-e is more for completeness and answers the question of whether it is effective of not. Otherwise, the question would remain unanswered. Since directions change with each Gauss quadrature order, faux abscissae, which remain invariant with order, provide the measure of convergence of the angular intensity. Convergence is based on the relative error of the faux quadrature set between consecutive orders, which are given a stride of 10.
5. BENCHMARKS
We now come to benchmarks-- the focus of this work. The benchmark format follows that of Ref. [6]. We consider two standard benchmarks—the HAZE L and Cloud C1 phase functions [13], representing the state of art in radiative transfer verification. APPENDIX C gives both phase functions for completeness. With the beam source normalized to ½, the tables of intensities presented are precise to a solid six places on rounding with many to seven places. The error measure is for uniform faux quadrature abscissae in the interval [-1,1]. Both benchmarks agree to all five places of Siewert’s pioneering ADO benchmarks [6]. The first author’s response matrix method (RM/DOM) will serve as the standard since RM/DOM should be precise to one unit in the 7 th place. To the authors’ knowledge, no other benchmarks exist to the precision presented here. The first benchmark is for a HAZE L conservative scattering atmosphere of one mfp thickness. The source is a perpendicular beam on the near surface ( τ = 0). APPENDIX C gives the 82 term scattering kernel. Table IIIa shows intensities uniformly distributed in the direction interval [-1,1] correct to one unit in the last place indicated by an emboldened ournal of Computational and Theoretical Transport digit on comparison with Ref 10. Table IIIb shows six places when compared to Table IIIa. Table IIIa: HAZE L /~ Seven-places/ N ~190/ h =2 -23 /~150 s / µ \ τ τ /20 τ /10 τ /5 τ /2 /4 τ -1.000E+00 3.614515 E-02 3.4339396E-02 3.2510866E-02 2.8812216E-02 1.7628611E-02 8.525890 E-03 0.0000000E-00 -9.000E-01 3.9781870E-02 3.7872320E-02 3.5920681E-02 3.1930313E-02 1.9620173E-02 9.457313 E-03 0.0000000E-00 -8.000E-01 4.2731263E-02 4.0840607E-02 3.8873442E-02 3.4767734E-02 2.1601856E-02 1.039585 E-02 0.0000000E-00 -7.000E-01 4.800514 E-02 4.6131929E-02 4.4130697E-02 3.9829198E-02 2.5247889E-02 1.2217079E-02 0.0000000E-00 -6.000E-01 5.582135 E-02 5.404317 E-02 5.2059351E-02 4.759858 E-02 3.118365 E-02 1.5361836E-02 0.0000000E-00 -5.000E-01 6.609422 E-02 6.4629636E-02 6.2844874E-02 5.8497071E-02 4.0273976E-02 2.056212 E-02 0.0000000E-00 -4.000E-01 7.8148081E-02 7.7440255E-02 7.625076 E-02 7.270487 E-02 5.3729974E-02 2.9128534E-02 0.0000000E-00 -3.000E-01 8.9968154E-02 9.077064 E-02 9.0878384E-02 8.9471128E-02 7.2964349E-02 4.3468799E-02 0.0000000E-00 -2.000E-01 9.7081540E-02 1.0042085E-01 1.0278927E-01 1.0550594E-01 9.837771 E-02 6.7994924E-02 0.0000000E-00 -1.000E-01 9.2932814E-02 9.9818714E-02 1.0519502E-01 1.1349749E-01 1.2403692E-01 1.0839912E-01 0.0000000E-00 0.000E+00 6.9877391E-02 8.4667310E-02 9.416629 E-02 1.0872694E-01 1.3576248E-01 1.4277947E-01 0.0000000E-00 0.000E+00 0.0000000E-00 8.4667310E-02 9.416629 E-02 1.0872694E-01 1.3576248E-01 1.4277947E-01 1.1480771E-01 1.000E-01 0.0000000E-00 2.9541820E-02 5.2434564E-02 8.456491 E-02 1.3509602E-01 1.5610649E-01 1.5697621E-01 2.000E-01 0.0000000E-00 1.649068 E-02 3.2281653E-02 6.0752687E-02 1.2435036E-01 1.5892546E-01 1.7681766E-01 3.000E-01 0.0000000E-00 1.2342100E-02 2.4848764E-02 4.9396789E-02 1.1481121E-01 1.5793652E-01 1.8830112E-01 4.000E-01 0.0000000E-00 1.1187938E-02 2.264499 E-02 4.5754673E-02 1.1226864E-01 1.6086203E-01 2.000186 E-01 5.000E-01 0.0000000E-00 1.179594 E-02 2.3790964E-02 4.8000306E-02 1.1907946E-01 1.7319074E-01 2.1963289E-01 6.000E-01 0.0000000E-00 1.4204907E-02 2.8458379E-02 5.6873102E-02 1.3905102E-01 2.0144487E-01 2.5598334E-01 7.000E-01 0.0000000E-00 1.9583294E-02 3.892484 E-02 7.6745368E-02 1.8200357E-01 2.5898644E-01 3.2512495E-01 8.000E-01 0.0000000E-00 3.1953231E-02 6.2942983E-02 1.2204484E-01 2.7718191E-01 3.8276705E-01 4.686577 E-01 9.000E-01 0.0000000E-00 6.8726703E-02 1.3391677E-01 2.5425935E-01 5.4460066E-01 7.1944669E-01 8.4608373E-01 1.000E+00 0.0000000E-00 3.6493954E-01 7.0026634E-01 1.2895497E+00 2.5225517E+00 3.0931861E+00 3.3809098E+00
Table IIIb: HAZE L/Six-places/ N ~190/ h =2 -23 /~60 s / µ \ τ τ /20 τ /10 τ /5 τ /2 /4 τ -1.000E+00 3.614516E-02 3.433940E-02 3.251087E-02 2.881222E-02 1.762861E-02 8.525891E-03 0.000000E-00 -9.000E-01 3.978187E-02 3.787232E-02 3.592068E-02 3.193031E-02 1.962017E-02 9.457313E-03 0.000000E-00 -8.000E-01 4.273126E-02 4.084061E-02 3.887344E-02 3.476773E-02 2.160186E-02 1.039586E-02 0.000000E-00 -7.000E-01 4.800515E-02 4.613193E-02 4.413070E-02 3.982920E-02 2.524789E-02 1.221708E-02 0.000000E-00 -6.000E-01 5.582135E-02 5.404318E-02 5.205935E-02 4.759858E-02 3.118366E-02 1.536184E-02 0.000000E-00 -5.000E-01 6.609422E-02 6.462964E-02 6.284487E-02 5.849707E-02 4.027398E-02 2.056213E-02 0.000000E-00 -4.000E-01 7.814808E-02 7.744025E-02 7.625077E-02 7.270487E-02 5.372997E-02 2.912853E-02 0.000000E-00 -3.000E-01 8.996815E-02 9.077064E-02 9.087838E-02 8.947113E-02 7.296435E-02 4.346880E-02 0.000000E-00 -2.000E-01 9.708154E-02 1.004209E-01 1.027893E-01 1.055059E-01 9.837771E-02 6.799492E-02 0.000000E-00 -1.000E-01 9.293281E-02 9.981871E-02 1.051950E-01 1.134975E-01 1.240369E-01 1.083991E-01 0.000000E-00 0.000E+00 6.987739E-02 8.466731E-02 9.416630E-02 1.087269E-01 1.357625E-01 1.427795E-01 0.000000E-00 0.000E+00 0.000000E-00 8.466731E-02 9.416630E-02 1.087269E-01 1.357625E-01 1.427795E-01 1.148077E-01 1.000E-01 0.000000E-00 2.954182E-02 5.243456E-02 8.456491E-02 1.350960E-01 1.561065E-01 1.569762E-01 2.000E-01 0.000000E-00 1.649068E-02 3.228165E-02 6.075269E-02 1.243504E-01 1.589255E-01 1.768177E-01 3.000E-01 0.000000E-00 1.234210E-02 2.484876E-02 4.939679E-02 1.148112E-01 1.579365E-01 1.883011E-01 4.000E-01 0.000000E-00 1.118794E-02 2.264499E-02 4.575467E-02 1.122686E-01 1.608620E-01 2.000187E-01 5.000E-01 0.000000E-00 1.179594E-02 2.379096E-02 4.800031E-02 1.190795E-01 1.731907E-01 2.196329E-01 6.000E-01 0.000000E-00 1.420491E-02 2.845838E-02 5.687310E-02 1.390510E-01 2.014449E-01 2.559833E-01 7.000E-01 0.000000E-00 1.958329E-02 3.892485E-02 7.674537E-02 1.820036E-01 2.589864E-01 3.251249E-01 8.000E-01 0.000000E-00 3.195323E-02 6.294298E-02 1.220448E-01 2.771819E-01 3.827671E-01 4.686578E-01 9.000E-01 0.000000E-00 6.872670E-02 1.339168E-01 2.542594E-01 5.446007E-01 7.194467E-01 8.460837E-01 1.000E+00 0.000000E-00 3.649395E-01 7.002663E-01 1.289550E+00 2.522552E+00 3.093186E+00 For doubling, l varied form l = 21 to l = 25 and the quadrature order was about N = 190. From experimentation, it seems that MREM is more appropriate for small slabs; while, the RM/DOM is best for larger slabs. The computational time for these table is 150 s and 60 s respectively. ournal of Computational and Theoretical Transport The second benchmark is for the Cloud C1 conservative scattering atmosphere of 64 mfp thickness. Again, the source is a perpendicular beam on the near surface ( τ = 0). The 300 term scattering kernel is given in APPENDIX C and represents a true challenge. Table IVa Table IVa: Cloud C1//~ Seven-places / N ~190/ h =2 -23 /~265 s / µ \ τ τ /20 τ /10 τ /5 τ /2 /4 τ -1.000E+00 1.0636984E+00 1.0062387E+00 9.6320640E-01 8.582423 E-01 5.2453336E-01 2.4600228E-01 0.0000000E-00 -9.000E-01 9.5309008E-01 9.956622 E-01 9.697242 E-01 8.693898 E-01 5.3598880E-01 2.5740482E-01 0.0000000E-00 -8.000E-01 9.5407647E-01 9.9828274E-01 9.777659 E-01 8.805282 E-01 5.474442 E-01 2.6883574E-01 0.0000000E-00 -7.000E-01 8.8254184E-01 9.885061 E-01 9.8351863E-01 8.915683 E-01 5.5889973E-01 2.8028148E-01 0.0000000E-00 -6.000E-01 8.2471232E-01 9.7909867E-01 9.889035 E-01 9.025562 E-01 5.703551 E-01 2.9173427E-01 0.0000000E-00 -5.000E-01 7.7260568E-01 9.697724 E-01 9.939905 E-01 9.134975 E-01 5.818106 E-01 3.0319021E-01 0.0000000E-00 -4.000E-01 7.1143850E-01 9.5800446E-01 9.983295 E-01 9.2437587E-01 5.9326601E-01 3.1464747E-01 0.0000000E-00 -3.000E-01 6.4031056E-01 9.4342529E-01 1.0017757E+00 9.351789 E-01 6.0472139E-01 3.2610518E-01 0.0000000E-00 -2.000E-01 5.5848173E-01 9.2583435E-01 1.0042140E+00 9.458945 E-01 6.161767 E-01 3.3756297E-01 0.0000000E-00 -1.000E-01 4.5873404E-01 9.045925 E-01 1.0054770E+00 9.565094 E-01 6.2763204E-01 3.4902064E-01 0.0000000E-00 0.000E+00 2.5158245E-01 8.795199 E-01 1.0054752E+00 9.670143 E-01 6.3908730E-01 3.6047812E-01 0.0000000E-00 0.000E+00 0.0000000E-00 8.795199 E-01 1.0054752E+00 9.670143 E-01 6.3908730E-01 3.6047812E-01 0.0000000E-00 1.000E-01 0.0000000E-00 8.5069625E-01 1.0042130E+00 9.774071 E-01 6.5054252E-01 3.7193540E-01 7.2039069E-02 2.000E-01 0.0000000E-00 8.1871840E-01 1.0018813E+00 9.877015 E-01 6.619976 E-01 3.8339248E-01 8.8948974E-02 3.000E-01 0.0000000E-00 7.8521516E-01 9.990155 E-01 9.979409 E-01 6.7345285E-01 3.9484937E-01 1.0414040E-01 4.000E-01 0.0000000E-00 7.545985 E-01 9.967883 E-01 1.008226 E+00 6.8490802E-01 4.0630611E-01 1.1838392E-01 5.000E-01 0.0000000E-00 7.351621 E-01 9.976028 E-01 1.0187710E+00 6.9636329E-01 4.1776272E-01 1.3199043E-01 6.000E-01 0.0000000E-00 7.376549 E-01 1.0059571E+00 1.030001 E+00 7.078188 E-01 4.2921921E-01 1.4513931E-01 7.000E-01 0.0000000E-00 7.779240 E-01 1.0300617E+00 1.042781 E+00 7.1927503E-01 4.4067560E-01 1.5795774E-01 8.000E-01 0.0000000E-00 8.870455 E-01 1.0859454E+00 1.058910 E+00 7.3073267E-01 4.5213189E-01 1.7054334E-01 9.000E-01 0.0000000E-00 1.1539390E+00 1.214126 E+00 1.082650 E+00 7.421935 E-01 4.6358809E-01 1.8295432E-01 1.000E+00 0.0000000E-00 8.074596 E+01 1.1786122E+01 1.260595 E+00 7.536634 E-01 4.7504419E-01 1.9523120E-01
Table IVb: Cloud C1//Six-places / N ~190/ h =2 -23 /~185 s / µ \ τ τ /20 τ /10 τ /5 τ /2 /4 τ -1.000E+00 1.063698E+00 1.006239E+00 9.632064E-01 8.582423E-01 5.245333E-01 2.460023E-01 0.000000E-00 -9.000E-01 9.530901E-01 9.956623E-01 9.697242E-01 8.693898E-01 5.359888E-01 2.574048E-01 0.000000E-00 -8.000E-01 9.540765E-01 9.982827E-01 9.777659E-01 8.805282E-01 5.474443E-01 2.688357E-01 0.000000E-00 -7.000E-01 8.825418E-01 9.885061E-01 9.835186E-01 8.915683E-01 5.588997E-01 2.802815E-01 0.000000E-00 -6.000E-01 8.247123E-01 9.790987E-01 9.889036E-01 9.025562E-01 5.703552E-01 2.917343E-01 0.000000E-00 -5.000E-01 7.726057E-01 9.697724E-01 9.939905E-01 9.134975E-01 5.818106E-01 3.031902E-01 0.000000E-00 -4.000E-01 7.114385E-01 9.580045E-01 9.983295E-01 9.243758E-01 5.932660E-01 3.146475E-01 0.000000E-00 -3.000E-01 6.403106E-01 9.434253E-01 1.001776E+00 9.351789E-01 6.047214E-01 3.261052E-01 0.000000E-00 -2.000E-01 5.584817E-01 9.258343E-01 1.004214E+00 9.458945E-01 6.161767E-01 3.375630E-01 0.000000E-00 -1.000E-01 4.587340E-01 9.045925E-01 1.005477E+00 9.565094E-01 6.276320E-01 3.490206E-01 0.000000E-00 0.000E+00 2.515825E-01 8.795200E-01 1.005475E+00 9.670143E-01 6.390873E-01 3.604781E-01 0.000000E-00 0.000E+00 0.000000E-00 8.795200E-01 1.005475E+00 9.670143E-01 6.390873E-01 3.604781E-01 3.926386E-02 1.000E-01 0.000000E-00 8.506962E-01 1.004213E+00 9.774071E-01 6.505425E-01 3.719354E-01 7.203907E-02 2.000E-01 0.000000E-00 8.187184E-01 1.001881E+00 9.877015E-01 6.619977E-01 3.833925E-01 8.894897E-02 3.000E-01 0.000000E-00 7.852151E-01 9.990155E-01 9.979409E-01 6.734528E-01 3.948494E-01 1.041404E-01 4.000E-01 0.000000E-00 7.545986E-01 9.967883E-01 1.008227E+00 6.849080E-01 4.063061E-01 1.183839E-01 5.000E-01 0.000000E-00 7.351621E-01 9.976028E-01 1.018771E+00 6.963633E-01 4.177627E-01 1.319904E-01 6.000E-01 0.000000E-00 7.376549E-01 1.005957E+00 1.030002E+00 7.078188E-01 4.292192E-01 1.451393E-01 7.000E-01 0.000000E-00 7.779241E-01 1.030062E+00 1.042782E+00 7.192750E-01 4.406756E-01 1.579577E-01 8.000E-01 0.000000E-00 8.870455E-01 1.085945E+00 1.058911E+00 7.307327E-01 4.521319E-01 1.705433E-01 9.000E-01 0.000000E-00 1.153939E+00 1.214127E+00 1.082650E+00 7.421936E-01 4.635881E-01 1.829543E-01 1.000E+00 0.000000E-00 8.074596E+01 1.178612E+01 1.260596E+00 7.536635E-01 4.750442E-01 1.952312E-01 gives intensities uniformly in the direction interval [-1,1] to several units in the last place indicated by emboldened digits. For doubling l varied form 21 < l < 25. A six-place ournal of Computational and Theoretical Transport benchmark is given in Table IVb confirmed to all digits by Table IVa. The computational time for these tables is 265 s and 185 s respectively. CONCLUSION
A new radiative transfer numerical solution enters into the elite category of an extreme benchmark. The method, as some are, is counterintuitive in that it takes a linear problem into a non-linear one. The centerpiece of the Matrix Riccati Equation Method (MREM) is the set of four non-linear matrix Riccati ODEs that result from the fundamental radiative transfer equation as a first order ODE. Any two of these equations provides the interaction coefficients of reflectance R and T transmittance. These coefficients are numerically found from a TS representation for a thin slab to initiate doubling. Forming the interaction principle with the interaction coefficients and doubling then enables the response for slabs of any thickness. With the slab response, the surface exiting and interior intensities become known. We have established a numerical method based on TS evaluation, doubling and adding wrapped in the Wynn-epsilon convergence acceleration to give a highly precise radiative transfer solution. Tables for the HAZE L and Cloud C1 benchmarks to nearly seven-place precision are included. All exiting intensities are precise to seven digits in comparison to RM/DOM. REFERENCES [1] Efremenko, D.S., et. al.,
A review of the matrix-exponential formalism in radiative transfer , JQSRT , , pp. 17–45, 2017. [2] Grant, I.P. and Hunt, G.E., Discrete space theory of radiative transfer I. Fundamentals, Proc. Roy. Soc. A, , pp.
Diffuse reflection of light by a foggy medium , Dokl. Akad. Nauk. SSSR , pp. 229–232, 1943. [4] Wiscombe, W.J., On Initialization, Error and Flux Conservation in the Doubling Method,
JQSRT , , pp. 637–658, 1976. [5] Davis, P.J. and Rabinowitz, P., Methods of Numerical Integration , Academic. Press, 1984. [6]. Garcia, R.D.M. and Siewert, C.E.,
Benchmark results in radiative transfer , TTSP , , pp. 437–483, 1985. [7] Irvine, W.M., Diffuse reflection and transmission by cloud and dust layers , JQSRT , , pp. 471–485, 1968. [8] Chang, H.W and Wu, T.L., Numerical solutions of matrix Riccati equations for radiative transfer in a plane-parallel geometry , Waves Random Media , pp. 147–168, 1997. [9] Chang, H.W. and Wu, T.L., Matrix Riccati equation formulation for radiative transfer in a plane-parallel geometry , Waves Random Media , pp. , pp. 29–45, 1997. [10] Ganapol, B.D.,
The response matrix discrete ordinates solution to the 1D radiative transfer equation , JQSRT , , pp. 72–90, 2015. [11] Waterman P.C., Matrix-exponential description of radiative transfer . J. Opt. Soc. Am ., ,410 − Practical Extrapolations Methods , Cambridge University Press, Cambridge, 2003. ournal of Computational and Theoretical Transport [13] Lenoble, J., Standard procedures to compute atmospheric radiative transfer in a scattering atmosphere , Radiation Commission , Intl. Assoc. Meteor , Atmos; Phys ., 1977. [14] Siewert, C.E.,
A concise and accurate solution to Chandrasekhar's basic problem in radiative transfer , JQSRT , pp. 109–130, 2000. [15] Ganapol, B.D., An Extreme Benchmark for Monoenergetic Scattering in Hydrogen , ANS Winter Meeting, Orlando FL., 2018. [16] van de Hulst, H.C.,
A New Look at Multiple Scattering , NASA Institute for Space Studies, New York, 1965. [17] Ganapol, B.D. and Patel,J.,
DPN Solution of the 1D Monoenergetic Neutron Transport Equation with Benchmarking , ANS Winter Meeting Wash. DC, 2017.
APPENDIX A Previous Benchmark Solutions to the 1D Radiative Transfer Equation
A.1 The response matrix discrete ordinates method (RM/DOM) [10]
The response matrix solution is similar to the Analytical Discrete Ordinates (ADO) method [14], with one major difference. While ADO generates a solution based on the discrete form of the singular eigenfunctions leading to exponential eigenfunctions, RM/DOM assumes the even/odd parity form of the transport equation leading to modal solutions of hyperbolic type in the space of solutions of second order ODEs. Matrix diagonalization then yields two independent solutions. One advantage, apart from compactness, is that there always exists two independent solutions to the parity equations even in the case of no absorption. Thus, special attention is avoided for ω = 1. In addition, the solution features the development of a response matrix (as the name suggests) over a homogeneous slab of any thickness dependent on material properties only. Responses, seamlessly added, represent contiguous heterogeneous media using the interaction principle as described above. Additional features include convergence acceleration in quadrature order through Wynn-epsilon acceleration; faux quadrature to assess precision with quadrature grid refinement also enabling a beam source to be specified in any direction; and a consistent solution at µ = 0. We have achieved seven-place precision for the Cloud C1 300 term scattering kernel with quadrature order 2 N of ~300. A.2 Method of doubling (MoD) [16]
This solution closely follows that developed by Van De Hulst [16] but applied to the fully discretized 1D radiative transfer equation in both space and direction. The slab is discretized and a response for a thin slice of the discrete slab is continuously doubled until the entire discrete slab is covered. The full slab is then covered by adding via the interaction principle. Richardson and Wynn-epsilon convergence acceleration applies to the initial thin slice size and to the quadrature order respectively. In addition, Richardson extrapolation is applied to the original slab discretization. The features found in RM/DOM are also included. Heterogeneous media are naturally accommodated with this approach. The Cloud C1 case is as precise as that of RM/DOM. The advantage of doubling is that ournal of Computational and Theoretical Transport no matrix diagonalization requiring eigenvalues is necessary and essentially a host of spatial discretization schemes to evaluate the matrix exponential as a solution to a first order ODE are possible. These include single and multisteps methods. The doubling solution demonstrates that the most basic of numerical schemes provide highly precise numerical solutions to the radiative transfer equation thus demonstrating simplicity. A.3 Double PN method (DPN) [17]
The DPN method, though theoretically equivalent to the discrete ordinates method, is not necessarily numerically equivalent. The usual half range Legendre polynomial expansion is assumed yielding moments equations closed by setting the last moment to zero. These equations form a set of ODEs solved in the even/odd parity form. Again, as in the RM/DOM, matrix diagonalization enables an explicit solution yielding a response matrix. With all the features included in the previous solutions, DPN also gives nearly seven places for the Cloud C1 case. What is particularly interesting about RM/DOM and DPN is that they both use the same solutions as found in solving the diffusion equation, which gives another connection between diffusion and transport theory. This is not surprising since they are all solving second order ODEs.
APPENDIX B
Indirect Taylor Series Representation B.1 Reflectance
One can linearize Eq(12c) with the following transformation: ( ) ( ) ( ) τ τ τ − ≡ R Q P (B1) as shown. Since from Eq(B1) in Eq(12c), by differentiation and grouping, one finds [ ]
11 1 1 1 1 + . d dd d τ τ −− − − − − + Q PP Q = P Q P QP QP QP β −α − α β (B2) If one defines dd τ ≡ Q P Q β −α , (B3) and pre-multiplies by P , then + dd τ − − = PP QP −α β . (B4) ournal of Computational and Theoretical Transport Noting d d d dd d d d τ τ τ τ − − − = = = + PP I P PP P , Eq(B4) becomes + dd τ − − − = P P QP −α β to give dd τ = P P Q α −β . (B5) If ≡
PY Q , (B6) Eqs (B3) and (B5) combine as dd τ ≡ Y AY (B7a) with initial condition from ( ) ( ) ( ) − ≡ = R Q P satisfied by ( )( ) ( ) N − ≡ ≡ QP = P I to imply ( ) N ≡ IY . (B7b) ournal of Computational and Theoretical Transport To find the components of ( ) τ Y , we impose the TS ( ) nnn τ τ ∞= = ∑ Y Y (B8a) to give ( ) ( ) n n n n n τ τ −= ≡ = Y Y AY . (B8b) and therefore from Eq(B7a) ( ) ( ) [ ] ( ) ( ) [ ] n n n nn n n n n nn n ττ ττ − −= − −= ≡ = −≡ = − P P P QQ Q P Q α ββ α (B9) For the TS coefficients for R of Eq(B1), we write ( ) ( ) ( ) τ τ τ≡ P R Q , (B10) and substituting TS nn n l l nl −= + ∑ P R P R = Q ; (B11) and therefore the indirect TS coefficients for the reflection coefficient nn n n l ll −= − ∑ R = Q P R . (B12)
B.2 Transmittance
For transmission, we begin with Eq() dd τ T = R T T β −α (B13) into which we introduce ournal of Computational and Theoretical Transport − T = WP , (B14) to find
111 1 1 . d d dd d d τ τ τ −−− − − + T W P= P W= QP WP WP β −α (B15) Since
11 1 1 1 d dd d τ τ −− − − − + W PP W= QP WP WP β −α (B16a) and from Eq(B4) above + dd τ − − − − = P P P QP − α β , (B16b) there results + . dd τ − − − −− − − = W P WP WP QP= QP WP WP − α ββ −α (B16c) From the second identity of Eq(12f) + − − − − − − = QP WP WP WP WP QP β −α − α β resulting in dd τ = W . (B17a) Thus, ( ) τ W = C (B17b) and since ournal of Computational and Theoretical Transport ( ) ( ) ( ) ( ) N N − = T = I = W P W I and ( ) ( ) N τ = W = W I . (B17c) The TS emerges when Eq(B14) is re-expressed as ( ) ( ) N τ τ P T = I (B18) to give the TS coefficients n n l l N nl δ −= ∑ P T = I and therefore
10 0 nn N n n l ll δ − −= = − ∑ T I P T . (B19) and − T = P . (B20) One finds the identical result if Eq(12b) were used as the initial equation rather than Eq(12d). Seems, the last equation is new and implies ( ) ( ) ( ) τ τ τ
R = Q T . (B21) Additional implications of this expression needs further study in particular for the eigenvalue based radiative transfer formulation [10].
APPENDIX C
Legendre Coefficients for Scattering Phase Functions
The following are the scattering coefficients for HAZE L and Cloud C1 phase functions represented to orders L =82 and 300 respectively as the Legendre series ournal of Computational and Theoretical Transport ( ) ( ) ( )
1, '2
L ll f P P µ µ β µ µ = ′ ≡ ∑ . (C.1) C.1 HAZE L