DASA-1257 The University of Michigan Department of Chemical and Metallurgical Engineering Thermal Research Laboratory EXACT SOLUTIONS FOR ANISOTROPIC, MULTIPLE SCATTERING BY PARALLEL PLANE DISPERSIONS by, ~ _ Stuart W. Churchill Professor of Chemikcal Engineering. Chiao-Min Chu":il".' Associate Professor of Electrical Engineering Larry B. Evans Li-Chiu Tien Sing-Chin Pang Assistants in Research This research was supported by the Defense Atomic Support Agency Contract No. DA-49-146-XZ-039 Annual Report September 1961 "Request for copies of this report should be submitted to ASTIA, Arlington Hall Station, Arlington 12, Virginia."

t7 LC

TABLE OF C0ITEIEIS LIST OF FIGURES vii LIST OF TABLES xi NOMENCLATURE xv ABSTRACT xix INTRODUCTION PART I —NUMERICAL SOLUTION OF THE TRANSPORT EQUATION 4 A. Mathematical Formulation 4 1. Geometry and Coordinate System 4 2. / and 4 Functions 7 3. Intensities of Reflected and Transmitted Radiation 10 4. Integrated Reflection and Transmission 12 B. Numerical Procedure 15 1. The General Numerical Problem 15 2. Specific Numerical Techniques Used in Computation 18 C. Solution to Specific Problems 21 1. Choice of Independent Variables 21 a) Choice of Phase Functions 23 b) Choice of Other Independent Variables 29 2. The q and q Functions 31 3. Sample Values of the Intensity of Reflected and 32 Transmitted Radiation 4. The Integrated Reflection and Transmission 42 60 5. Discussion of Numerical Techniques iii

PART II —EVALUATION OF APPROXIMATE MODELS 65 A. Approximate Models 65 1. Two-flux Model 65 2. Six-flux Model 66 3. Richard's Modified Diffusion Theory 67 4. Numerical Results 68 B. Interpolation Factors 76 PART III —DEVELOPMENT OF IMPROVED APPROXIMATE MODELS 89 A. Introduction 89 B. Approximate Differential Equation 90 C. Infinite Medium With Point Source 96 D. Boundary Conditions 101 E. Solution of One-Dimension Problem 107 1. Second Approximation 107 a) Medium Bounded by Parallel Plate b) Semi-infinitive Medium 110 2. First Approximation 110 F. Calculation of Albedo 113 1. First Approximation 113 2. Second Approximation 114 G. Numerical Results 116 H. Conclusion 124 125 SUMMARY AND CONCLUSIONS iv

APPENDIX A - Derivation and List of Equations 128 1. Derivation of Explicit Relations for bm(lr,llY4o) and 128'm ( 1n'' 1o) 2. Integrated Reflectance of Half Space for Rayleigh 131 Scattering 3. Solution for the Six-flux Formulation of the Slab 135 Problem APPENDIX B - Tables of Computed Functions 138 APPENDIX C - Computer Programs 160 1. Solution of the Integral Equations 160 2. Subroutine to Calculate Associated Legendre 172 Polynomials 3. Program to Calculate Integrated Reflectance and 174 Transmission REFERENCES 178 DISTRIBUTION LIST 180

LIST OF FIGURES I-1 Geometry of the problem 5 I-2 Coordinate system 6 I-3 Characteristics of a two term phase function 27 I-4 Phase functions used in calculations, F = 0.50, F = 0.75, 28 and F = 0.933, Parameters of P I-5 Intensity of reflected and transmitted radiation as 35 function of bearing angle, I/Io vs. f, parameters of Co' L = 0.50, )o= 0.90, U1 = 1.00, F = 0.9330, and P = 0.4821 I-6 Intensity of reflected and transmitted radiation as 36 function of bearing angle, I/Io vs. C, parameters of C',,o = 0.50, 0Wo= 0.90, T1 = 1.00, F = 0.9330, and P = 0.4821 I-7 Intensity of reflected and transmitted radiation as 37 function of bearing angle, I/Io vs. A, parameters of SO 0o = 0.50, L = 0.76923,'1 = 0.50, F = 0.9330, P = 0.4821 I-8 Intensity of reflected and transmitted radiation as 38 function of bearing angle, I/Io vs.', parameters of 1,' 1O = 0.50, I = 0.76923, Wo.= 0.90, F = 0.9330, and P = 0.4821 I-9 Intensity of reflected and transmitted radiation as function 40 of bearing angle, I/Io vs T, isotropic and Rayleigh scattering, i = 0.04691 and 0.76923, o0 = 0.50, QO= 0.90~, 1 = 1.00 I-10 Intensity of reflected and transmitted radiation as function 41 of bearing angle, I/Io vs. T, parameters of P, Lo = 0.50, = 0.76923, Qo= 0.90, 1l = 1.00, F = 0.750 I-11 Integrated reflectance as function of angle of incidence for 43 isotropic and Rayleigh scattering, R vs. So, parameters of T1, oo = 0.90 I-12 Integrated reflectance as function of angle of incidence 44 for isotropic scattering, R vs. o, parameters of 1, 0= 0.60 and 0.30 vii

I-13 Integrated reflectance as function of angle of incidence 45 for anisotropic scattering, R vs. po, parameters of T1l o0o= 0.90, F = 0.9330, P = 0.4821. I-14 Integrated reflectance as function of angle of incidence 46 for anisotropic scattering, R vs. 4, parameters of T1i 0o = 0.60 and 0.30, F = 0.9330, P = 0.4821. I-15 Total integrated transmission as function of angle of 47 incidence for isotropic and Rayleigh scattering, T vs. io, parameters of'1, c,= 0.90 I-16 Total integrated transmission as function of angle of 48 incidence for isotropic scattering, T vs. o0, Parameters of cWo, ~1 = 0.05 and 0.50 I-17 Total integrated transmission as function of angle of 49 incidence for anisotropic scattering, T vs..40, parameters of T1i C0o= 0.90, F = 0.9330, P = 0.4821 I-18 Total integrated transmission as function of angle of 50 incidence for anisotropic scattering, T vs.,u, parameters of ~oy,'1 = 0.05 and 0.50, F = 0.9330, P = 0.4821 I-19 Integrated reflectance as function of angle of incidence, 51 R vs. 40, parameters of F, c0W= 0.90 and,1 = 1.00 I-20 Total integrated transmission as function of angle of 52 incidence, T vs. o0, parameters of F, 0= 0.90 and T1 = 1.00 I-21 Integrated reflectance and total integrated transmission as function of optical thickness of dispersion, R and T vs. T1' parameters of 40 isotropic or Rayleigh scattering F = 0.50 I-22 Integrated reflectance and total integrated transmission 54 as function of optical thickness df dispersion, R and T vs.' 1' parameters of o, anisotropic scattering F = 0.9330 and P = 0.4821 I-23 Integrated reflectance as function of albedo for single 55 scattering, R vs. c30, parameters of F, ~o, and F1. viii

I-24 Total integrated transmission as function of albedo for r6 single scattering, T vs. o0, parameters of F, 40,.l' I-25 Integrated reflectance and total integrated transmission 57 as function of phase function, R and T vs. F, parameters of ~o U3o- 0.90,'1 = 1.00 II-1 Integrated reflectance for various solutions as function 69 of angle of incidence for isotropic scattering, R vs. 40o 0O0 = 0.90, T1 = 1.0 II-2 Integrated reflectance for various solutions as function 70 of angle of incidence for isotropic scattering, R vs. 40o parameter of T1, 3o= 0.90 II-3 Integrated reflectance for various solutions as function 71 of angle of incidence for anisotropic scattering, R. vs.,o' parameter of T1, c0.= 0.90, F = 0.9330, P = 0.4821 II-4 Total integrated transmission, various solutions as function 72 of angle of incidence for anisotropic scattering, R vs. [o, parameter of 0o, T1 = 0.50, F = 0.9330, P = 0.4821 II-5 Total integrated transmission, various solutions as function 73 of angle of incidence for anisotropic scattering, R vs. [o, parameter of P, Lo3 = 0.90,'1 = 1.0, F = 0.750 II-6 Integrated reflectance, exact, interpolated and extrapolated 77 as function of angle of incidence for isotropic scattering, R vs. 4OL parameter of'r1, ~o = 0.90 II-7 Integrated reflectance, exact, interpolated and extrapolated 78 as function of angle of incidence for isotropic scattering, R vs. o0, parameters of 3o0 and t1 II-8 Factor x in equation II-14 calculated from isotropic 81 Ro and R1, as function of angle of incidence, x vs. 40o parameter of Lu3 II-9 Factor x in equation II-14 calculated from isotropic 82 Ro and R1 as function of angle of incidence, x vs. 4o, parameters of T1

II-10 Integrated reflectance, exact and interpolated as 84 function of angle of incidence for anisotropic scattering, R vs. o,, parameter of T1, )o3= 0.90, F = 0.9330, P= 0.4821 II-11 Integrated reflectance, exact and interpolated as 85 function of angle of incidence for anisotropic scattering, R vs. iO, parameters of e& t1 III-1 4Er(r) as a function of r for an infinite medium 100 with point source III-2 Integrated reflectance as a function of incident angle, 117 R vs ~o, parameters of 300, and T1 = oo III-3 Integrated reflectance as a function of albedo for single 118 scattering, R vs h3, o = 0, and T1 = oo III-4 Integrated reflectance as a function of albedo for single 119 scattering, R vs. M3O, o = 0.5, and T1 = oo III-5 Integrated reflectance as a function of albedo for single 120 scattering, R vs.C oLo = 1.0, and T1 = oo III-6 Integrated reflectance as a function of incident angle, 121 R vs. ko0 parameters of ZO, and T1 = oo III-7 Integrated diffusive transmission as a function of 122 incident angle, TD vs. 4o, parameters of o., and T1 = 1 III-8 Total integrated transmission as a function of incident 123 angle, T vs uo, parameters of co, and T1 = 1

LIST OF TABLES I-1 Phase functions used in computation 26 I-2 Values of the parameters used in computation 30 I-3 Values of the coefficients used in equations I-21 and I-22 34 to calculate the specific intensities of reflected and transmitted radiation I-4 Comparison of'\(f9) functions for isotropic scattering 61 and 3~ = 0.90 with Chandrasekhar's H functions I-5 Comparison of the w(/.) and <(e*) functions for isotropic 61 scattering and,o = 0.90 with Chandrasekhar's X and Y functions I-6 Comparison of the ~ and $ functions of this report 62 with similar functions of reference (5) calculated using a 20-point Simpson's rule technique for al = 2.73458, a2 = 2.24153, Q0,= 0.90 and TL1 = oo I-7 Comparison of integrated reflectance for Rayleigh scattering, 62 T1 = oo and m<,= 0.90 with values calculated using the method described in Appendix A-2 II-1 Interpolation and Extrapolation of R for isotropic 79 scattering with respect to 4o and 1 II-2 Interpolation of R with respect to 40 for anisotropic 83 scattering II-3 Interpolation of T with respect to 4o for anisotropic 87 scattering III-1 The explicit values of l1, a2, ml, and m2 as a function 98 of so xi

A-1 The K-functions and related numerical values for 134 computing the integrated reflectance of half space, Rayleigh scattering B-0 The V and 4: functions, integrated reflectance, diffuse 139 portion of the integrated transmission, and total integrated transmission for 36 sets of parameters (Listed in Table I-2) Problem 1 139 Problem 2 139 Problem 3 139 Problem 4 140 Problem 5 140 Problem 6 140 Problem 7 141 Problem 8 141 Problem 9 141 Problem 10 142 Problem 11 142 Problem 12 142 Problem 13 143 Problem 14 143 Problem 15 144 Problem 16 144 Problem 17 145 Problem 18 145 Problem 19 146 Problem 20 146 Problem 21 147 Problem 22 147 Problem 23 148 Problem 24 148 Problem 25 149 Problem 26 149 Problem 27 150 Problem 28 150 Problem 29 151 Problem 30 151 Problem 31 152 Problem 32 152 Problem 33 153 Problem 34 153 Problem 35 154 Problem 36 154 xi1 Problem~~~~~~~~~~ 34

B-1 Comparison of exact and approximate integrated reflectances, 155 isotropic and Rayleigh scattering B-2 Comparison of exact and approximate total transmissions, 156 isotropic and Rayleigh scattering B-3 Comparison of exact and approximate integrated reflectances, 157 anisotropic scattering, F = 0.933, P = 0.4821 B-4 Comparison of exact and approximate total transmissions, i18 anisotropic scattering, F = 0.933, P = 0.4821 B-5 Comparison of exact and approximate integrated reflectances 159 and total transmissions, anisotropic scattering, F = 0.750, co3 = 0.90,'1 = 1.0 xiii

NOMENCLATURE ai angular distribution coefficients, i = O 1,...., N m i Ai coefficients in expansion of cos m = A0 + Z Ai cos A=j bi(T'ly4y o) coefficients in expansion 1-21 of angular dependence of specific intensity of reflected radiation. Defined in Appendix A-1 B fraction of radiation scattered backward. Defined by I-47 B6 six-flux representation of backward scattering component BP height of backward peak of phase function, p(-1.0) coefficients in expansion I-22 of angular dependence of diffuse portion of transmitted radiation. Defined in Appendix A-1 D diameter of scattering particle En(Z) defined by A dX f(g) angular distribution function, same as p(cos Q)/4L f(x) any function of x f*(x) any tabular function of x f*(n)(x) approximation to solution of I-39 after n iterations F fraction of radiation scattered forward, defined by I-46 F(Q) any function of i F6 six-flux representation of forward scattering component Fk(Ti,,' L) function defined by I-3 FP height of forward peak of phase function, p(1.0) function defined by I-4 H(~) Chandrasekhar's H function, the same as I, (oo, u) XCV

Tl'I 2'" discrete components of integrated radiant flux I(Z,p, t ) specific intensity of radiation at position Z in the solid angle defined by 1i and 0 0I0(-[ Io ) intensity of a parallel beam incident on plane parallel dispersion in the direction (-Io' y ) m index of refraction of scattering particle mi constants derived in equation (III-29, 30) from the approximation of the isotropic scattering integral equation n number of particles per unit volume ni coefficients of the approximate kernal derived in equation (III-18) N number of terms in the phase function p number of ordinates (excluding the zeroth ordinate) used in numerical integration with I=35 p(p) phase function P peakedness of the phase function in the forward direction defined by I-49 PB peakedness of the phase function in the backward direction defined by I-50 Pl(C,) Legendre polynomial P1M(,) associated Legendre polynomials ri i-th root of Pn(x) R integrated reflectance S(r) Virtual source function -S( (-r,; Vf) scattering function defined by I-11 c(m);',,o) function defined by I-13 xvi

t thickness of dispersion T total integrated transmission TD diffuse portion of integrated transmission 7 (ro;?,;Po,o) transmission defined by I-16 7 (T,;/,() function defined by I-18 W. weighting coefficients for use in numerical integration by I-33, i= O...., P x. specific values of ordinates to be used in numerical integration by I-33, i = O...., p x Ao-interpolation factor for integrated reflectance y 0o-interpolation factor for total integrated transmission X six-flux representation of sidewise scattering component X(A) Chandrasekhar's X function —the same as ~( Y(~) Chandrasekhar's Y function —the same as +0(t,,/) Z distance into the dispersion in optical units ta WD/A, circumference of particle expressed in wavelength positive roots of equation (III-26) Z i ~i positive roots of the equation 4-. + cO, 5 b=5 4 where a = 2, b=-, c 6SA Kronecker delta ( <,= 0 when m i n;,.= 1.0 when m = n) error in numerical integration with I-35 0 azimuthal angle of radiation 90 azimuthal angle of incident radiation xvii

A wavelength of incident radiation cos 0 direction of radiation cos 0 direction of incident radiation p(o energy density of isotropic scattering at a position r 1i optical thickness of dispersion ~' Total cross-section of a particle (scattering and absorption) functions defined by I-2 bearing angle of radiation Ifo bearing angle of incident radiation functions defined by I-1 CWo albedo for single scattering, ratio of scattered to intercepted radiation coefficients in Legendre polynomial expansion of phase function U(- m)! t4, It + 13) Ixviii xviii1

ABSTRACT The effect of anisotropic scattering on radiative transfer was investigated theoretically. Exact solutions were obtained for a plane source obliquely incident on a multiply scattering, parallel plane dispersion by numerical integration. The results are utilized to evaluate the various approximate models which have been proposed for radiative transfer and neutron diffusion. A new variable-order, diffusion-type model is proposed as an improved approximation. xix

INTRODUCTION The long range objective of this research is to develop methods for predicting the transport of thermal radiation through the atmosphere. The scattering of radiation by dispersed material presents a problem of considerable mathematical difficulty. Exact solutions have been completed only for highly idealized conditions: isotropic or Rayleigh scattering in a one-dimensional configuration (1). Accordingly many approximate models have been developed (2). Examples are the diffusion model and the various discrete flux models. In the initial phase of this research a discrete flux model was used to obtain approximate results for a finite, spherical source both above and inside a haze (3). The work was then redirected to the development of exact solutions for particular conditions as discussed below. The lack of exact solutions for general conditions makes it difficult to compare and evaluate the various approximate models which have been proposed. As a consequence of a general conference on radiant transport (4) and subsequent discussions with Dr. S. Chandrasekhar and AFSWP personnel in which this situation was stressed, the immediate objective of this research was shifted to the development of exact solutions for idealized but significant conditions. Particular attention

has been given to the development of exact solutions for anisotropic scattering. Solutions were first obtained for parallel-plane radiation, obliquely incident on a semi-infinite dispersion. These results are presented in the annual report for the previous year (5). The present report presents similar results for parallel-plane radiation obliquely incident on parallel-plane dispersions of finite thickness. The method of solution follows the lead of Chandrasekhar (1) and involves the reiterative numerical solution of integral equations on a high speed computer. (The University of Michigan IBM-704). The formulation of the method of solution and the numerical results are presented in Part I of the report. In Part II numerical values obtained from various approximate models are compared with the exact values presented in Part I. Suggested formulas for interpolation and extrapolation of the exact values are also included. None of the approximate models which had previously been proposed proved to be entirely satisfactory and a new variable-order, diffusion-type model is developed in Part III. During the period while this work was being done an approximate but exceedingly detailed incremental model for one-dimensional radiant transport with anisotropic scattering was developed by the Internuclear Corporation (6). Their final results have not been available for comparison with the values obtained in this investigation but comparison of preliminary values for a few cases indicates that their model may be quite accurate.

In the continuation of this research it is planned 1) to refine the exact method of solution to expedite the numerical calculations 2) to extend the new variable order diffusion-type model for anistropic scattering and 3) to develop exact or nearly exact solutions for two-dimensional problems utilizing numerical methods and/or the model.

PART I - NUMERICAL SOLUTION OF THE TRANSPORT EQUATION A. Mathematical Formulation 1. Geometry and Coordinate System. This investigation is concerned with the irradiation of a semi-infinite slab by a uniform parallel flux as indicated in Figure I-1. The Z = 0 plane represents the interface between region (1) consisting of free space containing the source and region (2) consisting of a uniform dispersion of multiply scattering and absorbing particles each characterized by the same absorption cross-section, scattering cross-section and phase function, p(cos @), for single scattering. The slab has an optical thickness T1 and the Z ='1 plane represents the interface between region (2) and region (3) consisting of free space. The objective of this study is to determine the intensities of transmitted and reflected radiation as a function of position, direction, and the characteristics of the source and dispersion. The coordinate system used in the analysis is shown in Figure I-2. An incident flux of intensity Io strikes the interface at an angle GO with respect to the normal (Z axis). The direction at Z = 0 of the reflected intensity leaving the dispersion is defined by the azimuthal angle @, measured from the normal, and the bearing angle j measured from the Y axis. For simplicity the coordinate system is

Region (3) II1~'i Rgion (2)T, Region (I) Figure I-1 Geometry of the Problem

,I,, X I l0 XY plane: media interphase Io: incident flux I: emergent flux under consideration eo =cos-I oe= cos'- I a.,I= Bearing angle Figure I-2 Coordinate eystem 6

oriented so that Po = O, i.e., I0 lies in the YZ plane. 2. The 4 and 0 Functions. The following integral equations for the two basic functions 1 and 01 have been derived by Chandrasekhar(l)*. ytm /i) m CJ) + ( - 1 ) @; L'P~ (z,u)j(ri' I -1 _ m (1" Z,,,'Y, t 2, ~l(/ ) L4 a m -0m, in.., where, p. = cos Q = optical thickness of media (number of mean free paths for scattering) = naot P",(,) = associated Legendre polynomials N = number of terms in the expansion of the phase function N p(p.) X: 1P (p.) 1=0 b01 = coefficients of phase function expansion m (1-m)! "'1::1 (l+m)! * Hereafter referred to as R. T.

The equations may be written in more convenient form by defining m C ~ ~I ~I m m s (tt, _, tzl,, )('I (CTr' (efA,,M))L I-4 The equations then become mp ) P,,,r') P ( )J I-5 IC 0, N, ~(,,tf)C= e ~ (A) + 2 t/ (,,?,' )P'(/')d4' I-6 m= O, i,', N =, mi, +,, N For any value of N, equations I-1 and I-2 yield N + 1 systems of simultaneous integral equations. Each system consists of 2(N+1), 2N,..., 2 simultaneous equations corresponding respectively to m = N, N-l,...O. As an example, using the notation F1 = and G = Gi(T1,,p1I'), the systems of integral equations for N = 2 are for m = 0 2?,L) = 1 F - FI dj4' I-7a'V^ ( - F, -F r F2 J'4F' I-7b t, (,,s)= 5j1 (3rk- 1) + [ FI - F - + F) LI I-7c +02~o =e2 1,[ (.r,,p/)= er//LI c: * "I-7 It''8

#:(z.,/),e te-t* * t/A; ( + e+ 1 1, -7e + (z,1 ) = ( (3/2 _ -)e -UC' +4 * C & 1.SC + C 0 + I-7f for m = 1 4,< (r),?)- F. (? _ /2 d,' I-8a, (?/) - 3 (O F )- 2t4Q tF-Fai 3 ( tA1 -8b ('1,~ 1 # c~4'' ],#')d' I-8c m =2 q'(r,)=3CI -/t2) + fu F3(G)3(1 T-9a 0 t'?~

3. The Intensity of Reflected and Transmitted Radiation The intensity of radiation reflected diffusely by a semiinfinite slab at the plane Z = 0 in the direction ( -, ~) may be expressed in terms of a scattering function* 1-10 I (o,?, p ) = lo (-?o04errj) where I(-.,o) = incident flux intensity in the direction (-,oO). I(o,,, ) = intensity of flux reflected in the direction (I,?) at Z = 0. S = Scattering function The scattering function is** N (m) (Z.; /,Q; \.) =y S. (V, b) f CoSm (~. - ) or for 0o = 0 m=Q S ws Pn I-12 where, si(c~.~?,?.) tf2 ~.,.~?/,. Z (-0"' S.(,,A 2 S,, )e in(' F (,,AHe I-13 * R. T. Page 161 ** R. T. Pages 180 and 177 10

The diffuse reflected intensity may then be written as Iwi_(~/o,+7, Cos rf 2 - So) (l) - I-14 sN m=Om In a similar manner the intensity of radiation transmitted diffusely through the slab at the Z = T1 plane may be expressed in terms of a transmission function* Z(,, )= ( t,+ Cr,; An Y; /Alt (fo) I-15 where, i.(-jOJ'~')= incident flux intensity in the direction (-Tok'yo) I(z,,-,yc)= intensity of flux at Z = 1 transmitted in the direction (-, j) = transmission function The transmission function is** J(,~;?, ~;? o,,))= E 1t (,, 4o) P.) -os0~6 M( (,-+) I-16 O;Im s C a /. 0 -7 w (e,;;",~; t,'i'.)= 7 (T,?, o)c mo I-17 where, r7m(r~:.,.4,?.) j,,,. (2z- Ss) ~ ~, ( z~,,A, o)i-18 The diffuse portion of the transmitted intensity may then be written as I*(RIo) T. age (2 - 6o ) co 1i (r v I-19 * R. T. Page 161 ** R. T. Pages 180 and 177 11

Equations I-14 and I-19 may be written in a slightly more convenient form by noting that m cos5mT -- A, + ~ Ai cos'f I-20 j=:1 where, for example, for m equal to 2, AO = -1, A1 = 0, and A2 = 2. Thus, new expressions can be derived for the reflected and the diffuse portion of the transmitted intensities of the following form N rI ( )o,, = bo (s,, pJo,) + Z bm (i,,,a o ) CoS"' I-21 and I (t:,,A,,) c. (Ir,, 4,?.) E C,(Zl,?,?o) Co~S 1 -22 This means that for a N term phase function and a given combination of!1,.,, and io, the intensities need be calculated at only N + 1 values of 4 in order to obtain the complete angular distributions. The explicit relations for bm(l,p,p.o) and m M Cm(Tl,~,4o ) in terms of *1(l,.,op) and 1(~l,4,.o0) for N = 2 are derived in Appendix A-1. 4. Integrated Reflectance and Transmission The integrated reflectance is defined as the ratio of reflected to incident power passing through a unit area of the interface. The flux, Io0, is defined as the radiant energy passing through a unit area normal to the direction of propagation per unit time. Therefore, the 12

incident power striking a unit area of the interface is Io x [o0 Similarly, the reflected power passing out of a unit area of the interface is I( Ie (O PIf)rud. The integrated reflectance is then expressed as R = ag i I 0( o t I f )p JUO fI-23a Replacing I(O,t,f) by I-14 and noting that 2Tr I F(/) cosStnJ 2 rTF() or o for n o or m~yo respect;velyI-24 F? = 2 2 t l~l 5 p Ft (T1,,,8,,h,) 1,"I-23b lgo In a similar manner the integrated diffuse transmission Tr is defined as the ratio of transmitted to incident power passing through a unit area of interface. Thus Ts = I <Ad Ads-5 9! I,, -/,f)d I-25a and substituting I-19 for I(1r,-1,1) TV'-T rip 2; D @ (Z.,, e)dpI-25b L=o The total transmission T is the sum of the direct transmission and the diffuse transmission. Thus T e + T I-26 It has been shown(5) that, since, N i 13

then x1 I) FO C(',, Pu, I -28 and R = i - for N 3 1,;o 0 I -29?,, Similarly it can be shown that, since ano Then g (r,/ LA). 1-31 1A ) e'/=2; j rT I-31 and I-32 For N i, S, >o The total transmission is simply &I, ( ~,,?,) T- - I-33 14

B. Numerical Procedure 1. The General Numerical Problem. As shown in the previous section, the intensities of reflected and transmitted radiation for a particular phase function for single scattering can be expressed in terms of the solutions of sets of simultaneous, non-linear integral equations. The numerical problem, then, is to solve these sets in order to obtain values for the intensity, I, in terms of its parameters. The most straightforward method of obtaining numerical solutions of integral equations is simple iteration using the original equations. For example, in the case of a single integral equation of the type f(x) = g(x) + A 1 f(x)f(xt)h(x')dx' 1-34 where f(x) is the function to be found g(x) and h(x) are known functions of x and A is some known constant, the method involves the following: 1. Replacing the continuous function, f(x), by a tabular function, f*(x) where f*(x) consists of a set of p values of f(xi), {fo,fl,..fpt, corresponding to a set of values of Xi, XoXl,...Xp I where p is the number of points. 2. Choosing an initial f*(O)(x). 3. Evaluating the right side of I-34 by some mechanical quadrature method —thus, obtaining a new f* (1)(x).

4. Testing to see if f*(l)(x) agrees with f*(0)(x) everywhere within acceptable limits and, if not, 5. Repeating the process using f*(i)(x) in place of f*(O)(x). 6. The process is continued until the sequence: f*(O)(x), f*(l)(x),... f*(n)(x) converges to some limiting value. Advancing one step in the iteration procedure to obtain f*(n+l)(x) from f*(n)(x) requires making p+l numerical integrations. Integration by mechanical quadrature consists of replacing the integral by a sum, or in other words (t f(x)dX = 4 wi 0(X;) + I-35 where ( is the error. Thus, selecting a specific method of numerical integration involves selecting the set of values of x,IXo,x, XP, and the set of weighting coefficients, {W,Wo...W}. For integration by the trapezoidal rule Xi = i/p I-36 W= Wp = 1/2p Wi = 1/p (i=l,...,p-l) For integration by Simpson's rule p = 2k where k is an anteger Xi = i/p 1-37 W = Wp = 1/3p Wi = 4/3p i = 1, 3,...,p-l Wi = 2/3p i = 2, 4,...., p-2 16

For integration using n-point Gaussian quadrature X; = (fit + 1)/2 I-38 where ri is the i-th root of Pn(x) I PA ( "(x)lx W; 0, (K) PA(X)dx ZP~(x) IX: —Thus, the definition of a specific numerical procedure involves: 1. Selecting the number of points, p+l, ixii and {Wi. 2. Selecting the initial function f*(O)(x). 3. Selecting a test to determine when the sequence f*(O)(x), f*(l)(x),.., f*(n)(x) has converged. It should be noted that the fact that the function f*(x) has converged to some limiting value merely assures that a solution has been found to the following equation: 4f(xi)' - h f(i) fx f x;) x;) -39 j=0 for j = 0, 1,..., p. This does not, however, guarantee that f*(x) = f(x). Convergence in this latter sense may be tested by comparison of the results with (a) results obtained by other methods (b) results obtained using more accurate methods of numerical integration, and (c) solutions for idealized and limiting cases. Extension of the preceeding comments to the case of N+1 simultaneous integral equations is straightforward. The general form of the equations is N f, (X) I C X) + ( A* (XtX') he (' x' ( l,1, N) I-40 17

Advancing one step in the iteration procedure involves obtaining f.(n+l) (n+l) (n) f*(n) 0 N 0 N 2. Specific Numerical Techniques Used in Computation. A computer program was written to solve equations I-1 and I-2 in their general form by the method of simple iteration. The programs were written in a form suitable for use on the IBM 704 or 709 and are described in detail in Appendix C. All calculations, however, were performed on the IBM 704. Every effort was made to develop general techniques suitable for solution of equations with an arbitrary value of N; arbitrary set of constants in the phase function; an arbitrary method of numerical integration; any thickness, t1; any albedo for single scattering, o,; and any choice of starting functions. The same procedure was used for all values of m and N to obtain all of the functions presented in this report. The only limitation on the size of the problem that can be solved is that of computer time required and number of storage locations available. In practice, almost all of the results were obtained using five point Gausian quadrature. For this case, using numerical values given by Lowan, Davids, and Levenson'7 W0 = 0.11846 = 0.04691 w1 = 0.23931 kt = 0.23077 w2 = 0.28444 >2 0.50000 W3 = 0.23931 = 0.76923 W = 0.11846 "4 = 0.95309 18

In addition some values were obtained for comparison using a 10-point Simpson's Rule technique. The starting functions used for all calculations were,(t) o) P(() I-41a 9 t~r A) ~ A~n (iY) I-41b After considerable experimenting, it was found that for m f 0 a simple numerical test was satisfactory to determine convergence. In other words, when MX +*(n (i) - f e(X0) < Tolerance I-42a the iteration was stopped. For all cases the tolerance was chosen as 0.001. For m = 0, however, it was found that a better test was to stop the iteration whenever [:'. (X;). ")(X,])[ W' (Xi) "f (*'.X;)] O I-42b for' or [ or whenever I-42a was satisfied. This effectively stops the iteration as soon as the incremental change in any function reverses sign. This gave the best agreement with known solutions for (8) special cases (i.e. the X(4) and Y(i) functions of Chandrasekhar which are the same as o~0(r) and ~[(Q) in this study), 19

Some difficulty occurs when evaluating c (rC, p,') for the case p = A'....._ [= I- + _ 1I-43 However, by L'Hospital's rule the limit is C., (2,,j,?) q.,.. d O),, dA (M) I-44 film ~ t dr' I d,4 In practice, H and d — were evaluated at a given value of t by fitting a second order polynomial to the three nearest values of the function, and differentiating analytically. In addition, a program was written for calculating the integrated diffuse reflection and transmission for any set of functions using equations I-23b and I-25b. A discussion of the effectiveness of the numerical procedures will be deferred to section C-5 of this report after the presentation of numerical results. 20

C. Solutions to Specific Problems 1. Choice of Independent Variables. The choice of specific values of the independent variables was based on the following: 1. To indicate quantitatively the effect of the independent variables for representative conditions. 2. To gain experience in techniques of numerical solution of the problems which would be useful in solving further problems. 3. To obtain solutions for specific problems which can be used as a standard for evaluating approximate techniques. In order to completely specify the problem the following parameters must be chosen: 1) number of terms in the phase function, N; 2) coefficients of the angular distribution function, ao, al,...aN; 3) albedo for single scattering,.; and 4) the optical thickness of the dispersion, T1. In order to obtain a reasonable picture of the effect of co. and T1 a minimum of three values of each variable must be used, making a total of at least nine solutions for each phase function. There is an almost endless number of possible phase functions. For instance, for only a two term phase function (a0 = 1) using 10 values of a1 and a2 would require 900 solutions. For a five term phase function, 900,000 solutions would be required. In addition, for a given value of m, the effort in solving one set of integral equations increases with the square of the number of 21

equations in a set. Also, the number of sets of equations increases with N. For example, in comparable units of computer time, the effort in solving the equations for N = 0, 1, 2, 3, 4, and 5 is N =0 2 (1) = 1 unit N =1 (1)2 + (2)2 = 5 units N =2 (1)2 + (2)2 + (3)2 = 14 units N = 3 (1)2 + (2)2 + (3)2 + (4)2 = 30 units N =4 (1)2 + (2)2 + (3)2 + (4)2 + (5)2 = 55 units N= 5 (1)2 + (2)2 + (3)2 + (4)2 + (5)2 + (6)2 = 91 units (Note: in the preceeding it is assumed that all coefficients in the phase function are non-zeroin estimating the labor required to solve a set of equations N may be taken as the number of non-zero terms in the phase function) Thus, it is apparent that the selection of representative phase functions is very important and that the number of non-zero terms in the phase function must be kept small in order to be able to vary 0o and T1 sufficiently to achieve the previously stated objectives. 22

a) Choice of the Phase Function The construction of a phase function for scattering by spherical particles depends on the wave length A of the radiation, the diameter of the scattering particle D and the index of refraction m. (9) The problem is fully described by Chu and Churchill In summary, the phase function, p(cos Q) or p(i), can be written ~?) - 1 + I anPc(?) I-45 n I where P(?)/4IT is the fraction of randomly polarized radiation scattered by a spherical particle into a unit solid angle in the direction G. Q is cos-1( A) an( e,m) angular distribution coefficients Pn(4) are Legendre polynomials A aT D/A The phase function (I-45) is for non-aborbing particles ( o = 1.0). For partially absorbing particles with complex indices of refraction, the angular distribution coefficients must be recomputed for each value of W.. As a good approximation, however, the phase function may be considered to be WOp (4) where p (~) is the phase function of a non-absorbing particle. Probably the most important single characteristic of any given phase function is the fraction scattered forward and backward, F and B. 23

where F= 2nrtpCpdI 1-46 B = -Tr tp(') d I-47 B= I - F. Since r7 Fp()drl I-48 A good measure of how "peaked" is the phase function is the ratio of the second moment to the zeroth moment in the forward and backward direction p = Z'it t ts(ptrd/ F I-49 and p 2 T p(t') J'/ B-50 where P = Peakedness in Forward Direction. PB = Peakedness in Backward Direction. In addition, other characteristics are FP = Height of the Forward Peak, i.e. p(1.0) BP = Height of the Backward Peak, i.e. p(-1.0) The equations for these parameters in terms of the distribution coefficients, ao,..., aN are F I a a a; Y7+ _ 21a,, al + 3. I-51 2 + 16 32; z-6 5 2 Zo46 4-oq B =1F 1-48? = + T + -T - a +' S.'II + L,)s] I-52 _~t~a, a s azt,a, & a, ~)3-1/ 8 i1s W 3B+ 12SO 307z 6114+ 32716 -[ -1Aa a, a a al 3, a-53 PB = +6 - -I-iF 4 + 3S4- 1250 +3 +.I-53 I 3012 614+ 3116 - B 24

FP = j t a, + a2 4 t I-54 BP = 1 - a, t a 3 - 3 + a - I-55 For N = 0 the phase function has only one term, aO = 1. Thus, all quantities are fixed F = 1/2 B = 1/2 P = 1/3 PB = 1/3 FP = BP = 1 For N = 1 the phase function is completely specified by choosing any one characteristic auch as F. For n = 2, the phase function is completely specified by choosing F and P (or any other pair of characteristics). As more terms are added, more characteristics of the function may be choosen independently. In order to consider a significant range of (, and x1 with the computer time available a two term phase function was used for all computations including a few special cases in which some of the coefficients were zero. A careful study was made of all possible two term phase functions satisfying the conditions that p(,u) be positive for all values of 1t and that FP be greater than BP. It was found that the phase function p(,)- 0 + a,P,Q) + -,PI(p)]

gave the greatest range of F and P while satisfying the conditions that p(i) o for 0 < i ~ 1 and FP, BP. Figure I-3 is a plot of FP vs F with P and BP as parameters. The shaded area represents all values of F and FP for which p(i) was always positive. From consideration of Figure I-3, the following phase functions were selected for computation. TABLE I-1 Phase Functions Used in Computation a 1 a2 F P Remarks 0.00000 0.00000 0.50000 0.33333 Isotropic Scattering 0.00000 0.50000 0.50000 0.40000 Rayleigh Scattering 1.00000 1.81650 0.75000 0.5504 o Arbitrary: same F, different P 1.00000 1.04904 0.75000 0.4821 1.00000 0.00000 0.75000 0.3899 1.73205 1.00000 0.9330 0.4821 Largest F for N = 2 1.73458 2.24153 0.9336 0.5821 Same Phase Function as Reference (5). The last phase function listed in Table I-1 is the same as that in reference (5). Computations with it were made only for 1 = aD, since p(i) was not always positive for this function. The results were used as a check on the accuracy of the numerical techniques. All of the phase functions are plotted in figures 4a, 4b, and 4 c. In these figures the function f(Q) is plotted which is equivalent to p(cos.)/41( 26

8.0 7.0 Il 60/ / -00'/'/ 6.0 -9950 aPBar o.0- ~69, 3.02 _I ___.Oo _-0-I~FP, FPC BPand 2.0,!P0_- / FP BP a,FP * / ~~ =0.3889 line on the top. 1. 0 otropc FPShaded area:' p(M)always positive Pl/ 0 F2 + I/ /4 27 0~,-5 0.00.008009010

0.20 4a | - Rayleigh 0. 1 5 -- - Isotropic f(6) O. I _ 0.05. l i. 0.00..... 0 20 40 60 80 100 120 140 160 180 f ( @ )0.2 5 __'_ 1 n —=9 P =0.4821 0.40 ____ __2568 0.35 0...889 0.2 N 0. P0 0.1 0 0.15 0. 00 000 20 40 60 80 100 120 140 160 180 0.40, 0 F 0.35a 09=0.93360, t0.5821.. 0. 35 ~0 - FO. =9330, P=0.4821 0.25 0.20 0.115 O. lO 0.05 0.00AV, AP I. -'-, — ~ 28

b) Choice of Other Independent Variables The independent variables other than the phase function (i.e. Wo and wl) were selected to give as good an indication of their effect as possible. For most of the calculations Ao was chosen as 0.90. This was done to be able to compare results for T1 = oo with previous solutions. For every phase function results were obtained at least for T1 = 1.00. For isotropic scattering, Rayleigh scattering and for the function corresponding to F = 0.9330 and P = 0.4821, results were obtained at several values of T1. A total of 36 problems were solved; the values of the parameters for each are listed in Table I-2. 29

Table I-2 Values of the Parameters Used in Computation Problem Number F P 1 di 1 0.500 0.333 0.05 0.90go 2 0.500 0.333 0.25 0.go 3 0.500 0.333 0.50 0.90 4 0.500 0.333 1.00 0.go 5 0.500 0.333 1.50 0.go 6 0.500 0.333 2.00 0.go 7 0.500 0.333 o0 0.90 8 0.500 0.333 0.05 0.60 9 0.500 0.333 0.50 0.60 10 0.500 0.333 00 0.60 11 0.500 0.333 0.05 0.30 12 0.500 0.333 0.50 0.30 13 0.500 0.333 00 0.30 14 0.500 o.400 0.05 o. go90 15 0.500 o.400 0.25 o.90go 16 0.500 0.400 0.50 0.go 17 0.500 o.400 1.00 o.go 18 0.500 o.400 1.50 o.90 19 0.500 o.400 2.00 0. 90 20 0.500 o.400 oo o. 9go 21 0.750 0.5504 1.00 0.90 22 0.750 0.4921 1.00 0.90 23 0.750 0.3889 1.00 0.90 24 0.9330 0.4821 0.05 0.90 25 0.9330 0.4921 0.25 0.90 26 0.9330 0.4921 0.50 0.90 27 0.9330 0.4821 1.00 0.90 28 0.9330 0.4821 2.00 0.90 29 0.9330 0.4821 0 o0.90 30 0.9330 0.4821 0.05 o.6o 31 0.9330 0.4821 0.50 o.6o 32 0.9330 0.4821 00 o0.60 33 0.9330 0.4821 0.05 0.30 34 0.9330 0.4821 0.50 0.30 35 0. 9330 0.4821 00 0.30 36 0.9336 0.5821 oo O..90go 30

2. The r and 0 Functions For each set of parameters as indicated in Table 1-2 the complete set of 1m and 01 functions are calculated at the 1 1 Gaussian quadrature points. These functions are tabulated in Appendix B. The number of iterations required to evaluate the functions for each value of m is given. In every case except Problem 36 for m = 0, the starting functions are those given by Equations I-41a and I-41b. 31

3. Sample Values of the Intensity of Reflected and Transmitted Radiation The normalized intensity of reflected radiation and the diffuse portion of the transmitted radiation are, in general, a function of two angles 0 and 0 and four parameters: Go, 0o' T1, and the phase function. The effect of the variables is illustrated by considering a sample problem with 1) -0 = cos-1 0.5 (i.e. %0 = 0.50) 2) Q = cos-1 0.76923 (i.e. p. = 0.76923) 3) Co = 0.90 4) X! = 1.00 5) A two term phase function such that F = 0.9330 and P = 0.4821 The angular distribution of the specific intensity is obtained as a function of the bearing angle ( as well as these six parameters. The effect of each variable is indicated by varying QO, Q, ), r1, and the phase function, one at a time, while holding the remaining variables at the above, "standard values". (There are exceptions which are noted.) The specific intensities for each case were calculated on a desk calculator using I-21 or 1-22 and the equations in Appendix A-1. The apparently arbitrary values of p and lo which were used correspond to the Gaussian points for the five-point Gaussian quadrature utilized for integration. The values of the coefficients bo, bl, b2, co, cl, 32

and c2 which appear in equations I-21 and 1-22 are tabulated in Table I-3. The effect of Q0 (or 4o) is shown in Figures I-5a and I-5b. In this case [ was held at 0.50 and four different values of 0o were used. It is observed that there is a rather complicated dependence of the intensity of reflected radiation on /Pe, with some of the curves overlapping. At (p = 1800 the maximum intensity of reflected radiation occurs for So between 0.23079 and 0.76923 which is not surprising since the cosine of the angle of incidence is 0.50. The effect of 0 (or A) is shown in Figures I-6a and I-6b with all other variables at the "standard" values. In this case the dependence of the intensity of both reflected and transmitted radiation on p produce overlapping curves. The maximum intensity of reflected radiation at r = 1800 apparently occurs at a viewing angle almost the same as the angle of incidence. The effect of "0 is shown in Figures I-7a and I-7b. In this case T1 is 0.50. It is apparent that decreasing W. simply decreases the intensity of the reflected and diffuse portion of the transmitted radiation without affecting the nature of the dependence on (. The effect of T1 is shown in Figures I-8a and I-8b. The effect of decreasing T1 is similar to that of decreasing 0,: the intensity is decreased but the dependence on 9 is not changed very much. 33

Table I-3 Values of the Coefficients Used in Equations I-21 and I-22 to Calculate the Specific Intensities of Reflected and Transmitted Radiation. F P bw r, b b~ b 2 Co C1 c 0.500 0.333 0.90 1.0 0.5000 o.o469 0.1033 0.0000 0.0000 0.0326 0.0000 o.oooc "t 1 IIt It It 0.7692 0.0560.o0000.o0000 0.0429.o0000 o.OOOC "t It It 0.25 " I 0.0219 0.0000 0.0000 0.0174 0.0000 0o.00C "I IT 0.60 0.50 i 0.0205 0.0000 0.0000 0.0183 0.0000 0.000C' o0.400 0.90 1.0 " 0.0o469 0.0691 0.0020 0.0765 0.0268 0.0035 0.013C It IT i" t it 0.7692 0.0485 0.0091 0.0133 0.0388 0.0061 0.0091 0.750 0.550 I" It I 0.0078 0.0186 0.0581 0.0348 0.0443 O.042C t 0.482 t" I t 0.0247 0.0014 0.0300 o0.0o459 0.0285 0.020S It o.389 " I' IT 0.0423 0.0195 0.0000 0.0571 0.0131 0o.oooc 0.933 0.482 " "t It 0.0367 0.0225 0.0297 0.0000 0.0000 0.000C It it "i 1.0 " o.o469 0.0214 0.1189 0.1597 0.0304 0.0329 0.028[ "' " I " "t 0.2308 0.0230 o.o845 0.1163 0.0432 0.0502 0.040 II'I iIt " " 0.7692 0.0132 0.0201 0.0284 0.0601 0.0467 0.0197?' II WIt " " 0.9531 0.0126.o63 0.0055 0.0719 0.0245 o.oo4c t I I " 0o.o469 0.5000 0.0020 0.0112 0.0150 0.0028 0.0031 0.0027 I it it 0.2308 t 0.0106 0.0390 0.0537 0.0197 0.0232 0.018E " IT IT I"t 0.7692 " 0.0204 0.0310 0.0436 0.0924 0.0718 0.0301 "I It It It 0.9531 it 0.0241 0.0120 0.0105 0.1371 o.o467 0.0077 it It IT 0.50 0.5000 0.7692 0.0044 0.0142 0.0232 0.0392 0.0425 0.021C it it 1 o0.25 IT i7 0.0003 0.0076 0.0157 0.0245 0.0298 o.o015C it iIt o.60 0.50 ot 0.0007 0.0076 0.0149 0.0218 0.0254 0.013~ It It 0.30 " It o.ooo6 0.0030 0.0072 o.oo93 0.0117 0.00o6 34

0.1 6 5a 0 0.1 4 Z 0.1 2 u 0.10 w 0.0 8 WL 0 0.0 6 Z I 0.0 2_ 0 20 40 60 80 100 120 140 160 180 0.22 I 5b 0.2 0 i For both figures: 00.~18 _~' /~~= 0.04691 0 0.1 8 ------ /,0=0.23077,- 01O =0.76 923 0Z 6 "IN C, O=O.95309 O un 0.14,~ 0.12 t 0.08 0.06. 0.04 z 0.02 _ 0.00 0 20 40 60 80 100 120 140 160 180 Figure I-5 Intensity of reflected and transmitted radiation as function of bearing angle, I/I0 vs. f, parameters of Io.) - = 0.50 c (3= 0.90,'1 = 1.00, F = 0.9330, and 35

0.30 6a 6a 0.28 0.26 __ For both figures: _ = O. 04691 0.24 - / - 0.23077 ~~~0.22 ____/~C -= 0. 76923 0.2.. = O. 95309 0~~~~~,- 0.20 ~ 0.18 0 O 0.16 W w 0.14 LL 0. 12 0 0.1 0 z w 0.08 z 0.06 0.04 0.02 0.Q0 20 40 60 80 100 120 140 160 180 0.14 C0.1Q2 -w O, 0.062 LL 0 0. 06 0.04 z - 0.02 Z 0.00 0 20 40 60 80 100 120 140 160 180 Figure I-6 Intensity of reflected and transmitted radiation as function of bearing angle, I/Io vs. f, parameters of C~ ~o = 0.50, o3o= 0.90, 71 = 1.00, F = 0.9330, and P = o.4821

0.08 o 7a'-I 0.07 I-4 z 0.06 0 I — U 0.05 J I w 0.90 w 0.04 b. o 0.035 o,= 0.60 0.02 (f) w 0.01 0.00 0 20 40 60 80 100 120 140 160 180 0.1 I wo=0.90 7b 0.1 0 c — 0.09 0.08' 0.07 0 06 =0.60 0.06 o 0.04 Vi 0.03 w00.30 z 0.02 z~ 0.01 0.00. 0 20 40 60 80 100 120 140 160 180 If Figure I-7 Intensity of reflected and transmitted radiation as function of bearing angle, I/I0 vs. e, parameters of r3o ~o = 0.50, i = 0.76923, T1 = 0.50, F = 0.9330, P = 0.4821 37

0.09 oa 0.08,0.07 ~~~~~~~o N'. 0.06 z 1_ 0.05.J0.04( ~.,, - 0. IJ La' 0.04 w LL 0.03 \ >- 0.02 H C,) Z 0.01 w z 0.00 0 20 40 60 80 100 120 140 160 180 0.1 3 8b 0.12 0.1 I l=1.0 010 I 0.09 o 0.08 C, 0.07 C,) z 0.06 0.05 0 >- 0.04 C, Z 0,03 w Z 0.02 0.01 0.00 0 20 40 60 80 100 120 140 160 180 Figure I-8 Intensity of reflected and transmitted radiation as function of bearing angle, I/I0 vs. f, parameters of'r1, ko = 0.50, Ci = 0.76923, co= 0.90, F = 0.9330, and P = 0.4821

The effect of phase function is indicated in Figures I-9a and I-9b. Results are given for isotropic and Rayleigh scattering in which F is the same (0.50), but P is different (0.3333 and 0.40). An additional set of results was obtained for 4 = 0.04691. It can be seen that the two phase functions produce approximately the same angular distribution of intensity except for reflected radiation with a near grazing viewing angle (I = 0.04691) in which the "dumbell" characteristic of the Rayleigh phase function is accentuated. Another illustration of the effect of phase function is provided in the polar plots of Figures I-lOa and I-lOb. In this case the three phase functions for F = 0.750 and P = 0.5504, 0.4821, and 0.3899 are used with the "standard" values of the other variables. It can be seen that the different values of P at the same value of F produce noticeably different angular distributions for the reflected and transmitted intensity. 39

0.15 90 0.13 0.1 32 0.12 /L=0.04691 00.11 H H 0.10 _o 0.09 w cr 0.07__ _ _ _ LI r 0.06 ____ ___ "' 0.05 LI.. z LU.: 0.04 _ —0.03 >- u.0 = O 76923~ 0.05 2 0.04] 0.03. 0~~...... Isotropic 0.0 I - Royleigh 0.00... 0 20 40 60 80 100 120 140 160 180 0.06 0.06__9 b____ 0 9b,) 0.05 _ v) L=0.76923 cn 0.04 z cr 0.03 -. H p.=0.04691 LU.. O~~~ 0.02 1- ----— l~-Isotropic U 0.01,, Rayleigh z w 0.00 z 0 20 40 60 80 100 120 140 160 180 Figure I-9 Intensity of reflected and transmitted radiation as function of bearing angle, I/Io vs.,, isotropic and Rayleigh scattering, i = 0.04691 and 0.76925, po = 0.50, o3o= 0.90, t1 = 1.00 40

I a 600 /200 900 / / F=0.750,,wo=0.90,T/,=tO 9 0o =0.50, =076923 /50 4 -180 ~ 0,, ~ ~.2-.02-. 4 06-.08-./0 -./2 - 0o 270 = 0.76923, o= 0. 9001 =.00 150 4

4. The Integrated Reflection and Transmission The integrated reflectance, the diffuse portion of the integrated transmission, and the total transmission were calculated for each set of conditions and the results are tabulated with the < and Al functions in Appendix B. In addition, these values are presented graphically in Figures I-11 to I-25. The integrated reflectance for isotropic scattering is shown as a function of angle of incidence in Figures I-11, I-12a, and I-12b for t. of 0.90, 0.60, and 0.30 with 1 as a parameter. Similar plots are presented for anisotropic scattering with F = 0.9330 and P = 0.4821 in Figures 1-13, I-14a, and I-14b. The total integrated transmission, including both diffuse and direct components, for isotropic scattering is plotted vs. So in Figure I-15 for W. = 0.90. In Figures I-16 the total transmission for isotropic scattering is presented for T1 = 0.05 and 0.50 with parameters of O.. Similar plots are given for anisotropic scattering in Figures I-17 and 1-18. Values of the integrated reflectance and transmission were also obtained for Rayleigh scattering for which F is the same (0.50) as for isotropic scattering, but P is different (0.40 compared with 0.333). The results as shown in Figures I-l1 and I-15 are almost indistinguishable from those for isotropic scattering. Integrated reflectances for 42

0.70 = 0.90 0 Isotropic scattering A Rayleigh scattering 0.60 0.60 a t~~~~ SNot indicated when very close to isotropic case 0.50 L) ~- ~0 w 0 a:_ 0.30 LLCD LLI 0.20 Or 0. 09 43

0.40 12a 0nr"~~~120 uu W. = 0.60 t0 n' 0.20 W 0.10 = 0.05 0 0.20 0.40 0.60 0.80 1,00 0 cc 0. 30 04',W 12b - o =0.30 IW _0.20 _ Isotropic scattering LL 0::: <z ~ ~ ~ ~._... 005 0 0,20 0.40 0.60 0.80 1.00 pMo Figure 1-12 Integrated reflectance as function of angle of incidence for isotropic scattering, R vs. ~o, parameters of ~1,'io = 0.60 and 0.30 44

0.70 0.60 CWo =0.90, F=0.933 0 P= 0.482 1, This investi - gation CrI 1k \\'\1 0 P=0.5821 WL 0.50 Reference (S) 0 w0 0.40 0Q~0 ( 0.30 0.20 \ x 0 0 = 0.290, F =0.94 0 0.60 0.80 1.004821 Figure 4-13 Integrated reflectance as function of angle of incidence for anisotropic scattering. R vs. ~oi parameters of Ti~ 45

0.40 O 14 a 0.30 (g i wo= 0.60 LU -J F = 0.9330 LL LU n- 0.20 04 P=0.4821 0 0 0.20 0.40 0.60 0.80 00 2 I I 14b H 0.20 UF=0.9330 Cr LU 0 0.20 0.40 0.60 0.80 1100 /Lo Figure I-14 Integrated reflectance as function of angle of incidence for anisotropic scattering, R vs. L0,O parameters of 71y 46

0.90 0.80 0.70 o 0.60 0.20 Ico = 0.90 Isotropic scattering icdnef0.r iorp10i and Rayleigh scattering Not indicated when very close to isotropic case 0 0Q20 040 0.60 0.80 1.00 0.50 ~ ~ 4

0.90O..r=0.05 z-~.=0.50 0.80 0.70.,o H ~ ~ ~ ~ ~ -- J 0.60 0... (/) //0.60 -. coII'll- / / 1 0 ~30/,, o.0. / d' LU 0.50 AH / ~~~~~~// Q~~~~ CE o/__ __ < / / 040/ / H~ / // 0 O~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 0.30 0.10 /d 0 0 0.20 0. 40 0.60 0.80 1.00 1-o Figure 1-16 Total integrated transmission as function of angle of incidence for isotropic scattering, T vs. Poe parameters of uo, T1 = 0.05 and 0.50 48

6t1 TIVo = J'o0~6o = a'o6'o =o06'To jo saaq0amTiavd c'BA JG'I UTa~q.;SDS odoaqosius Ioj aourapToui JO GTCuv JO UIOIIouUJ Ss UOTSSLISuuzPV paleLUBWWT. Tu-O~O Li-I GanD1Jd ~W 00'1 0_ 0 09'0 0~'0 OZO 0 0 0~'0 -H ~~~~~~~~~~~~o'o — I -~~~~~~~~~~~~~ > 0 H 18t7'0= d 09'0 H "Y ~~~z o~'o-~'/ / /H~ O ~~~~ 6'0 = j O~OL'.0 06'0 = Om b~~~06 — I 001 Z ~~~~C3 ~ ~ ~ ~ ~ O 09'O ~~~~I/ I 0 I /J ~~' 0'0 -,~~~~~~~~ 00'1

IZ8g'0 = d C'oE6'o =' 05O P0 puB 0'0 = L ~)3 jo sLaG[~%ead ~Tt *SA'TsG I, uIa s oDdoa;osuL aoj aouGapTDuT Jo sITuB jo uotoiunj se uoTss'msu~[; psG;DaluIT Tieow,, 81-I ianI,2 Or7 00'1 08'0 090 0" 0 0 Z'0 0'! I tv; 00 / 01'0 /.. / / 9 / / 09'0 / / P p 0-,/ / 0'0 / H0 /' /- /i /! /.1' /c.z / /.0' O

00 - = 1 pu o606'0 =oC)'j jO sla9acU d R l *T SA { auGpTDoUT Jo aGlTu Jo uoiTounJ se aoun;DGTJGa pGa%;aDnuI 61-I alnD!: 00'1 08'0 09 0 O O 0 ZO O 0 z -I m Ot,*O 0'0 m z 68~'0= d V IZ'B8'0: d 0 _o t0cgg0=d 0:aAJno GIPPIVI 09'0 0'1=' 06'0 = Om

0.80 0.70... 00/ o 0.50 z 0.40 03 w=0.90.r, 1I.00 0.20 Middle curve: o p=05504 o p = 0.4821 A P = 0.3889 0.10 | Bottom curve: o Isotropic _a Rayleigh 0 0 0.20 0.40 0.60 0.80 1,00 Figure I-20 Total integrated transmission as function of angle of incidence, T vs. HO parameters of F,.= 0.90 and - 1.00 52

0.8 0 Ii 21 a W"= 0.90 _sotropic or Rayleigh scattering z " 0.6 - =0046 U~~~~~~~~~~'09' 0.2 T Lw H~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ZD j LU r- 0.2 0 i z 0 0 0.25 0.50 0.75 1.00 0.75 050 0.25 0 L. _i____________........... z' -, —- -— _~'n~ T,~ ~ ~ - 1.0 21 b wo= 0.90 Isotropic or Rayleigh scottering o 0.8 U) 0~~~~~~ HO9. 0.6 LU02 0; rr' Z 06"0~ 0 0 — 1~~- 0o-~69/ - H~~~~~~~~ _ 0.O Q o~5 O 0 0.25 0.50 0.75 1.00 0.75 0.50 0.25 0 I_~~~r Ti L I _! ~ —— ~I- -- T Figure 1-21 Integrated reflectance and total integrated transmission as function of optical thickness of dispersion, R and T va. Ti~ parameters of Po isotropic or Rayleigh scattering F = 0.50 53

TOTAL INTEGRATED TRANSMISSION,T INTEGRATED REFLECTANCER ~It 0 0 0 9 9 0 0 0 H'. 0 b:0D 0 0 N o C d O- 0 (iD 0.. c+ H(D H iu O 3d 0CD 1o CO~0 C C+ c H 0' FJ (.7 ~' I 1 o~ c t, I (D D 0 0 o ct+D H,7 o o 0~~ ~~~o 03 H~ b 0 0 H' Ct.0- 00 01~~~ 9 m (D C c+~-.C NO III ( 0ct~~~~ ct~~~ ~ -0 0o CDe H'P I ~ CtcD Oa o - 0H 1-10 _o_ 0 Id i-b J F-J- O co C/ o ~~~ ~ ~ ~ o W (D c+ (D FJ- 0 W 0 0 r U Qj cn 0 O F.o ~ o (z) v, O~~~~~~~~~~~~0 ~~~~~u ~~~~~~~~~~~~~~ ~~~~~~~CX) -- 00 OD 0~~~~~~~~~~~~~~~~~~~~~~~~~~~~0 L~~~lO~,, I I I I\ \ I I I'F- 0 - -- ~~~~~~~~~06

I SOTROPIC 1.0 1.0 1.0.0,=0.04691 / LO=0.5000 1,u O=0.95309 0.8 0.8 0.8 0.6 0.6 0.6 nr 1,59, ~o I:::~ — 0.4 04 040.2 0.2 0.2 0 0.5 1.0 0 0.5 1.0 0 0.5 1.0 Wo Wo W ANISOTROPIC 1.0 d 1.0 1.0 d e f /Lo=0.04691 /1o0=0.5000 1 o=0.95309 0.8 0.8 0.8 0.6 -0.6- 0.6 0 0.4 04 02 5 0.2 0 0.5 1.0 0 0.5 1.0 0 0.5 =0.01.0 0 0 0 Figure I-23 Integrated reflectance as function of albedo for single scattering, R vs. 00, parameters of F, ~,o and 1 55

ISOTROPIC r a 1l~lb c 1.0 1.0 Ib' =0.05 0.8 004691 0.8. 0.8 Lo:0.04691 /'Lo=000 /0=0.9530 9 0,6 0.6 0.4 0.4 0-.4 0.2 0,2 0~2 0 0 - 0 0 0.5 1.0 0 0.5 1.0 0 0.5 1.0 Wo W ao Wo ANISOTROPIC d e If 0.05 1.0 0.0O5 1.0 L0= 0,04691 U, = 0.5000 /.L =0.95309 0.8 0.8 0.8 0 0.6 0.6 - 0.6 04' 0.4i / 0.4 1 0.4 0~ 0.2 $ 0.2 0.2 O l i o I O I 0 0.5 1.0 0 0.5 1.0 0 0.5 1.0,o Wo Wo Figure I-24 Total integrated transmission as function of albedo for single scattering, T vs..o, parameters of F, 4o,1'

25a 0.802 0.80 1 G=w.e 0.90, TI 1.0 0.60 bI-0=0.04691 w0 /o 0.23077 -J I_ C 0.40 0 /=::Q6922 0.20 0 I I I 0.50 0.60 0.70 0.80 0.90 1.00 z 25 b z 0.80 <oo=.90, T, = 1.0 I — t 0.40 Io i 0.50 0.60 0.70 0.80 0.90 1.00 Figure I-25 Integrated reflectance and total integrated transmission as function of phase function, R and T vs. F, parameters of To' 0.9= 09,0 1 = 1.00 57

the anisotropic phase function of reference (5) for ~1 = co in which F was 0.9336 fall on the same curve as the results for this investigation in which F is almost the same (0.9330), but P is different (0.4821 compared with 0.5821). The effect of anisotropy is best indicated in Figures I-19 and I-20 where the integrated reflectance and the total transmission are plotted versus angle of incidence with 1 = 1.00 for each of the five phase functions considered in this investigation. Again, it is noted that the results for the three different phase functions with F = 0.750, but with different values of P are almost identical. The integrated reflectance and total transmission are cross-plotted versus T1 with 40 as a parameter for isotropic scattering in Figures I-21a and I-21b and for anisotropic scattering in Figures I-22a and I-22b. This should provide a satisfactory method of interpolating for different values of 71. As is expected, the dependence on T1 is almost linear for small values of T1. The effect of wo on the integrated reflectance is indicated in Figures I-23a and I-23f for isotropic and anisotropic scattering. The effect of t, on the total transmission is shown in the same manner in Figures I-24a and I-24f. The dependence on Qo. is almost linear for small values of 4.. For small values of T1 it is linear over the whole range of O..

The effect of phase function —in particular, the effect of F-is illustrated by cross plots of R and T versus F in Figures I-25a and I-25b where 0o is a parameter. This indicates the approximate error of neglecting anisotropy and provides a means of estimating the reflection and transmission for an arbitrary dispersion on the basis of the value of F for the phase function. As can be seen from Figures I-19 and 1-20, the effect of P is not large. In fact, except for the integrated transmission at the two largest angles of incidence, the complete variation of reflectance and transmission with the range of P considered in this study is indicated by the size of the circles used to plot the points in Figures I-25a and I-25b. 59

5. Discussion of Numerical Technique It is very difficult to estimate the error involved in the numerical solution of the integral equations. Probably the best indication is comparison with previous solutions and with other methods of solution. The following solutions are available for comparison. 1. Chandrasekhar's H functions (l) for, = 0.90, 0.60, and 0.30. The H(CI) functions are the same as the ~([) functions of this report for 1 = oo and isotropic scattering. The comparison is given in Table I-4. The values generally agree to within 0.1%. 2. Chandrasekhar's X and Y functions (8) for 0o= 0.90 and T1 = 0.05, 0.25, 0.50, and 1.00. The X(>) and Y(Ci) are the same as the,(~) and 0o(i) functions of this report for isotropic scattering. The comparison is given in Table I-5. The values generally agree to within a few percent. 3. The o(t), *(0), f(Q), 41(p), *(p) and 2(p) functions calculated by Churchill, et.al.(5) for a0 = 1.000000, a1 = 1.73458, a2 = 2.24153, wo= 0.90 and T1 = oo.with iteration using a 20-point Simpson's rule technique. The comparison is given in Table I-6. The values generally agree to within 1%. 4. The integrated reflectance calculated for Rayleigh scattering, SO= 0.90 and T1 = oo using the method described in Appendix A-2. The comparison is given in Table I-7. The agreement is very good. 60

TABLE I-4 Comparison of io8( i) Functions for Isotropic Scattering and with Chandrasekhar's H Functions. () oo0.90 o.6o 0.30 C1 H(4) 00(4) H(4) 0~([1) H(0) O( ) 0.04691 1.0937 1.0936 1.0520 1.0520 1.0233 1.0233 0.23077 1.3233 1.3236 1.1584 1.1587 1.0661 1.0663 0.50000 1.5560 1.5554 1.2485 1.2456 1.0976 1.0975 0.76923 1.7299 1.7289 1.3011 1.3009 1.1160 1.1159 0.95309 1.8273 1.8259 1.3291 1.3288 1.1249 1.1248 Table I-5 Comparison of thefo(G) and 0o(4) Functions for isotropic Scattering and 4o= 0.90 with Chandrasekhar's X and Y Functions. T1 = 0.05 x ~0'o 0.04691 1.0560 1.0546 0.3914 0.3838 0.23077 1.0778 1.0726 0.8805 0.9049 0.50000 1.0814 1.0783 0.9852 0.9824 0.76923 1.0826 l.0800 1.0190 1.0138 0.95309 1.0831 1.0806 1.0314 1.0252 T1 = 0.25 x IX Y 00 x 0 0.04691 1.0810 1.0797 0.0441 0.0467 0.23077 1.1967 1.1993 0.5017 0.4990 0.50000 1.2440 1.2413 0.8285 0.8317 0.76923 1.2610 1.2611 0.9672 0.9501 0.95309 1.2668 1.2694 1.0228 1.0001 T1 = 0.50,, x Y 0 0.04691 1.0861 1.0842 0.0269 0.0275 0.23077 1.2553 1.2530 0.2769 0.2636 0.50000 0.3528 1.3465 0.6506 0.6497 0.76923 1.3948 1.3896 0.8636 0.8462 0.95309 1.4122 1.4089 o.g9583 0.9249 T = 1.00 I. X I, y 0 0.04691 1.0902 1.0888 0.0171 0.0165 0.23077 1.2933 1.2915 0.1888 0.1231 0. 50000 1.4574 1.4541 0.4120 0.4012 0.76923 1.5472 1.5427 0.6646 0.6476 0.95309 1.5880 1.5852 0.8000 0.7553 61

Table I-6 Comparison of the r and 0 Functions of this Report with Similar Functions of Reference (4) Calculated Using a 20-point Simpson's Rule Technique for a1 = 1.73458, a2 = 2.24153, Wo = 0.90 and T1 = co. Ref. This Rept. Ref. This Rept. 0.04691 1.1186 1.1185 0.0187 0.0185 0.23077 1.3538 1.3550 0.1182 0.1179 0.50000 1.4865 1.4868 0.3121 0.3121 0.76923 1.4693 1.4699 0.5413 0.5412 0.95309 1.3946 1.3948 0.7094 0.7094 $2 Ref. This Rept. Ref. This Rept. 0.04691 0.5368 0.5370 1.0515 1.0516 0.23077 0.5110 0.5124 1.0888 1.o896 0.50000 0.2078 0.2079 0.9564 9.0561 0.76923 0.3782 0.3771 0.6578 0.6583 0.95309 0.9347 0.9348 0.2822 0.2905 1 2 K1 Ref. This Rept. Ref. This Rept. 0.04691 0.1084 0.1080 3.1399 3.1376 0.23077 0.5908 0.5911 3.1955 3.1955 0.50000 1.2591 1.2596 2.6409 2.6406 0.76923 1.5187 1.5179 1.4697 1.4695 o.g95309 0.8840 0.9173 0.3330 0.3331 Table I-7 Comparison of Integrated Reflectance for Rayleigh Scattering, 1 = oo and wo= 0.90 with Values Calculated using Method Described in Appendix A-2. K0 Appendix This Report 0.04691 0.660 0.656 0.23077 0.581 0.580 0.50000 0.507 0.507 0.76923 0.453 0.453 0.95309 0.424 0.424 62

It must be noted that in all the preceeding comparisons (except for the case of p = 0.50) the results from other investigators were linearly interpolated to get values for comparison with the results at the Gaussian quadrature points of this report. From the comparisons, however, it is estimated that the functions are generally accurate to within +0.02 which in most cases in an error of less than two to three percent. The intensities of reflected and transmitted radiation as well as the integrated reflectance and transmittance are of the same general accuracy as the functions. This is evident from equations 1-15, I-19, I-29 and 1-32. The use of Gaussian quadrature seems preferable to Simpson's Rule or to the trapezoidal rule in that five point quadrature afford the same accuracy as the use of ten equally spaced ordinates. In addition, a better approximation may be had by adding only one additional point where with Simpson's rule, in order to obtain a better approximation, it is necessary to double the number of points. The disadvantage of the use of quadrature is that results are obtained at odd values of.. On the other hand ^in using the results they are to be integrated again, this may be an advantage, since quadrature can be used to perform the integration. The time in seconds required for one iteration on the IBM 704 is Time = 0.025 (no. Points)2x(no. equations)2x(no. iterations) I-56 63

For example, using five-point quadrature; for m=0, N=2; performing four iterations Time = 0.025 (5)2(3)2(4) I-57 = 22.5 seconds Considering the cost of 704 time at $5.00 per minute Cost = t2.08 x 10-3 (no. points)2x(no. equations)2x(no. iterations) I-50 To scale the cost to some other machine, simply multiply the factor $2.08 x 10-3 by (cost per operation on new machine)/(cost per operation on the IBM 704). The 704 performs 40,000 operations on the average per second so its cost per operation is $5/(60x40,000) = $2.08 x 10-6 per operation It is estimated that the $2.08 x 10-3 factor in I-58 could be reduced by a factor of five through the use of efficient hand-programming. 64

PART II EVALUATION OF APPROXIMATE MODELS Exact solutions in closed form are not attainable by analytical methods for problems of multiple scattering. Exact solutions are attainable by numerical methods, as in Part I, but the required calculations are lengthy and expensive, and the results are for specific cases. Simple, approximate solutions therefore have considerable value if reliable. In this section the reflectance and transmission obtained from various approximate models are compared with the "exact" values obtained in Part I. A. Approximate Models 1. Two-flux Model(11, 12) The two-flux formulation of the slab problem with an obliquely incident plane-parallel source is'Jr - (I F) re + w, I: II-la,/,,dr2= (I- T.)cF)Is - II-lb Where Il is the integrated radiant flux in the forward direction and 12 the integrated radiant flux in the backward direction, and F and B are the forward and the backward components respectively of the phase function for single scattering.

Equations II-la and II-lb can be solved with the boundary conditions I 1 at Z = o, and I2 o at Z = T1 to give (Xi)z=o/4. (I -, F -pXI - e r ) II-2 R =: i-wI —2 I -,oF p e and T(=-(a ) Pz? =11.. _II-3 T (I,)zoo I4F - oF + - F - oFp) e_'Pt'/ where pJ-o 11 F) - ( II-4 2. Six-flux Model (11,12) The six-flux formulation of the slab problem with an obliquely incident parallel-plane source is Ja = - C, seceI, + C2 secej t C3 5eCeo( 3 + 1)5a dI2 -c, secOI1 + C SecG.x, + C SecO.(I3 + I) I-3. - C, c CS,3 + C cscI + C3 CSCe0(, + I,) I-) xr I ( + I t 13 + 14) II-6 I - to0 F6 - 0 66

Where Ii. I2, I3, I4, 15 and I6 are the components of the radiant flux in the forward, backward and sidewise directions, and c1 = |46 - F C2oX e II-7a ~CI~ - co Ft -0o, oF6 - 2X II-b Ca = So B6 + II-7b I-,F6 8- B 6 C Q X + 2 -Xz II-7c I - )F6- 0oB6 where F6, B6 and X are the six-flux representations of forward, backward and sidewise scattering components. F6 = (F)(P), B6 = (B)(P) X = (1 - F6 - B6)/4. The solution for Equations II-5'a to II- 6 with the boundary conditions I1 = 1 and I3 = o at Z= o, and I2 = 4 =o at Z = T1 is reproduced in Appendix A-3. 3. RicIard's Modified Diffusion Theory (13) The transport equation is approximated as dog _ 1(2g X _ ( 1_ K )t>j II-8 dz,. - ( 1-8(- K )ed where Z is normal distance in mean free paths for scattering, k = V3(CI e)o/(3.2'o); 0d is the density of unscattered photons, p is the photon density and q is defined by g= (i,- -o.)e - e'd II-9 The source strength is pd = e 67

Equation II-8 can be solved for the conditions (o = o in vacuum, p finite, and q and EJ both continuous at 1 = o and T1 Giving I -(X )e (Iz K)1". II-1Q Rz II-la ~ XZO ( +K )Ze'-(- K)'e-'' K and Tz-~I, -, + /, _ 21+.) ~-(,+K I + K)e +, - XI-K/)e -'/ K + e / - (!+ ~0, )e, i ( i - 2)-' I'-zIThe above two equations are for c.< 1. The solutions for the case Qo, 5 1 can be found easily and are omitted here. 4. Numerical Results Numerical values were computed from the above approximate solutions for comparison with the "exact" results obtained in Part I. The computed values are given in Tables B-1 to B-5 in Appendix B and are illustrated in Figures II-1 through II-5. Table B-1 gives reflectances for isotropic scattering as obtained from Richard's model, the two-flux model and the six-flux model. The exact values are included for comparison. The six-flux results for Rayleigh scattering are also included. The exact values for Rayleigh scattering did not vary significantly from the values for isotropic scattering and hence are not reproduced here. The two-flux model

0.70 __i_ 0.60 0. 10, - -_,! z j 0.40 U] 0.30 z 0.20 Isotropic scattering T1' = 1.0 A Exact.-. ~Six - flI-ux 0.10 —. —Two-flux v Richards' 0 0.20 0.40 0.60 0.80 1.00 Mo Figure II-1 Integrated reflectance for various solutions as functions of angle of incidence for isotropic scattering, R vs. ~o' 3o = 0.90, n1 = 1.0 69

0.70 0.60 rrz~~~~~~ t << T, = Co 0.5 -- z _1 GO LL i sr 0.2 5 l w0.0 _ - " —-- - O 0.30 Figure II-2 Integrated reflectance for various solutions as function of angle of incidence for isotropic scattering, R vs. 40 parameter of T1' "~-= 0.90 70o

0.70 =o = 0.90 F = 0.9330, P = 0.4821 A Exact 0.60 Six- flux -— Two-f lux O 0.50 Ui LiC 0 LI upr of,} o -= 0.30 Li 0.20 0.10 xr = 0.25 Ti: 0.25 0 0 0.20 0.40 0.60 0.80 1.00 /Lo Figure II-3 Integrated reflectance for various solutions as function of angle of incidence for anisotropic scattering, R. vs. o parameter of T1, (oo= 0.90, F = 0.9330, P = 0.4821 71

WT(a'O = d'o0~6,o = A -0 0 =Ti'O jo cGai;mmua "OTI SA a'9UaaG 3 Tdoa;OSTLU JOJ GDUGpTOtIT JO G{DUa JO uoTpounf sj suo$~nTos snofTJ UoTssTUs1U~ paGs~a u$ snoOaA T-II aGn.Ta 0r/ 00'1 080 09'0 0v0 0O'0 xr;-oMnl: -O.-L- /Q xnl1-x!S 43DX3' 010 ~~I~ 01' o g-o:~7 2~~~0, I /I / 0 / i /I m / r O t' I ~m — I -~~ I OG'O.~ I 3z / U=(rO) I L09'0 7~~~~~~ /3 yI 5) 3/~~~~ ~~~~~~~, ~/~~ 06'0 =o / / -* —---- /08'0.. --- 06'0

0.70 0.60 z / / /' 5 Six-flux soln. cn 0.30 Upper curve: 2 0.50 p=0.5504 ~ z/ Middle curve: rrI P~~P p0.4821 o Q4Q _______ jLower curve: 40 P 0.3889 0.0 / / Six-flux soln. 0.30,(~;(Upper curve: p=0.3889 /- V4 Middle curve: p=0.4821'< Lower curve: p=0.5504 0.20 - -,/ F=0.750, wo=0.90,0rlI.O I Exact 1 0~o P:0.5504 0.10 / p=0.4821 ~/ ~o p=0.3889 ~/ ~ ---— Six-flux soln. I.. —-— Two-flux soln. Any p 0 0.20 0.40 060 0.80 1.0 /Lo Figure II-5 Total integrated transmission, various solutions as function of angle of incidence for anisotropic scattering, R vs. ~o, parameter of P, CA)= 0.90, T1 = 1.0, F - 0.750 73

obviously gives identical values for Rayleigh and isotropic scattering. Representative values are plotted in Figures II-1 and II-2. From the tabulated and plotted values it can be observed that the six-flux solution is generally better than the two-flux solution which is in turn generally better than the Richard's solution. The simple diffusion solution which is not given is poorer than Richard's solution. Exceptions are for normal incidence where the six-flux solution falls off and for infinite thickness for which the two-flux reflectance is invariant with angle of incidence and inferior to Richard's solution. All of the approximations improve as the albedo, co., increases and as the optical thickness, r1, decreases. This behavior can be attributed to the decreased contribution of multiple scattering to the reflectance. Table B-2 gives transmissions for the same conditions and models as in Table B-l, except for the omission of the infinitely thick medium. Again the values for Rayleigh scattering are included for the six-flux model, are identical for the two-flux model and are essentially equal for the exact method. These values are not illustrated graphically but study of the tabulated values indicates clearly that 1) the six-flux values are better than Richard's values which are this time better than the two-flux values. The two flux values are two low for small io but improve for o~ 0.5. Richard's values and the six-flux values are remarkably equal for normal incidence. The greater the albedo, cO., the 74

better are the approximate values. No clear cut trend in accuracy is apparent with thickness. The predictions of transmission by Richard's and the six-flux model are better than the predictions of reflectance but the opposite is the case with the two-flux model. Table B-3 and B-4 give reflectances and transmissions respectively for anisotropic scattering with F = 0.933 and P = 0.4821 as computed from the six-flux and two-flux models. The exact values are also reproduced for comparison. Representative values are plotted in Figures II-3 and II-4. It can be observed that the six-flux model is more accurate than the two-flux model for both reflection and transmission. The six-flux model is most reliable near to = 0.5 and for large 0( and appears to give consistently low reflectances for normal incidence. Table B-5 gives reflectances and transmissions for anisotropic scattering and with F = 0.750 and P = 0.5821, 0.4821 and 0.3889 as computed exactly and from the six-flux and two-flux models. The computed transmissions are illustrated in Figure II-5. The six-flux model yields a dependence or P which is in accordance with the exact values. The twoflux model of course indicates no dependence on P.

B. Interpolation Factors Examination of the approximate solutions and the trend of the exact values provided the basis for the development of several simple relationships for approximate interpolation and extrapolation of the computed values. The integrated reflectance and the total transmission are observed to be strong functions of 10o' o, F and T1 but only weakly dependent on P. Interpolation equations for to and T1 were developed for some conditions as follows. For isotropic scattering by a half-infinite medium simple diffusion theory suggests the interpolation formula s R. II-12 where R and R1 are the integrated reflectance for grazing (40 = 0) and normal (I = 1) incidence, respectively. This equation is found to be satisfactory for interpolation of R with respect to 40 even though the absolute values of R indicated by diffusion theory are seriously in error. The success of this interpolation is indicated by the top curve (for tl = oo) in Figures II-6 and II-7, and in Table II-1. Richard's model does not give accurate values of R for isotropic scattering by finite slabs but the indicated dependence of R on 1 is 76

0.70 w0, =0.90 Isotropic scattering o Exact 0.60 o A. Interpolation. -oo a T Extrapolation 0.50 () 0. 1.0. lI -J 0.40 LLu 0.50 bLJ <0.30 Lu I Z, 10.25 0.20 0.10......_ O, I, 0 0.20 0.40 0,60 0.80 1.00 Figure II-6 Integrated reflectance, exact, interpolated and extrapolated as function of angle of incidence for isotropic scattering, R vs. [o, parameter of T1, A)= 0.90 77

0.40 Isotropic scattering o Exact 0.35 4 0.Lo Interpolation A rI Extrapolation 0.30. 0.25 0 _ _ oao of an -J LL LdJ z' =0.50 0.10 Z1 =00 0.05 0 0 0.20 040 0.60 0.80 1.00 Figure II-7 Integrated reflectance, exact, interpolated and extrapolated as function of angle of incidence for isotropic scattering, R vs. [o, parameters of (o and T1 78

Table II-1 Interpolation and Extrapolation of R for Isotropic Scattering with Respect to.o and -1 cO Q) 0.90: e, = ~o;, =.o = O.5, = 0.25O Exac%) Interp. *t* * 0 Exact Interp. Exact Extrap. Exact Extrap. Exact Extrap. 0.00000 0.684 o.684 0.624 0.625 0.584 0.586 0.571 0.556 0.04691 0.664 o.664 0.598 0.606 0.548 0.561 0.498 0.516 0.23077 0.582 0.595 0.504 0.511 0.417 0.433 0.300 0.316 0.50000 0.508 0.517 0.394 0.401 0.279 0.293 0.173 0.185 0.76923 0.452 0.456 0.316 0.322 0.209 0.215 0.125 0.129 0.95309 0.423 0.423 0.279 0.284 0.179 0.184 0.105 0.105 1.00000 0.415 0.415 0.271 0.274 0.170 0.176 0.100 0.100 o4 = o.60: No = 0.30:, =- oo O,= 0.50 T, = oo =r, O.5o C10 Exact Interp. Exact Extrap. Exact Extrap. Exact Extrap. o.ooooo 0.368 0.368 o.345 0.352 o.163 0.163 0.156 o.161 0.04691 0.334 o.346 0.318 0.321 0.144 0.150 0.142 0.141 0.23077 0.267 0.279 0.232 0.236 0.108 0.114 0.100 0.098 0.50000 0.212 0.218 0.154 0.152 0.082 0.085 0.066 0.061 0.76923 0.177 0.179 0.114 0.109 0.066 0.067 0.048 0.045 0.95309 0.159 0.159 0.097 0.091 0.059 0.059 0.041 0.036 1.00000 0.155 0.155 0.093 0.087 0.057 0.057 o.o40 0.035 * Extrapolated by Lagrange's formula for 4 = 0 and o = 1.0 79

excellent suggesting the equation Re = Err, Rpa, II-13 /1 ro,, T, Ose n ~R;chorJ; for the extrapolation of exact values of R for 1 = oo and any,o and o, to finite values of T1. The success of this formula is indicated in Figures II-6 and II-7 and in Table II-1 for the finite values of'1. The ratios of R for Richard's solution were taken from Tables B-l and B-2. The success of the interpolation equations indicated that exact calculations for isotropic scattering could be limited to normal and grazing incidence on a half-infinite medium. For anisotropic scattering the expression R = XR, + ( I -X)Ro II-1 where x is a function of co, ~1 and po but not of F, can be used to interpolate for o0. Values of X can therefore be determined from computed values of R for isotropic scattering. A plot of X as a function of ~o for several values of T1 and'0. are shown in Figures II-8 and II-9. These figures were prepared from the isotropic values. (The six-flux solution erroneously indicates independence of X from ~1 and 40 as well as from F and P.) Values of R for intermediate values of o were computed from the exact values of Ro and R1 and the values of X in Figures II-8 and II-9 using Equation II-14 and are given in Table II-2 and in Figures II-10 and II-11 (triangles). The agreement with the exact values is seen to be reasonably good. 80

1.00 z o:= 0.90 w -0.80.. CP C) U E040 00.60~ 0 IL C/) _0.40 LL 00.20 pIL_ _' _I, 0 0 0.20 0.40 0.60 0.80 1.00 ~.o Figure II-8 Factor x in equation II-14 calculated from isotropic R and Ri, as function of angle of incidence, x vs. o, parameter of (OZ

o 1.00 z nL rI oo w IH~~~~~~~~~~~~~~~~~~~~~~~~ o 0.80 U)\O 0 0 n, ~ rr a H 0.60 I iI LL U) U~~k......i.. 0 0.20 0.40 0.60 0.80 1.00 Uo Figure 11-9 Factor x in equation II-14 calculated from isotropic R and R1 as function of angle of incidence, x vs. ~o, x~~~~ 0~~~~~ parameters of T 82

Table II-2 Interpolation of R with Respect to, o for Anisotropic Scattering F =0o.9330, p=o.4821, oo= 0.90:',=oo = 1.0'0. =?0.257 Exact+ Interp. Exact* Interp. Exactt Interp. Exact* Interp. 0.00000 0.628 0.628 0.567 0.567 0.547 0.547 0.552 0.552 0.04691 0.595 0.599 0.530 0.530 0.500 0.502 0.466 o.468 0.23077 0.489 0.478 0.396 0.397 0.333 0.338 0.238 0.241 0.50000 0.377 0.370 0.239 0.241 0.163 O.166 0.097 o.o96 0.76923 0.292 0.299 0.131 0.131 0.078 0.079 0.043 0.041 0.95309 0.244 0.245 0.079 0.078 0.041 0.041 0.018 0.019 1.00000 0.233 0.233 0.067 0o.o067 0.030 0.030 0.013 0.013 F= 0.9330, P = 0.4821, oA)= O. 60: W0= 0.30: _,= o To 0F0', = o00'= 0. o 4lo Exact* Interp. Exact* Interp. Exact' Interp. Exact Interp. 0.00000 0.336 0.336 0.330 0.330 0.157 0.157 0.151 0.151 0.04691 0.303 0.289 0.293 0.296 0.144 0.139 0.131 0.134 0.23077 0.206 0.197 0.180 0.189 0.108 0.105 0.075 0.082 0.50000 0.124 0.121 0.085 0.092 0.082 0.081 0.034 0.040 0.76923 0.073 0.073 0.040 0.042 o.o66 o.o66 0.015 0.018 0.95309 0.048 o.o49 0.021 0.021 0.059 0.059 0.008 0.009 1.00000 0.043 0.043 o.o16 0.016 0.058 0.058 0.007 0.007 F=.75'00, P= o.4821, w0o 0.90; I= 1.0 410 Exact* Interp. 0.00000 0.598 0.598 0.04691 0.564 0.566 0.23077 0.445 0.451 0.50000 0.310 0.317 0.76923 0.218 0.221 0.95309 0.175 0.176 1.00000 0.166 0.166 * Extrapolated by Lagrange's formula for 0o = 0 and Ko = 1.0 83

0.70 wt0O.90 0 F=0.9330, P = 0.4821 Exact 0.60 l o F=0.7500,P=0.4821 Exact A jo Interpolation 0.50 0.40 L O0 I _o 0.20. 0..1.0 0.20 0 _ 0 0.20 0.40 0.60 0.80 1.00 I-o Figure II-10 Integrated reflectance, exact and interpolated as function of angle of incidence for anisotropic scattering, R vs. 4o, parameter of h1' o~= 0.90, F = 0.9330, P = 0.4821 84

0.40 F= 0.9330,P=0.4821 o Exact A.l0-Interpolation 0.35 0.30 z 0.25 U_ LLJ Li U)A,-," 0,20 p u~iI ~ ~~ 8 0Q O~~=0*60 H 0.15 0.10 " i=0.50 0.05 0 0 0.20 0.40 0.60 0.80 1.00 /1o Figure II-11 Integrated reflectance, exact and interpolated as function of angle of incidence for anisotropic scattering, R vs. I1o' parameters of mo 85

The same type of interpolation, i.e., T=y ~r * (r~y)T, II-1F T = I To + ( i -,) To II-15 was also successful for the total transmission as indicated in Table II-3. 86

Table II-3 Interpolation of T with Respect to [o for Anisotropic Scattering F = 0.9330, P = 0.4821, CO= 0.90:', = 1.0 = 50, = 0.25 p. Exact* Interp. Exact* Interp. Exact* Interp. 0.00800 0.761 0.761 0.761 go4 o.904 0.991 0.991 0.04691 0.750 0.754 0.885 0.886 0.979 0.983 0.23077 0.692 o.698 0.833 0.841 0.933 0.951 0.50000 0.553 0.551 0.723 0.733 0.853 0.881 0.76923 0.372 0.365 0.507 0.505 0.669 0.689 0.95309 0.259 0.256 0.331 0.328 0.396 0.402 1.00000 0.236 0.236 0.291 0.291 0.297 0.297 F= 0.9330, P=o.4521, F=-o.75oo, P=o,-481 u,,= 0.60: J0 o= 0. o 4o0 0.90: r, = o.50', = 0.50 VI, = 1.0 Co0 Exact* Interp. Exact* Interp. Exact* Interp. 0.00000 0.783 0.783 0.692 0.692 0.652 0.652 0.04691 0.761 0.761 o.668 0.665 0.645 0.646 0.23077 0.701 0.701 0.601 0.600 0.598 0.598 0.50000 0.565 0.561 0.452 0.433 o.476 0.473 0.76923 0.319 0.309 0.198 0.189 0.318 0.313 0.95309 0.165 0.161 0.063 0.059 0.222 0.221 1.00000 0.143 0.143 0.052 0.052 0.203 0.203 * Extrapolated by Lagrange ts formula for 40 = 0 and po = 1.0 87

I

PART III - DEVELOPMENT OF IMPROVED APPROXIMATE MODELS A. Introduction It is well known in classical physics that multiple scattering processes can be described approximately by means of the classical diffusion equation. The major difficulty in the application of classical diffusion equation to the problems of multiple scattering however, is in the specification of the correct boundary condition. Various authors have investigated the use of the simple diffusion equation in multiple scattering problems, and formulated modified (13,14,15) versions of the diffusion equation. In this investigation it is shown that the various modified diffusion equations can be derived in a unified manner from the integral equation of transport, and therefore can be not only refined indefinitely, but also, in principle, can have appropriate boundary conditions derived corresponding to each step of refinement. On this basis a new modified diffusion equation is proposed and solutions compared numerically with known exact solutions. Better agreement is found than with previous models. The boundary conditions are obtained from purely mathematical reasoning and a good physical interpretation is yet to be found. To simplify the mathematical formulation, isotropic scattering has been considered as an example, but the extension to anisotropic scattering is trivial, and will be carried out in the continuation of this work. 89

B. Approximate Differential Equations For isotropic scattering, the integral equation of transport for the photon density ft~) is P(Y) =,_r _v Ae Cat/') d-' III-1 p) 5 -t —-~I3 7,/ 3 V 4niQ of dCe rs; Or) where S(O) is the source term, WD is the albedo of single scattering, and y represents the space coordinates measured in mean free paths. Equation (1) is a linear integral equation of the Fredholm type in three dimensions. SC co fffK(Y ) f c d) III-2 For some special forms of the kernal K(rr ) this equation can be reduced to a simple differential equations. In this section, it is shown that the various forms of diffusion equation can be obtained by such reduction after approximating the kernal K(Y,r') = m e by comparatively simple forms. Consider the class of K('r,i ) whose Fourier transform can be represented by ratio of two polynomials*: _(___ a0 a 21k. D tKk) o ( blt b A k... * The same derivation applied to a two dimensional case was given in Ref. (16) 90

so that, K(V r ) (2r)3 I J D (*2) e d -k~x l do III-3 Then, since Equation (3) may be reduced to: O (-r -) t D) ) III-5 Substituting Equation (5) into (2) yields?() S- (r) -+'M /(-v2) (t)' u -(/ ye III-6 Now, operating on both sides of (6) by D(-Vt), (- V 2)>y) D ( )+ w ( -V2)'( s4,( c. d/ e6 III-7 ofvdiSU,'n It is recognized that f7r)2 -3 JJJ k J the three dimensional 6 -function. Since Equation (7) may be reduced to: Lg(-v)- w-/l(-v&1p(r) = PF- V)) 1(r III-8 Equation (8) cannot be used directly for Equation (1), because the kernal in Equation (1) cannot be exactly expressed in the form given by Equation (3). In fact, it is easy to see that 7Y7j -O 7/2/Xt 9 1 91

Thus a whole class of approximate solutions of Equation (1) can be obtained from approximations of the function Ad 1- (_1 He1' 4 III-9 in terms of ratio of two polynomials. It is interesting to note that although the classical diffusion (14) (13) equation, the modified diffusion equation proposed by Richards (15) and of the modified diffusion equation proposed by Grosjean, were obtained from various physical and/or mathematical arguments, they can also be obtained in a unified manner by different approximations of Equation (9). These are respectively;: F(~) ~ it - e + 11: i1. III-o10 so that Equation (8) takes the form 75(t - 3(v-wo)e(r) - - (3- V7>) 3 () III-11 )F'1)1+ III-12 so that Equation (8) takes the form: CCV r (1_w3 XD)5trJ 1_3 Ad 5(r ) III-13 )(j) ^ > Z. k2 III-14 92

so that Equation (8) takes the form: Cr ) X (3- 2 S (3) III-15 The left hand sides of Equations (11), (13) and (15) are the equivalent terms of the classical diffusion equation, the modified diffusion equation proposed by Richards and the modified diffusion equation proposed by Grosjean respectively. It is mathematically difficult to estimate the errors involved in these approximations, perhaps the only valid criterion is to compare the results obtained in various simple problems with exact solutions. Intuitively, however, it seems to be reasonable to require any approximation of F(k) to be as close as possible to the exact function at least in the limiting case of k = O, and k-4co. Inspection of Equation (9), (10), (12) and (14) shows that as k-O all three approximations are exact to the coefficient k2 of the expansion of F(k). However, as k-oo, and only the simple diffusion equation has the correct limit. Based upon the above observations it is clear that a class of differential equations based on diffusion type of approximation can be systematically obtained as follows: (i) First approximation: FrkV) a lii jL 3e 93

For which, the approximated kernal is, _1-, A et~/ ~A~~ ) III-16 and the resulting differential Equation is: v:(Y) 3(, -o)FP_(L - (V-3)S(7) (ii) Second approximation For best approximation at small values of k, 8 a = 21 5 b= 7 and 4 c = 105 so that / -AZ i 72 * III-17 corresponding to this approximation. I - III-18 where 31 = 1.2344 2 = 4.1505 are the positive roots of the equation cpa- 4- 4 I = O 94

while n,i, - 10/3 and The resulting differential equation is: In principle, c + can be extended indefinitely. In principle, the approximations can be extended indefinitely. But numerical calculations given in the subsequent sections indicate that with the correct form for the boundary conditions, the first approximation yields very good results for the photon density, and the second approximation is almost indistinguishable from the corresponding exact solution.

C. Point Source in Infinite Medium For the problem of multiple scattering involving sources in an infinite medium, reduction of integral Equation (1) into differential Equations (11) and (19) enables one to find the approximate solution by use of Green's functions. It is well known that the solution of the equation v (jr 9 ~- K G J- -1t) III-20 subject to the condition that ~i r is zero as r->oc is,....L _ # e I cx (I, y") -- 4 rr 1?~- A'III-21 This Green's function can be used to solve the equation of the lst approximation Let K= 3 (6 — ) multiply Equation (20) by -(3-;7)S(B3 and integrate. It is easy to see that the solution of (11) is * means the appropriate differentiation with respect to the primed coordinate variables 96

This yields* i ) C= s () -+(3- t<): J/f c ) s() dcU9 III-22 or more explicitly:?(; s(? t- TJ J? s 1r III-23 Similarly, the solution of the equation % ~oo-bv _ +.. ~; =W _o ), -== III-24 subject to the condition that;(r;'-O as r —>oo is, (,rJ >1-t 3'L c r I' L 31'2 III-25 where cal, ac2 are the positive roots of the equation** ~( Z (1 C III-26 * Equation (22) is obtained by using the second Green's identity, so that jJf/i CT (P)' P jd' -) I %llc'fJrd1 But in view of Equation (19) f SfP) - 9K'/if/sf r 7 ) - Cr?)'" ** For the present problem, wo, alc, ca2 exist, and are given in Table III-t. 97

Using the Green function Y(r,I') it is possible to write down formally the solution of the equation in an infinite medium. [q t- =; _2 bs () III-27 The result is: (r) s(~9 Mre qlr )le 4 - 2 III-28 where' (-iz 111- III-29 and I-2a m ad, 5 e: C III-30 Table III-1 The explicit values of al1, m2, m 2 as a function of 0 1.2344 4.1505 1.7013 9.2987 0.1 1.2039 4.0374 0.7916 9.2084 0.2 1.1683 3.9224 0.8987 9.1013 0.3 1.1264 3.8054 1.0264 8.9736 0.4 1.0765 3.6866 1.1793 8.8207 0.5 1.0159 3.5662 1.3632 8.6368 0.6 O.9408 3.4446 1.5847 8.4153 0.7 0.8447 3.3221 1.8517 8.1483 o.8 0.7161 3.1996 2.1721 7.8280 o.g9 0.5264 3.0778 2.5532 7.4468 1.0 0 2.9580 3.0000 7.0000 98

Equation (23) and (28) can be directly applied to the simple problem of an isotropic point source of unit strength at the origin in an infinite dispersion. The source term may be taken as the unscattered radiation, -r S(r = ~) III-31 Substituting Equation (31) into Equation (23) and (28) yields, respectively, for the first approximation: C[f) —-. +- U)- -- C 333 F-, J3Ci, 7 J TjjIII-32 t X: {) 1(1YC- rn and the second approximation, where E,@) is the exponential integral defined by ~E,~(X) ~ = ~ 0 ~ z~ta~ ~III-34 Some numerical results calculated from Equation (32) and Equation (33) are compared with the known exact solution in Figure III-1 for u )0= 0.3. Results of both the first and the second approximations are very close to the exact solution whereas the solution of the simple diffusion differs significantly. 99

1.00 0.80 - 0.3 ~o 0.3 0.60 - Exact soln. --- Simple diffusion \ zI ~a Ist Approx. 0.40 \ o 2 nd Approx. 0.20 ~ I0.10 0.08 0.06 0.04 _.... 0.020 _ _ _ _ 0 1.0 2.0 3.0 4.0 5.0 Figure III-1 4tt>'(r) as a function of r for an infinite medium with point source 100

D. Boundary Conditions The excellent numerical results obtained from the approximate solutions in the case of an infinite medium suggests the possibility of applying the same approximation to a finite region of dispersion. It is well known that the classical diffusion approximation yields very poor results near the boundaries of a dispersion. This inaccuracy probably results from poor representation of the actual boundary conditions. In this section, instead of using physical arguments, a rigorous mathematical derivation for the boundary conditions appropriate to each order of approximation is derived. It will then be shown that the resulting boundary conditions yield results that are very close to the exact solutions. Mathematically, for a finite region, the integral equation to be solved is for) = y) dir III-35 where V is the volume of the dispersion. In one dimensional cases such equations are known as Weiner -Hopf equations and has various application in electrical engineering problems(l7)(18) For the present problem, the kernal is of the exponential type, and may be represented by III -36

The integral Equation (2) may be transformed to the differential Equation Note that if c = a -o, and b = 1/3, Equation (36 ) and (19) Note that if c = a = o, and b = 1/3, Equation (36) and (19) represents the first approximation, so that the discussion that follows is applicable to both the first and second approximation. The complete solution of Equation (19) is in general made up of a particular solution depending on 5(Y) and a complementary function which is the general solution of the homogeneous Equation, i.e. Equation (19) with S(-= 0. The complimentary function would of course involve arbitrary constants, which are to be determined by the boundary conditions. Since Equation (19) is actually derived from the integral Equation (36), the boundary conditions can be derived directly by mathematical manipulation without resorting to physical arguments. Let V be the volume of the dispersion, r the boundary of V (assume convex) and 1m be the outward normal to f from V, then the value of p(v) at the boundary can be directly obtained by use of Green's identity. Values of p(v ) at the boundary 7 must satisfy the condition 102

=- 0fjOu~Thv( vZYb'-}~III-37 for all o inside V. The proof of Equation (37) is straight forward, and is outlined in the following steps: (i) Multiply Equation (19) by K(ir-'l) and integrate over v yielding: JJfK(-',) E l- b cVt 7(?) a LJ~ - j Y/(lr-' V < 7$3 5(r ) d " - rA)0 JtJ K~lt-tl) 01-7 Jcr) d C7 o III-38 (ii) Using Green's identity, Equation (38) is reduced to: 103

If lif~~ is 4 vx - ~ - E /f S gt) - Ad'73 vv, ~,~. (,, \;,r -.'I;d — bJ(r Elt- Ovbr3 -p: C:T/ (I- z i. -+. c ff -' V2"&) <r) V vj j rA i J + C'- ) v bury - o(t T g( r- L/slC t 4 + f L: -t') v god - J.c-) - cff fTs( ('o v v,/c(/e-e,/.J. ~ 4 c ffE:c'~OGy) v 2, (I - sOl-' )?: i F ther(: exnt0vpr:o ( given in Equai 3 III-39 (iii) From the expression of K (Jr-r/) given in Equation (36), it can be shown that the left hand side of Equation (39) is identically zero. Collecting the terms given in the right hand side of Equation (39) produces the boundary condition given by Equation (37). Equation (39) is the general boundary condition that can be derived. It can, in principle, be applied to any geometry and source distribution. However, in applying it to simple problems some simplified version will be considered. Before considering the simplified version of the boundary conditions, it is to be noted that Equation (35) is a linear integral 104

equation of Fredholm type, and therefore can be solved in terms of orthogonal expansion without resorting to boundary conditions In principle, for the particular type of K(i-r ~' it is possible to find a complete set of eigenfunctions (r) corresponding to eigenvalues X, satisfying the monogeneous equation t;(r>'~ if~ (( <Erl,,') rn 111ii-40 V If such a set can be found, then, it can be shown that the orthogonal to each other in the sense that This immediately leads to the solution of Equation (35) in a series form and the result is: yCrc Z vi III-41 where z AX h~~~- bJ o,1 III-42 and Bi is related to the source terms by t- -- -Y..... III-43 The approximated kernal given by Equation (36) enables one to find and ( in a systematic way. Since r(rJ must satisfy the homogeneous differential equation rcv tCb) (+ C —n) r III-44 the forms of Cr-) corresponding to each >%, are known in terms of arbitrary constants and ). Substituting the expression for 0,,(' ) 105

thus obtained into Equation (40) yields the characteristic equation for X;, and the constants for 4(t. In this one dimensional case, such a procedure has frequently been used e.g. (18). The extension to three dimensional case is straight forward. In the present work, the primary purpose is to test the accuracy of the approximation, so that only simple problems such as parallel dispersion of finite and infinite thickness are considered. For such one dimensional cases, both the integral Equation and the boundary condition can be greatly simplified. The boundary condition used in the present work therefore combines the idea of the above two methods to minimize the numerical calculations. In the next section, such simplified boundary conditions are illustrated. 106

E. Solution of One-Dimensional Problem 1. Second Approximation a) Medium bounded by parallel planes To test the accuracy of the proposed approximation consider the problem of a parallel plane flux of unit intensity obliquely incident on a parallel dispersion bounded by two planes Z = 0 and Z = L. If the direction of the incidence is inclined at an angle 8 = Cyl O from the Z-axis, then C(r is a function of Z only so that Equation (19) is reduced to: Lda-4 tz ~ (( I -t(I-2 7 {' s(d where S(Z) is obviously of the form S(z - e Y~ The general solution of Equation (19) can be written in the form: Z)-A, + 4,' - + (1+ D)) III-45 where;Oo~(.tct(- III-46 The constants Al, A2, B1, B2, can be obtained from the original integral equation. In the present one dimensional case, this integral equation takes the form: CZe No+ o LA t1 - )dr' III-47 107

where K(-Vl) is given by Equation (18) as KY4l) =t i 71- 2 me By carrying out the integrations that are independent of Z, Equation (47) is reduced to the form: - -(Z dh u III-48 Now, substituting Equation (45) into Equation (48), and collecting terms involving the same exponential forms of Z, the following boundary conditions for A and B are obtained (i) For W,4l +t — _ -- 0 dA,j + t2 il e ha 4.2 I-+-D The set of equations can be solved for the four constants, _+ -- 0 - I10~8~ O -+ - - -+ Z- --'4 + 2The set of equations can be solved for the four constants, lO8

The above equations are valid for. (ii) For u0o = 1, (See table 1) d&O 0= i0. 80 So that the solution of Equation (19) is' — A- A, 5 42 -t, t- +t C e (i ID)1-49 The "boundary conditions" for determining Al, A2, B1, B2 are slightly modified to: - _ A2. - 5- +0 Ai Az pi i L 1 4X pe -h~s tL-, 0.+,? 4, + d~,..i''+ r ~C = 0 i_ [I~ ~ d -0- - Al. CtPIt1,0 9t 3>e 5' _8. (I-P)~ 5 109

b) Semi-infinite medium For a semi-infinite medium, C,-co the results can be simplified considerably. Following the same procedure as in the finite case, (i) For 0o +# 1 pec- A= le +, -_D (IDS A~ O III-50 where D is given by Equation (46) and Al and B1 are determined from: AI e, l\(ii) For ct o 1. 3 _, (1ID) 4 ~ III-51 where Al and B1 are determined from: 3 __ _,I jo A, - +i 0 2. First Approximation The formula above were derived for the second approximation. The corresponding formulas for the first approximation can simply be obtained by letting a = c = o and b = 1/3. The explicit results for the one dimensional problem, using the first approximation are given below: 110

(i) Integral Equation 7Zi 111 e + < eUz~d III-48' (ii) Differential Equation,? iEa -(, ) t9(Z) - O I ) ef III-19' (iii) General solution of the differential equation D(Z)- A,1-dl'-+ Azc'J Z (D III-45' where an D (-)- 830l III-46' and a, = C3(1- WOt (iv) Semi-infinite dispersion, Wo0 I p(z) 41e16 i (+)e III-50' where A - 111 d)III-52 (v) Semi-infinite dispersion, a = | e(z) = A I -+ (+i ) e / III-51' A III-53 (vi) Parallel plane dispersion o# $ i gIi= Az4 Aze ifA- p)e III-45' where A1 and A2 are determined from, 111

(vii) Parallel plane dispersion, 1i4= 1 ez) A- +- A2z - (- zJ III-49' where A1 and A2 are determined by 1I,__A1 2-tA, 3'LlS 112

F. Calculation of Albedo Available results on the exact solutions of the parallel plane problem are generally expressed in terms of the reflectance of a half space, and the transmittance T or reflectance R of a parallel plane dispersion. One may calculate these parameters by direct integration of the expression for p(Z). The resulting expressions are given below: 1. First Approximation (i) Semi-infinite medium, III-50' fPzj3z- AG~l'-~ (:fl4GJ7Ao ~~~~Z) C 4 -ditt00 \D ~-~ R~i, ~~~~~~~ C ~ JZ_'~IA oJ 11-541 (ii) Semi-infinite medium, J = f(z4Y A+ (41) III-51' R =~ A!-f- P) _I 111, r III-55' 4/A0 ~~2. (iii) Parallel plane dispersion, u%* I e(z)== A&2e - Ae III4-45' Then R -A,-+A,ndl + K-L III-561 and'U 113

where and III-58 + e -Cs -,'', ( — ^ (iv) Parallel plane dispersion, - (z) A,-X Az z C+9 ] III-49' WV III-60' IT- A~4 wp 1 3(n]) + 42' L3u'4m -t + (\tDfii III-61' IZI-54 2. Second Approximation (i) Semi-infinite medium 04o= /-G A -+ EBie a- b-+1 III50 e Xo8A, -4 1- lAtlt iT+D +) v F III-541 (ii) Semi-infinite medium LC0= *() -- AL 1- III-51 I- B-'i i-+.....

(iii) Parallel plane medium )z>= A, —A'': -""'. c1:C III-45 -R -,-1 --, - C,. Bz + (I L)s 111-56 T>= A7T-I AT- T.E - P2T7 f1)T0 III-57 where R,. and Toa are defined as Equations (58) and (59) respectfully. (iv) Parallel plane dispersion h-t (HAi)TJ (z)f~~~- - Az - i 4,, 3T __?LZ3'..'.,

G. Numerical Results Some typical numerical results calculated by using the appropriate approximate formulas derived in the previous section are plotted in Figures III 2-8. In Figure III-2, the integrated reflectance R for a semiinfinite medium is plotted versus So for cn= 0.9 and 0.5. The results of the second approximation are indistinguishable from the exact values. In Figures III 3-5 the integrated reflectance R is plotted versus -W for several values of po. The numerical values were calculated from the first and second approximations. Again, the second approximations are indistinguishable from the exact values and the first approximation are acceptably accurate. In Figure III-6, the reflectance R for a finite medium is plotted against 0o for several it. Again the results of the second approximation are very good. The deviation of the first approximation from the exact solution becomes bigger as So increases. In Figures III-7 and 8, the diffusive transmission and the integrated transmission are plotted against So for several 0,. Only second approximations are compared with the exact solutions. From all these figures, it is obvious that the results of the second approximation are practically indistinguishable from the exact solution for all conditions, and the simple diffusion equation with correct boundary conditions (i.e. the first approximation) is an acceptable approximation. 116

0.70 T1 = o0 \ \ ----- Exact soln. --— Ist Approx. 0.60 0.60 \ % 0 2nd Approx. W0.50 z <) w) =0.90 Lu Lw 0.40 Lu o 0.30 z 0.20 O. I0 o 0.50 0 0 0.20 0.40 0.60 0.80 1.00 /Lo Figure III-2 Integrated reflectance as a function of incident angle, R vs ~o, parameters of ~o, and ~1 =0 117

1.00 --- 0.90 TI = a Tt =cO I bo= 0 Exact soln. 0.80. Ist. Approx. o 2nd Approx. -- Simple diffusion 0,70 Ix LU 0.60 / LU 0.50 Q / LUJ 0.40:i I_! LLI z/ 0.30 / r/ 0.20 0.0 / 0.01 / / 0 0.20 0.40 0.60 0.80 1.00 Co Figure III-3 Integrated reflectance as a function of albedo for single scattering, R vs (A) O~, = 0, and T1 = 0o 118

1.00 Lo = 0.50 0.90 Exact soln. 0 2nd Approx. Ist Approx. 0.80 - Simple diffusion a: 0.70 z 0.60 0.50 t9 0.40 -1 0.30..1. 0.20 2/.,o/~/ 0 0.20 0.40 0.60 0.80 1.00 Wo Figure TI1-4 Integrated reflectance as a function of albedo for single scattering, R vs. &., po = 0.5, and T1 = 0 119

1.00 0.90 [L =1.O Exact soln. 0 2nd Approx. 0.80 -— I st Approx -— Simple diffusion 0.70 z / < I LU LLJ I 0 0.50 LUJ 0 <E I 0.40 O z I // 0.30 // 0.20 0.10 0 0.20 0.40 0.60 0.80 1.00 Figure 111-5 Integrated reflectance as a function of albedo for single scattering, R vs.& )o0 - 1.0, anT Tc - xo 120

o00 = 1 putv' OCQ JO saaGui'0d` si 9A'9aTU Lueapiou! jo uoL4ounfJ e Se auuaeDGaJG pa;UaDG@uI 9-III Gan2FTa OT/ 00 1 08'0 09'0 0'0 OZO O 0 I l'O: Om01'0 O' = O\ O - 06i 0 0' I OM = Ovl \ O O rn C-) o'o- o.1 09'0'xojddV is I v'UlOS IODX3 o ~xoiddv pul - 0'1 ='l OZ'O

0.380 0.360 0.340 - o-1.0 0.320 o 0.300 0.280 Ho~ 0.2 60 re I/ Io i i a a f0.90 L 0.24o0 U- 0. 00..20 0 > 0.200 0.180 0.20 0.40 0.60 0.70 1.00 ILo Figure TTT-7 Integrated diffusive transmission as a function of

0.80 0.70 -- = 1.0 - Second Approximation o Exact Solution o 0.60 H 0.50 0 w 0.40 zn 0 0.30 0.20 0.10 0 0.20 0,40 0.60 0.80 1.00

H. Conclusions From the excellent numerical check with the exact solutions of one dimensional problems with isotropic scattering, it is reasonable to believe that the diffusion-type approximation, if used with correct boundary conditions, can yield very good results which are adequate for all practical purposes. In the present work a systematic approach to the diffusion-type approximation has been developed. A method of determining the correct boundary conditions appropriate to each order of approximation has also been developed. Since the exact solutions of the multiple scattering problem are generally difficult or impossible to obtain, the need for a comparatively simple, yet reasonably accurate approximate method of solution is obvious. It therefore seems worthwhile, to exploit further the possibilities of using the diffusion-type approximation. This extended investigation should proceed along the following two directions: (i) use of the suggested method to find approximate solutions for geometries where exact solutions do not exist, such as a point source in front of a parallel dispersion, (ii) extension of the method to anisotropic scattering, for which some exact solutions are now available for comparison. 124

SUMMARY AND CONCLUSIONS 1. A satisfactory method was developed for numerical integration of the transport equation for parallel plane radiation obliquely incident on a parallel-plane dispersion. The important advance is in the allowance for strongly anisotropic scattering. Previous results were limited to isotropic or Rayleigh (nearly isotropic) scattering. 2. The solutions are exact except for errors in the numerical procedure itself. These errors arise from numerical roundoff, from termination of the iterative procedure at an arbitrary point and from integration by quadrature. The error due to the use of quadrature arises from the use of a finite number of terms, i.e., from division of the integration into a finite number of intervals. The numerical results agree very well with previous results for special cases (isotropic and Rayleigh scattering, and anisotropic scattering by semi-infinite dispersions)indicating that the net numerical errors are not serious. 3. The numerical and graphical results illustrate the effects of a) the phase function for single scattering b) the optical thickness of the dispersion c) the angle of incidence d) the albedo for single scattering (ratio of scattered to scattered plus absorbed radiation) 125

4. The numerical results are limited to finite absorption since the iterative procedure did not converge satisfactorily for zero absorption. It appears feasible but unnecessary to develop a separate numerical procedure for zero absorption since the results for small finite absorption can readily be extrapolated to zero absorption. 5. The reflectance and transmission for strongly anisotropic scattering was found to differ decisively from that for isotropic and Rayleigh scattering. A similar result was previously found for semi-infinite dispersions. It is therefore concluded that solutions for isotropic scattering do not provide a reliable guide for the reflectance or transmission for dispersions of liquid or solid particles. 6. The transmission and reflectance were found to be dependent primarily on the fraction of the radiation scattered into the backward (or forward) hemisphere. Only slight dependence was found on the other characteristics of the phase function. This suggests that the complicated phase functions associated with large circumference to wavelength ratios can be approximated by simple functions without serious error if the back-scattered fraction is preserved. This result is important from a computational point-of-view since the extent and cost of the computations goes up rapidly with the number of terms used to represent the phase function. 126

7. Many approximate models have been proposed for multiple scattering. The solutions based on these models were compared with the "exact" values for both isotropic and anisotropic scattering. The six-flux model was found to be the best of these models but to be a fair approximation, Richard's modified diffusion model and the two-flux model were poorer in all cases and the simple diffusion model was seriously in error in all cases. 8. Despite their inability to provide good absolute values for the transmission and reflectance the above approximate models were successfully utilized to provide simple expressions for interpolation and extrapolation of the exact values. 9. A new variable-order diffusion-type model was developed which yields very accurate results for isotropic scattering. It is proposed to develop this model further for anisotropic scattering. 10. All accurate solutions for multiple scattering are currently limited to one-dimensional radiant transport. However both the numerical method and the variable order diffusion-type model hold promise for the computation of two-dimensional radiant transport. 127

APPENDIX A Derivation and List of Equations 1. Derivation of Explicit Relations for bm(li, i, 1o) and cM (~l~ 1, 0)' The diffuse reflected intensity is N N + I(o,/ ~() = ~2 (2_ $ o,)Oc mX (-,) Ft (Z,/o) -14 substituting relation I-3 for F1'(Z~,#,;o) and letting 274 = +' __ Al-l gives I. (c- o,) c?,, i o (-I) At [,,, 9' t; c,, -,/,, ~,. I* () 4oy (pM) (L )A A1-2q1(l,A3m and for N = 2 co5s Z 2cosf' - A. 1 IC-,,,~ = 4e /. It: c),,/ (/o) Q,/ (I)]-A, AI(.)/,(I)-I(),( -As 3 [ C, (?.)~'(?) _ 4 (,), c(?)] A1-3 128

Therefore, let A= XL-e- o 1-A1-4 + A L (Kz )+() - (e) 2 A1-5 t( )?Al.) Q, = [,c: ~o: c?)- 9, M mo>9,C)] -A, [p(?o)?} -~,?o),(p)] A1-6 2 4 JOf? A4 f Al9 Q (# = A5 t)+() s)'(A1-7 and the expressions for bm (h l,d,/) become bo- K (Qo -Q ) A1-8a b, = KQ, Al-8b b I Y /Qz Al-8c The diffuse reflected intensity may be written as r b(o + ib cos (f + b. cos Al-9 Proceeding in the same manner to derive the relations for Cm(-[,l, o0), the expression for the diffuse transmitted intensity can be written as _ _ (o(s_) 4 T I-19 Equations I-14 and I-29 are of the same form with the term ~ -eI m (-I) 1 (Z,,/,M) in I-14 replaced by T,, I A, h.) in I-19. The only difference between F; and Gl is that the positions of *If(} and (/M are interchanged and the term r +,o in F1 is replaced by -~o in. 129

Therefore, if 1. All signs are made positive in front of the Al s 2. ~ + 0o is replaced in K by [ - 0o 3. c1(f) and 1A(r) are interchanged The expression for the diffuse transmitted intensity may be written directly. K' = ______ Al-12 ~ 4, —- rA-, (?~)' (~)' (?U), (,)] +A'?(~)',(?-'~,C.) (%)J + A (2?2. - (U)3 Al-11 Q A,[ (V()#,Cf) #Q'ilQ')) + 3 ( A1-12 A2 -6 + (r')92 (r) 3: (t~) If~ (,)] Al-13 6 with the expressions for cm CO= k[.' - Q1 Al-14a C, (I, Al-l14b C2,ZK' Q Al-14c and the diffuse transmitted intensity may be written as I (,-/,, ) l1 130

2. Integrated Reflectance of Half Space for Rayleigh Scattering There is a section on integral formulations in reference (4). For Rayleigh scattering (o, =o), Equation I-151 in that report then reduces to, for k = 0 and 1, ( o A.O -I)I U t U, + t Ao. U, = O A2-1 (, + L.? Aoi) U - 3/1U, (Z 3 I, A,,p) UA = 0 A2-2 where A and U are defined as tI (') = U4 (1A) K(Q) A2-3 A,* P (-,P,)')'IL4 a A2-4 Ar = UV,(U')P(-? )K(/')du' A2-4 K is equivalent to the H function in Chandrasekhar's Radiative Transfer (1) and P is the Legendre polynomial. According to Equation I-157 in reference (5). ULJ(/)= P, (') when?,wo A2-5 By this additional condition Equations A2-1 and A2-2 can be solved for F-o' 41 and 42' Uo I *t4(5 4l)zA2I/ t'O;t 4Alo/ A2-6 u.- (, 3. tt2, A.,,O A) o,. A2-7,,_L (,.W.,o^.., A,* +o)/,. [ 3 A W Ax (, - W-~o o.).,''."^.,M? A2-7 U, - eAo + 3 (2 )M D)Ie A2-8 The K function is,,sI _ zf'c', K(. J, A2 -9? +?, and N (' =1 ~ (-Q'Or Ur (/U'1 A2-10'. 131

The latter can be written in terms of A's using Equations A2-6, A2-7 and A2-8 i (,M) ='o + 4 + + (3 + o 0t Ao. - ) JA; 296OA. +'0 4)o4 A2o Wb *N 2" t2 9 Ai _' G A.. + 6 Xo2 eA2-11 The A's can be cancelled out by the use of Equations (I-149) and (I-152) in the reference (5). N r I" U. (,,) K (,) = Pi (I,) + I T Z (-1) (Jr t t,, ( -/',,'-.o,, A2-12 u, cy, K4rctl P~ (~ + z +'f1 rzo For k = O, multiplying A2-12 by Po(r) then integrating from 0 to 1 gives Apo - 1 t ao Ao,' F A2, A2-13 Similarly for k = 1 and 2. A = - IA., A2-14 3 + x@, As,, = - (,,, Ao. +0 A,0 ) + ( (.o A., +, A,, ) A2-15 Then A2-11 can be simplified with A2-13, A2-14, and A2-15 to give Tk (I)Lr Io + 4 -. (Ai22 - oI t+ j (I I0.) a/c t A2-16 The numerical solution of R was carried out as follows. The K function was first evaluated by the computer using an iterative procedure. Then the moments of K are computed. (, K( )/ ()l /g A2-17 Then A's are computed from the values ofc'as and A2-6, A2-7, and A2-8. However it appears more convenient to use A2-6 and A2-8 to derive the following. 132

AO.. -= 42,41.A2 + )(o2A, A 1 A2-18 A., =.(, +WCt2 A,1 +$O10,jA20 A2-19 Ate = 0(._l coase Aoi + (2 -ef*Ae )01z A2-20. =' )oW *2 ^d, ( -., -,o, ^ (2-21o A.), Finally R can be found by R=| = I _ 1 u, (/,Y) K (,.) =I _1 f -4.A*o + 4-~A~. ( Ao Az,+ Ox.OAA I Ao - 3, Ap. A,)?.]) < (?.) A2 -22 Sample values of the results for c4o= 0.9, 0.6, and 0.3 are tabulated in Table A-1. 133

Table A-1 The K-functions and Related Numerical Values for Computing the Integrated Reflectance of Half Space, Rayleigh Scattering o =0.90,o = o.60 q)o = o..3 p0 K(.0{) R K(o.0) R K(4o) R 0.00 1.0000 o.687 1.0000 0.369 1.0000 0.164 0.10 1.1863 0.630 1.0996 0.308 1.0436 0.128 0.20 1.3144 0.591 1.1578 0.273 1.0667 0.110 0.30 1.4214 0.559 1.2018 0.247 1.0832 o.og98 0.40 1.5147 0.531 1.2371 0.227 1.0960 0.088 0.50 1.5976 0.507 1.2665 0.211 1.1063 0.081 0.60 1.6722 0.485 1.2915 0.197 1.1148 0.075 0.70 1.7400 0.466 1.3131 0.186 1.1220 0.071 0.80 1.8020 0.449 1.3319 0.176 1.1282 0.067 0.90 1.8591 0.433 1.3486 0.168 1.1335 o.o64 1.00 1.9117 0.418 1.3634 0.161 1.1382 0.061 oo 4of all c C3 0A9o -AoI A2o -A2, 0.90 1.5582 0.8493 0.5868 0.4486 1.5198 0.8223 0.0295 0.1138 0.60 1.2442 0.6481 0.4391 0.3321 1.2252 0.6351 0.0203 0.1165 0.30 1.0969 0.5579 0.3744 0.2817 1.0889 0.5525 0.0093 0.1211 134

3. Soli*tin for the Six-flux Formulation of the Slab Problem The following solution of equations II-7a to II-8 subJect to the boundary conditions I1 = 1 and I3 = 0 at Z = 0 and kI = I4 = 0 at Z = v1 is reproduced directly from a previous report to avoid errors in transcription. For this reason the symbolism is slightly different. In equations A3-1 to A3-31 T is used in place of v1 for the optical thickness of the dispersion, Q is used in place of 0O for the angle the incident radiation makes with the normal to the slab, but as in this report Z refers to the normal optical distance through the slab. All of the other symbols such as the A's, a's, b's,a, p, etc. are defined as used: the definitions of these s bols appl only to equations A3-1 to A3-31 and have no relation to the definition of the symbols in the NOMENCLATURE or in the rest of this report. I1 = All e-PZ + A2 e-Pl(T-Z)+ A13 e-P2Z + A e-P2(T-Z) A3.-1 A A A A I = e-PlZ+ 5 e-Pl(T-Z)+ A3: e-PPZ+ A-2 e 14 =~ A4 e P1Z+- A1 e-pl(T-Z)+ _ eZ-PZ+4 ep2(T-Z) A3-4 A A A wrhere: A = (a1b3-a3bl )2 (ab3-a3b2 )2 e 2p1T_(alb4,a4b1 )2e -22T +2(alb-aI bl ) (a3b4-a4b3 )e-(P+P2 )T+ (a2b4-a4b2 )2e2 (P1+P2 )T A3T. All s alb3(a)b3-a3b ) -(ab4b(ab-ab)ePT b2(a3ba4b3 )e (P+P )T A3.6 Al= a2b3 ( alb3-a3b -a2b4(abb4-az4b )e -2P2T+a2b2 (a3b4-a4b3)e - (P+P2)T A5.7 A31 - blb3(aBlb3sa3bl )-b b4 (lb,-a,,b, )"2PaBT+blb2(a3b4-a4b3 )e-(p+p2)T A3"8 * bab3 (ab3 a3bl ).babP4 ( b 4bl o -2P2Tab2 (a 3b4 _a~bs )5e(P1+P2 )T A 135

Aiz = -a2b3(a2b3-a3b2)e PlT-a2b,(a3b4-a4b3)eiPzT +a2b 42b4-a4b2)eo (P1+2P2)T A-10 A22 = -alb3(a2b3-a3b2)e -p1T-albl(a3b4 —aib3)eP2T+-b4(a2b4-a4b2)e-((PL+2P2)T 51 A32 = -b2b3(a2b3-a3b2z)e P1T-b2bl(a3b4-a4b3)e-P2T+b2b4(aab4-a4b2 )e (Pj+2p2)T 13-12 A42 = -blb3(a2b3-a3b2)e-piT-bl2(a3b4-a4b3)e-P2T+blb4(avb4-a4b2 )e (pl+2-Pa)T 3-13 A13 = -a3bl(alb3-a3bl)+a3b2(aa2b-a3bz)e-2PiT+a3b4(alb2-a2b1 )e'(Pl+P2)T A-314 A23 = -a4bl(alb3-a3bl )+a4b2(a2b3-a3b2)e 2P1T+a4b4(alb2-a2b1)e(P1+P2)T 3-15 A33 = -b3bl(albs-a3bl)+b3bz(a2b3-a3b2 )e'2PT+b3b4(alb2-a2b~)e (P~+P2)T A35-16 A43 = -b4bl(alb3-a3sb)+b4b2(a2b3-as3b)e-2plT+ b4=(alb2-abi )e-(Pl+P2)T 0A-17 A14 = -a4b3(alb2-a2b1)eiP1T+a4b (alb4-a4bl )e-ePT-a4ba (a2b4-a4b2)e-(2PI+PaTr A_1 8 A24 = -a3b3(alb2-a2bi)e PlT+a3bl (alb4-a4b1 )e'P2T-a3b2(a2b4-a4bz)e'(2p1+p2)T A3-19 A34 = -b4b3(alb2-a2bl )e'PT+b4bl(ajb4-a4bl)e'P2T-b4b2(anb4.a4b2)e (2Pi+PT A3-20 A44' -b32(alb2-a2bl)e'PT+b3bl (a&b4-a4bl)e'P2T-b3b2(a2b4-a4b2)e'(2pl+p2)T A3-21 in which, Pi (J[- J -s'2;2 e4 J J C' sec, ase. 3-22136

r'/2 (c -cf~l2 sec 24 (Cl -C2 )2 a (Cl+C2) sec 9 + p (Cl+C2) sec 0 a's (C1+C2) sec 9 - P A3-25 (C1+C2) sec 9 a3 (Cl+c2) sec 9 + P2 A5-26 (C1,+C2) sec 9 (C1.+C2) sec o - p2 A3-27 (Clz C2) sec 9 b (C12-C22) sec2 - _ p2 (cl+c2) csc 9 + Pl A~28 2(Cl+C2) C3 seca (C1+C2) csc 9.b (C12-C22) sec2 8 - Pl2 (C1+C2) csc 9 - P A3-29 2(C1+C2) C3 sec 9 (C1+C2) csc G b = (ClZ'Ca) sec2) - pz2 (CC) csc 6 + P2 A3-30 2(C1*C2) C3 sedc2 (C1+C2) csc 9 and b (C2'-C22) sec2-9 p22 (C1+C2) CSC 9 - P2 A3-31 2(Cl+C2) C3 sec2 9 (Cl+C2) CSC 9 137

APPENDIX B Tables of Computed Functions The tabulated results obtained from the computer are presented in this appendix. In the first portion of the appendix the complete set of *~t) and c;(9) functions, the integrated reflectance R, the diffuse portion of the integrated transmission TD, and the total integrated transmission T are given as a function of [i for each of the 36 sets of parameters used. The results are presented as one long table (Table B-O) which is subdivided into 36 problems. The values of the parameters corresponding to each problem are listed in Table I-2. In the last portion of the appendix (Tables B-1 to B-5) tabular results are given which compare the exact solution with results obtained from approximate models. 138

Table B-O The [ and ~ Functions, Integrated Reflectances, Diffuse Portion of the Integrated Transmission, and Total Integrated Transmission for 36 sets of Parameters (Listed in Table I-2). PROBLEM 1 TAU = 0.050 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.0000 Al = 0.000 A2 = 0.0000 F = 0.500 P = 0.3333 NUMBER OF ITERATIONS WAS 2 01 0 TR T 0*04691 1.0546 0.3838 0.299 0.284 0.628 0.23077 1.0726 0.9049 0,078 0.098 0.903 0.50000 1.0783 0.9824 0.042 0.039 0.944 0.76923 1.0800 1.0138 0.028 0*025 0.962 0.95309 1.0806 1.0252 0,023 0.019 0,968 PROBLEM 2 TAU = 0.250 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 10000 Al = 0.0000 A2 = 00000 F = 0,500 P = 0.3333 NUMBER OF ITERATIONS WAS 4 o o 1 to 0O R TD T 0*04691 1.0797 0.0467 0,498 0.360 0.365 0.23077 1,1993 0.4990 0O300 0.268 0*607 0.50000 1.2413 0,8317 0,173 0.161 0.768 0.76923 1.2611 0.9501 0.125 0.104 0*827 0.95309 1.2694 1.0001 0.105 0,085 0.854 PROBLEM 3 TAU = 0.500 ALBEDO FOR SINGLE SCATTERING = 0*90 AO = 1.0000 Al = 0.0000 A2 = 0*0000 F 2 0.500 P = 03333 NUMBER OF ITERATIONS WAS 4 to 00 R TD T 0.04691 1.0842 0.0275 0.548 0,286 0.286 0.23077 1*2530 0*2636 0.417 0,312 0.426 0.50000 1.3465 0o6497 0o279 0.239 0,607 0.76923 1.3896 0.8462 0.209 0,171 0.693 0.95309 1.4089 0,9249 0.179 0,136 0.728 139

PROBLEM 4 TAU = 1.000 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.0000 Al = 0.0000 A2 = 0.0000 F = 0.500 P = 0)3333 NUMBER OF ITERATIONS WAS 6 to 0 R TD T 0.04691 1.0888 0.0165 0.598 0.197 0-197 0.23077 1.2915 0.1231 0*504 0.258 0o271 0.50000 1l4541 0.4012 0.394 0.263 0.398 0.76923 1.5427 0.6476 0.316 0*225 0.498 0.95309 1*5852 0.7553 0.279 0.186 0.536 PROBLEM 5 TAU = 1.500 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.0000 Al = 0.0000 A2 = 0.0000 F = 0.500 P = 0.3333 NUMBER OF ITERATIONS WAS 9 0 0 to'o R TD T 0.04691 l.0910 0.0123 0.622 0.152 0,152 0.23077 1.3063 0.0845 0.537 0.205 0.207 0O50000 1.5009 0.2769 0.444 0.248 0.298 0.76923 1.6202 0.5158 0,370 0.253 0.395 0.95309 1.6774 0.6565 0.330 0.242 0.449 PROBLEM 6 TAU = 2.000 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.0000 Al = 0.0000 A2 = 0.0000 F = 0.500 P = 0.3333 NUMBER OF ITERATTONS WAS 9 0t 0 R TD T 0.04691 1.0919 0.0098 0*632 0.129 0.129 0.23077 1.3122 0.0665 0.552 0.174 0.174 0.50000 1.5204 0.2145 0.465 0o233 0.251 0*76923 1.6567 0.4370 0.395 0.271 0.345 0.95309 1.7224 0.6056 0.355 0*290 0.413 140

PROBLEM 7 TAU = INFINITY ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1*0000 Al = 0*0000 A2 = 0.0000 F = 0.500 P = 0.3333 NUMBER OF ITERATIONS WAS 16 04 10 003R TD T 0.04691 1.0936 0.0000 0.654 0.000 0.000 0*23077 1.3236 0.0000 0.581 0.000 0.000 0.50000 1.5554 0.0000 0.507 0.000 0.000 0. 76923 1.7289 0.0000 0.452 0.000 0.000 0.95309 1.8259 0.0000 0.422 0.000 0.000 PROBLEM 8 TAU = 0.050 ALBEDO FOR SINGLE SCATTERING = 0.60 AO = 1.00O A1 = 0.0000 A2 = 0.0000 F = 0.500 P = 0.3333 NUMBER OF ITERATIONS WAS 3 t 00~ R TD0R T 0.04691 1.0355 0.3696 0,1931 0*1817 0.526 0-23077 1.0477 0*8694 0*0529 0.0665 0.871 0O50000 1.0516 0.9557 0.0273 0.0271 0.932 0.76923 1.0527 0.9877 0.0182 0.0173 0.954 0995309 1.0531 0.9991 0.0148 0,0132 0.962 PROBLEM 9 TAU = 0,500 ALBEDO FOR SINGLE SCATTERING = 0.600 AO = 1000 Al = 0.0000 A2 - 0.0000 F = 0.500 p = 0.3333 NUMBER OF ITERATIONS WAS 5 I \o 00' R TD T 0.04691 1.0507 0.0145 0.318 0.154 0-154 0*23077 1*1448 0O1939 0.232 0.168 0.283 0.50000 1.1958 0*5221 0.154 0.134 0.502 0.76923 1l2186 0.7046 0.114 0.101 0.623 0.95309 1.2283 0.7827 0.097 0.083 0.675 141

PROBLEM 10 TAU = INFINITY ALBEDO FOR SINGLE SCATTERING = 0.60 AO = 1.0000 Al = 0.0000 A2 = 0.0000 F = 0.500 P = 0.3333 NUMBER OF ITERATIONS WAS 7 0 Co 0~ R TD T 0.04691 1.0520 0.0000 0.335 0.000 0.000 0.23077 1.1587 0.0000 0,267 0.000 0.000 0.50000 1.2456 0.0000 0.212 0.000 0.000 0.76923 1.3009 0.0000 0O177 0,000 0.000 0.95309 1.3288 0.0000 0.158 0.000 0.000 PROBLEM 12 TAU = 0.500 ALBEDO FOR SINGLE SCATTERING = 0.30 AO 1.O0000 Al = 0*0000 A2 = 0.O0000 F = 0.500 P = 0.3333 NUMBER OF ITERATIONS WAS 4 0oB R TD T 0.04691 1.0232 0.0059 0.1416 0.0636 0.0636 0.23077 1*0632 0.1469 0.0999 0.0696 0.1842 0.50000 1.0844 0.4323 0.0658 0.0569 0.4248 0.76923 1.0936 0.6000 0*0484 0*0439 0.5659 0.95309 1.0974 0.6743 0,0410 0.0366 0*6284 PROBLEM 11 TAU = 0.050 ALBEDO FOR SINGLE SCATTERING = 0.30 AO = l.0000 Al = 0.0000 A2 = 0.0000 F = 0.500 P = 0.3333 NUMBER OF ITERATIONS WAS 2 0o 00 R TD T 0.04691 1.0172 0.3564 0.0935 0.0874 0.431 0.23077 1.0237 0.8365 0,0265 0.0337 0.839 0.50000 1.0252 0.9307 0.0132 0.0139 0.918 0.76923 1.0256 0.9624 0.0087 0.0087 0.945 0.95309 1.0258 0,9740 0.0071 0.0067 0.956 142

PROBLEM 13 TAU = INFINITY ALBEDO FOR SINGLE SCATTERING = 0.30 AO = 1.0000 Al = 0.0000 A2 = 0.0000 F = 0.500 P = 0.3333 NUMBER OF ITERATIONS WAS 5 t\o 00 R TD T 0.04691 1.0233 0.0000 0.1438 0.0000 0.0000 0.23077 1.0663 0.0000 0.1079 0.0000 0.0000 0.50000 1.0975 0.0000 0.0817 0.0000 0.0000 0.76923 1.1159 0.0000 0.0663 0.0000 0.0000 0.95309 1.1248 0.0000 0.0588 0.0000 0.0000 PROBLEM 14 TAU = 0.050 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.o000 Al = 0.0000 A2 = 0.5000 F = 0.500 P 0.400 NUMBER OF ITERATIONS FOR M = 0 1, AND 2 WERE 2, 2, AND 2 0 0 0 1 1 2 L 3 *1 *2 o1 i2'22 0~ 01 0.04691 1.0590 0.0000 -0.5158 0.0000 0.1409 3.0162 0.3864 0.0000 0.23077 1.0771 0.0000 -0*4452 0.0000 0.6760 2.8703 0.9132 0.0000 0*50000 1*0799 0.0000 -0.1493 0*0000 1.3040 2.2751 0.9835 0.0000 0.76923 1.0756 0.0000 0.3670 0.0000 1.4802 1.2387 1.0093 0*0000 0.95309 1.0705 0.0000 0*8457 0.0000 0.8688 0.2799 1.0155 0.0000.04691 0.1817 0.0000 0.0488 1.0452 0299 0283 0627 0*04691 -O.1817 0.0000 0*0488 1.0452 0.299 0.283 0.627 0.23077 -0.3774 0.0000 0.5447 2.3280 0.077 0.100 0.905 0.50000 -0.1370 0.0000 1.1801 2.0618 0.042 0.038 0.943 0.76923 0.3433 0.0000 1.3871 1.1615 0.028 0.026 0.963 0.95309 0.8015 0*0000 0.8247 0.2639 0.023 0.020 0.969 143

PROBLEM 15 TAU = 0.250 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 10000O Al = 0.0000 A2 = 0.5000 F = 0.500 P = 0.400 NUMBER OF ITERATIONS FOR M = 0, 1, AND 2 WERE 4. 2. AND 3 o 1 1 2 o C1 *1 32 1 *2 *2 00 1 0.04691 1.0857 0.0000 -0.5214 0.0000 0.1410 3.0222 0.0487 0.0000 0.23077 1.2102 0.0000 -0.4690 0.0000 0.6807 2.8997 0.5063 0.0000 0-50000 1.2461 0.0000 -0.1759 0.0000 1.3161 2*3058 0.8379 0.0000 0O76923 1.2525 0.0000 0O3446 0.0000 1.4954 1.2569 0*9439 0*0000 0.95309 1.2471 0.0000 0.8287 0.0000 0.8781 0.2822 0.9849 0.0000 o 1 1 2 II 02 01 02 02 R TD T 0.04691 -0.0099 0.0000 0.0010 0.0257 0.500 0.356 0.361 0O23077 -0.1737 0.0000 0.2340 1*0025 0.299 0.267 0.606 0.50000 -0.1237 0.0000 0.8030 1.4149 0.173 0.163 0.770 0.76923 0.2418 0.0000 1.0843 0.9141 0.125 0.107 0.830 0.95309 0.6315 0.0000 0.6787 0.2182 0.105 0.090 0.859 PROBLEM 16 TAU = 0.500 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1. 000 Al = 0.0000 A2 = 0.5000 F = 0*500 P = 0.400 NUMBER OF ITERATIONS FOR M = O, 1, AND 2 WERE 5, 2. AND 3 o 0 o 1 1 2 o to r1 2 1 2 2 00 1 0*04691 1.0908 0.0000 -0.5217 0*0000 0.1411 3.0222 0.0284 0.0000 0*23077 1.2680 0.0000 -0.4723 0.00 0*6822 2.9054 0.2680 0.0000 0.50000 1.3570 0.0000 -0.1809 0.0000 1*3226 2*3157 0.6554 0.0000 0 76923 1 3863 0.0000 0. 3414 0.0000 1.5049 1.2638 0.8365 0.0000 0.95309 1.3909 0.0000 0.8285 0.0000 0.8841 0.2839 0.9033 0.0000 0 1 1 2 p02 01 02 02 R TD T 0.04691 -0.0018 0.0000 0.0002 0.0049 0.553 0.280 0.280 0-23077 -0.0584 0.0000 0.0831 0.3492 0.419 0.308 0.422 0.50000 -0.0838 0.0000 0.4960 0.8722 0.281 0.236 0.604 0o76923 0.1671 0.0000 0.7948 0.6696 0.212 0.170 0.692 0.95309 0.4789 0.0000 0.5304 0.1700 0*183 0.138 0.729 144

PROBLEM 17 TAU = 1.000 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.0OO00 Al = 0.0000 A2 = 0.5000 F = 0.500 P = 0.400 NUMBER OF ITERATIONS FOR M = O, 1, AND 2 WERE 6, 3, AND 3 0 o o 1 1 2 o I4 to V1 2 i1 *2 2 o 01 0.04691 1.0949 0.0000 -0.5215 0.0000 0.1411 3.0223 0.0162 0.0000 0.23077 1.3044 0.0000 -0.4711 0.0000 0.6826 2.9065 0.1226 0.0000 0.50000 1*4614 0.0000 -0.1772 0.0000 1.3263 2.3196 0.4017 0.0000 0.76923 1.5334 0.0000 0.3485 0.0000 1.5116 1.2671 0.6412 0.0000 0.95309 1.5594 0.0000 0.8387 0.0000 0.8888 0.2848 0.7341 0.0000 02 01 02 02 R TD T 0.04691 0.0010 0.0000 0.0001 0.0010 0.601 0.190 0.190 0*23077 -0.0014 0.0000 0.0122 0.0459 0.504 0.253 0.266 0.50000 -0.0205 0.0000 0.1900 0.3272 0.394 0.259 0.394 0.76923 0.0917 0.0000 0.4247 0.3550 0.317 0.221 0.493 0.95309 0.2876 0.0000 0.3220 0.1019 0.282 0.180 0.530 PROBLEM 18 TAU = 1.500 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.0000 Al = 0.0000 A2 = 0.5000 F = 0.500 P = 0.400 NUMBER OF ITERATIONS FOR M = O, 1i, AND 2 WERE 7, 3. AND 3 o o 1 1 2 o 0 0*04691 1*0968 0.0000 -0.5213 0O0000 0.1411 3.0223 0.0113 0.0000 0.23077 1*3178 0.0000 -0.4699 0.0000 0.6826 2.9066 0.0800 0.0000 0*50000 1.5056 0.0000 -0*1731 0.0000 1.3270 2.3200 0.2678 0.0000 0O76923 1.6080 0,0000 0.3556 0.0000 1.5132 1.2677 0.4952 0.0000 0.95309 1.6491 0.0000 0.8476 0.0000 0.8901 0.2850 0.6242 0.0000 00 01 02 0 R TD T 0.04691 0.0013 0.0000 0.0001 0.0003 0.623 0.142 0.142 0.23077 0*0066 0.0000 0*0027 0.0076 0.537 0.194 0O196 0.50000 0.0076 0.0000 0.0737 0.1223 0.443 0.238 0.288 0.76923 0.0656 0.0000 0.2260 0.1870 0.371 0.241 0.383 0.95309 0.1937 0.0000 0.1945 0.0608 0.334 0.226 0.433 145

PROBLEM 19 TAU = 2.000 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.0000 Al = 0.0000 A2 = 0.5000 F = 09500 P = 0.400 NUMBER OF ITERATIONS FOR M = O, 1, AND 2 WERE 11, 3, AND 3 r-'3o i1 *2 *1 2 2 0 01 0.04691 1.0981 0.0000 -0.5212 0.0000 0.1411 3.0223 0.0092 0.0000 0.23077 1*3261 0.0000 -0.4692 0.0000 0.6826 2.9066 0.0632 0.0000 0.50000 1.5304 0.0000 -0.1709 0.0000 1.3271 2.3201 0.2048 0.0000 0.76923 1*6521 0.0000 0.3598 0.0000 1.5137 1.2678 0.4145 0.0000 0.95309 1.7034 0.0000 0.8528 0.0000 0.8905 0.2850 0.5620 0.0000 o 1 1 2 c t 02 01 02 02 R TD T 0.04691 0.0012 0.0000 0.0000 0.0001 0.636 0.118 0.118 0*23077 0.0077 0.0000 0,0009 0.0018 0.554 0.161 0.162 0.50000 0.0177 0.0000 0.0291 0.0458 0.468 0.215 0o233 0.76923 0.0562 0.0000 0.1200 0.0982 0.400 0.248 0.322 0.95309 0.1438 0.0000 0.1170 0.0361 0.363 0.259 0.382 PROBLEM 20 TAU = INFINITY ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.0000 Al = 0.0000 A2 = 0.5000 F = 0.500 P = 0.400 NUMBER OF ITERATIONS FOR M = 0, 1, AND 2 WERE 16, 3, AND 3 o o 1 1 2 o o o "1'2 *1 " 2',2 0O 1 0-04691 1.0996 0*0000 -0.5210 0.000Q0 0.1411 3-0223 0.0000 0.0000 0.23077 1.3362 0.0000 -0.4677 0 0000 0.6826 2.9066 0.0000 0. 0000 0-50000 1.5623 0.0000 -0.1661 0.0000 1.3271 2.3207 0.0000 0.0000 0.76923 1*7192 0.0000 0.3701 0.0000 1.5138 1.2679 0,0000 0.0000 0.95309 1.7999 0.0000 0.8677 0.0000 0.8906 0.2850 0.0000 0.0000 o 1 1 2 102 01 02 02 R TD T 0.04691 0.0000 0.0000 0.0000 0.0000 0.656 0.000 0.000 0.23077 0.0000 0.0000 0O0000 0O0000 0.580 0.000 0.000 0.50000 0.0000 0.0000 0.0000 0.0000 0.507 0.000 0.000 0.76923 0.0000 0.0000 0.0000 0.0000 0.453 0.000 0.000 0.95309 0.0000 0.0000 0.0000 0.0000 0.424 0.000 0.000 146

PROBLEM 21 TAU = 1,000 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.0000 Al = 1.0000 A2 = 1.8165 F = 0.750 P = 0.5504 NUMBER OF ITERATIONS FOR M = O, 1, AND 2 WERE 4. 11, AND 5 o o o 1 1 2 4 to 0 41 42 ti *2 2 01 0.04691 1.1075 0.0211 -0.5329 1.0266 0.1242 3.1069 0.0164 0.0089 0.23077 1.3101 0.1331 -0.5012 1.0228 0*6458 3.1140 0.1332 0.0665 0.50000 1.3954 0.3549 -0.2024 0.8819 1.3221 2.5446 0.4580 0.2319 0*76923 1*3537 0.6124 0.3586 0.6120 1.5497 1.4059 0,7536 0*4677 0*95309 1.2821 0.7956 0.8892 0.2757 0.9220 0.3175 0.8755 0.6381 02 01 0 02 R TD T 0.04691 0O0006 0O0027 0*0044 0.0052 0.563 0.203 0.203 0.23077 -0*0034 0.0397 0,0494 0.0788 0.438 0.290 0.303 0.50000 -0.0162 0.2019 0.2957 0.4097 0.303 0.337 0.472 0*76923 0.1277 0.3034 0.5826 0.4311 0.213 0.333 0.605 0.95309 0.3629 0.2116 0.4639 0.1220 0.173 0.307 0.657 PROBLEM 22 TAU = 1.000 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.0000 Al = 1.0000 A2 = 1-0490 F = 0.750 P = 0.4821 NUMBER OF ITERATIONS FOR M = O, 1, AND 2 WERE 5, 5, AND 4 o 0 0 1 1 2 o u *to t1 2 *1 *2 42 0 01 0*04691 1*0984 0*0207 -0*5263 1.0276 0.1228 3.0559 0.0185 0.0102 0.23077 1.2944 0.1295 -0.4847 1.0351 0.6288 2*9865 0.1422 0.0727 0.50000 1*3958 0*3481 -0.1923 0.9181 1.2706 2.4050 0.4702 0.2396 0.76923 1.3865 0.6065 0.3426 0.6611 1.4793 1.3195 0.7712 0.4682 0.95309 1l3458 0.7926 0,8450 0.3063 0.8782 0.2971 0.9051 0.6317 I 0 01 01 02 R TD T 0.04691 0.0011 0.0023 0.0036 0.0024 0.564 0.222 0.222 0o23077 -0.0002 0.0334 0.0392 0.0574 0.445 0*305 0.318 0.50000 -0.0134 0.1777 0*2572 0.3569 0.310 0.341 0.476 0.7692'3 0,1145 0.2583 0.5119 0.3825 0.218 0.325 0.598 0.95309 0.3290 0.1647 0,3853 0.1092 0.175 0.295 0.645

PROBLEM 23 TAU = 1.000 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.0000 Al = 1.0000 A2 = 0.0000 F = 0.750 P = 0.3889 NUMBER OF ITERATIONS FOR M = 0 AND 1 WERE 6 AND 4 o o o 1 2 o o p \o *1 V2 4i 2 2 0 01 0.04691 1.0856 0.0208 0.0000 1.0288 0.0000 0.0000 0.0198 0*0112 0*23077 1.2689 0.1282 0.0000 1.0498 0.0000 0.0000 0.1464 0.0767 0.50000 1.3856 0.3445 0.0000 0.9607 0.0000 0o0000 0.4720 0.2433 0.76923 1.4151 0.6033 0.0000 0.7185 0.0000 0.0000 0.7817 0.4675 0.95309 1.4126 0.7913 0.0000 0.3424 0.0000 0.0000 0.9396 0.6300 1' 02 01 02 02 R TD T 0.04691 0*0000 0.0018 0.0000 0.0000 0.559 0.238 0.238 0.23077 0.0000 0.0269 0.0000 0.0000 0*448 0.316 0,329 0.50000 0*0000 0.1540 0.0000 0.0000 0.315 0.345 0.480 0.76923 0.0000 0.2180 0.0000 0.0000 0.220 0.324 0.596 0.95309 0.0000 0.1312 0.0000 0.0000 0*174 0.296 0.646 PROBLEM 24 TAU = 0.050 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1*0000 Al = 1.7321 A2 = 1.0000 F = 0*9330 P = 0.4821 NUMBER OF ITERATIONS FOR M = 0, 1, AND 2 WERE 2, 4, AND 3 o o o 1 1 2 o 0.04691 1*0622 0*0336 -0.5187 1*0385 0.1182 3.0397 0.3903 0.0300 0*23077 1.0745 0.2166 -0.4485 1.0213 0.6486 2.9005 0.9285 0.2142 0.50000 1.0643 0.4890 -0.1510 0.9086 1.2790 2,3066 1.0007 0.4807 0.76923 1.0438 0.7628 0.3689 0.6676 1.4635 1.2528 1.0309 0.7542 0*95309 1*0262 0.9499 0.8511 0.3151 0.8916 0.2812 1.0390 0.9415 4 02 01 02 02 R TD T 0.04691 -0.1837 0.3718 0.0703 1.0601 0.287 0.295 0.639 0.23077 -0.3821 0.8605 0.5889 2.3706 0.060 0.117 0.922 0*50000 -0.1373 0.8387 1.2127 2.0874 0.023 0*055 0*960 0.76923 0.3473 0.6440 1.4139 1.1753 0.009 0.045 0o982 0o95309 0.8097 0.3108 0.8396 0.2670 0.004 0.040 0.989 148

PROBLEM 25 TAU = 0.250 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.0000 Al = 1.7321 A2 = 1.0000 F = 0.9330 P = 0.4821 NUMBER OF ITERATIONS FOR M = 0, 1, AND 2 WERE 6, 10, AND 3 uto 2 o21 i 2 2 01 0.04691 1.0890 0.0251 -0.5256 1.0522 0*1068 3.0524 0.0541 0.0186 0*23077 1.1976 0.1758 -0.4792 1.0859 0.5918 2.9633 0.5391 0.1543 0.50000 1.1860 0.4514 -0.1851 0.9716 1.2777 2.3655 0.9148 0.4269 0.76923 1.1337 0.7363 0.3424 0.7104 1.4325 1.2914 1,0700 0.7177 0*95309 1.0806 0.9348 0.8344 0.3328 0.8507 0.2901 1,1597 0.9336 4 02 0012 R TD T 0.04691 -0*0116 0.0303 0.0248 0*0378 0*466 0.394 0.396 0.23077 -0*1771 0.4320 0*3258 1.0475 0.238 0.335 0.669 0.50000 -0.1199 0.6730 0.9280 1.4696 0.097 0.256 0.853 0.76923 0.2664 0.5879 1.1875 0.9455 0.043 0.226 0.933 0.95309 0.6743 0.3072 0.7458 0.2255 0.018 0.230 0.979 PROBLEM 26 TAU = 0.500 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.0000 Al = 1.7321 A2 = 1.0000 F = 0.9330 P = 04821 NUMBER OF ITERATIONS FOR M = 0, 1, AND 2 WERE 5, 10, AND 4 o o o 1 1 2 0to'r1 t2 1 t2 o 2 o 0*04691 1.0926 0.0235 -0.5262 1*0532 0.1057 3.0527 0.0337 0.0156 0.23077 1.2466 0.1545 -0.4854 1.1044 0.5701 2.9763 0.3067 0.1173 0.50000 1.2615 0.4199 -0.1968 1.0005 1.1953 2.3878 0.7653 0.3632 0.76923 1.1958 0.7115 0.3311 0.7317 1.4104 1.3067 1.0156 0.6422 0.95309 1.1267 0.9164 0.8262 0.3421 0.8419 0.2939 1.1318 0.8434 C1 22 01 02 0 TD T 0.04691 -0.0023 0.0014 0.0160 0.0098 0.500 0.331 0.331 0.23077 -0.0588 0.1901 0.1720 0*3764 0.333 0.392 0.507 0.50000 -0.0779 0.4782 0.6520 0.9228 0.163 0.355 0.723 0*76923 0,1919 0.4925 0.9454 0.7034 0*078 0.311 0*833 0.95309 0.5234 0.2768 0.6291 0.1781 0.041 0.293 0.885

PROBLEM 27 TAU - 1.000 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1*0000 Al = 1.7321 A2 = 1*0000 F = 0*9330 P = 0.4821 NUMBER OF ITERATIONS FOR M = 09 1, AND 2 WERE 5, 10, AND 4 o o 1 1 2 to I1 t2 1 t2 2 00 0.04691 1.0951 0.0223 -0.5262 1.0537 0.1051 3.0528 0.0212 0.0119 0.23077 1.2736 0.1408 -0.4859 1.1107 0.5618 2.9790 0.1630 0.0852 0.50000 1.3347 0.3832 -0.1988 1.0178 1.1727 2.3969 0.5354 0.2784 0.76923 1.2749 0.6723 0.3281 0.7470 1.3913 1.3145 0.8972 0.5418 0.95309 1.1966 0.8816 0.8239 0.3493 0.8334 0.2959 1.0802 0.7334 o 1 1 21 4 012 01 02 02 R TD T 0.04691 0.0016 0.0053 0.0080 0.0023 0.530 0.259 0.259 0.23077 0.0028 0.0544 0.0677 0.0563 0.396 0.359 0.372 0.50000 -0.0043 0.2286 0.3185 0.3540 0.239 0.418 0.553 0.76923 0*1298 0*3207 0.5855 0.3758 0.131 0.419 0.692 0.95309 0.3496 0.2072 0.4387 0.1085 0.079 0.400 0.750 PROBLEM 28 TAU = 1.500 ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.0000 A1 = 1.7321 A2 = 1.0000 F = 0.9330 P = 0.4821 NUMBER OF ITERATIONS FOR M = 0, 1, AND 2 WERE 7, 11, AND 3 o o 0 1 1 2 o uto t1 t2 ti t2 0 0.04691 1.0972 0.0212 -0.5262 1.0537 0.1050 3.0528 0.0171 0.0100 0.23077 1.2882 0.1331 -0.4855 1.1115 0.5606 2.9792 0.1231 0.0702 0.50000 1.3974 0.3589 -0o1968 1.0214 1.1674 2.3981 0.4058 0.2253 0.76923 1.3398 0.6369 0.3314 0.7513 1.3850 1.3159 0,7702 0.4588 0.95309 1.2661 0.8436 0o8278 0.3516 0.8302 0.2963 0.9945 0o6437 U 0~ 0 1 02 02 R TD T 0.04691 0.0020 0.0027 0.0043 0.0007 0.551 0.215 0*215 0.23077 0.0118 0.0234 0.0341 0.0118 0.426 0.303 0.305 0.50000 0.2515 0.1108 0.1606 0.1347 0.286 0.397 0.447 0.76923 0.1036 0.1991 0.3554 0.2022 0.176 0.446 0.588 0.95309 0.2542 0.1496 0.3044 0.0653 0.118 0.456 0.663 150

PROBLEM 29 TAU = INFINITY ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.0000 Al = 1.7321 A2 = 1.0000 F = 0.9330 P = 0.4821 NUMBER OF ITERATIONS FOR M = 09 1, AND 2 WERE 21, 6, AND 4 o o o 1 1 2 o o ~to t1 t2 1 2 2 / 01 0.04691 1.1010 0.0190 -0.5257 1.0538 0.1049 3.0528 0.0000 0.0000 2.23077 1.3142 0.1180 -0.4825 1.1116 0.5604 2.9792 0,0000 0.0000 0.50000 1.4608 0.3115 -0.1869 1.0222 1.1659 2.3983 0,0000 0.0000 0.76923 1.4964 0.5451 0.3511 0.7530 1.3823 1.3163 0,0000 0.0000 0.95309 1.4740 0.7213 0.8547 0.3527 0.8283 0.2965 0,0000 0.0000 00 01 1 2 R T T Yt 2 01 02 02 TD T 0*04691 0.0000 0.0000 00000 0 0000 0.595 0.000 0.000 0.23077 00000 0000 000 0 0.0000 0.489 0.000 0.000 0.50000 0.0000 0.0000 0.0000 0.0000 0.377 0.000 0.000 0-76923 0.0000 0.0000 0.0000 0.0000 00292 0.000 0.000 0O95309 0.0000 0.0000 0.0000 0.0000 0.244 0.000 0,000 PROBLEM 30 TAU = 0.050 ALBEDO FOR SINGLE SCATTERING = 0.60 AO = 1o0000 Al = 1-7321 A2 = 1.0000 F = 0.9330 P = 0.4821 NUMBER OF ITERATTONS FOR M = O 1, AND 2 WERE 2, 3, AND 3 o o o 1 1 2.0 0 ~ ~2 ~1 *2 *2 0o 0 0.04691 1.0002 O0Oi83 -0.5109 1.0247 0.1261 3.0240 0.3734 0.0249 0 23077 1,0492 0.2212 -0.4387 1.0051 0.6568 2.8803 0.8846 0.2046 0.50000 1.0421 0.4927 -0.1419 0.8937 1.2862 2.2836 0.9689 0.4716 0.76923 1.0287 0,7649 0.3755 0.6575 1.4676 1.2134 0.9992 0.7431 0.95309 1.0172 0.9509 0.8552 0.3108 0.8630 0.2791 1.0083 0.9290 1 1 2 02 0)1 02 02 R TD T 0.04691 -0.1790 0.3618 0.0626 1.0501 0.184 0.188 0.531 0.23077 -0,3659 0.8338 0.5738 2.3420 0.042 0.080 0,887 0.50000 -0.1291 0.8214 1.2013 2.0703 O.O15 0.038 0.943 0.76923 0.3528 0.6286 1.4027 1.1661 0.006 0.029 0,966 0.95309 0.8127 0.3025 0.8327 0.2649 0.003 0.026 0.974 151

PROBLEM 31 TAU = 0.500 ALBEDO FOR SINGLE SCATTERING = 060 AO = 1.0000 Al 1.7321 A2 = 1.0000 F = 0.9330 P = 0.4821 NUMBER OF ITERATIONS FOR M = 0, 1, AND 2 WERE 4, 4, AND 3 o 0 0 1 1 2 0 to t1 2 1 2 2 0 01 0.04691 1.0558 0.0332 -0.5145 1.0331 0.1189 3.0322 0.0166 0.0077 0.23077 1.1385 0.1894 -0.4584 1.0523 0.6123 2.9283 0,2128 0.0737 0.50000 1.1435 0.4574 -0.1660 0.9453 1.2398 2.3389 0.5810 0.2827 0.76923 1.1049 0.7390 03567 0.6924 1.4399 1.2776 0,7966 0.5394 0.95309 1.0662 0.9338 0.8444 0.3252 0.8533 0.2871 0.8951 0.7262 o 1 1 2 I I202 01 02 0R TD T 0O04691 -00010 0O0077 0.0091 0O0063 0O293 0.165 0,165 0.23077 -0.0526 0,1566 0O1321 0.3578 0O180 0.204 0.319 0.50000 -0.0593 0.4118 0.5804 0.8883 0.085 0.197 0.565 0O76923 0O2020 0O4266 0O8720 0O6804 0.040 0.178 0.701 0.95309 0.5251 0.2327 0.5743 0.1726 0.021 0.169 0.761 PROBLEM 32 TAU = INFINITY ALBEDO FOR SINGLE SCATTERING = 0.60 AO = 1.0000 A1 = 1.7321 A2 - 1.0000 F = 0.9330 P = 0.4821 NUMBER OF ITERATIONS FOR M = 0, 1, AND 2 WERE 8, 5, AND 3 o o0 1 1 2 u Jto t* i11 22 *2 0t 01 0.04691 1.0568 0.0327 -0.5148 1.0332 0.1187 3.0322 0.0000 0.0000 0.23077 1.1499 0.1833 -0.4583 1.0551 0.6084 2.9300 0.0000 0.0000 0.50000 1.1795 0*4381 -0.1652 0.9540 1o2279 2.3452 0.0000 0.0000 0.76923 1.1517 0.7134 0.3589 0.7003 1.4298 1.2833 0,0000 0.0000 0.95309 1.1128 0.9074 0*8483 0*3286 0.8493 2.8867 0,0000 0000 0 1 1 2. 02 01 02 02 R TD T 0.04691 0.0000 OO0000 00000 0.0000 O 303 0.000 0.000 0-23077 0.0000 0.0000 0.0000 0o0000 0O206 0.000 0.000 0O50000 00000 00000 0.0000 O00000 O 0124 0.000 0.000 076923 0.0000 0.0000 0.0000 0.0000 0.073 0.000 0.000 0o95309 00000 00000 00000 00000 0.048 0.000 0.000 152

PROBLEM 33 TAU = 0.050 ALBEDO FOR SINGLE SCATTERING = 0*30 AO = 1.0000 Al = 1.7321 A2 = 1.0000 F = 0.9330 P = 0.4821 NUMBER OF ITERATIONS FOR M = 0, 1, AND 2 WERE 2, 2, AND 2 o 0 o 1 1 2 o o 4 to 0 1 12'1 i2'2 o 1 0*04691 1.0195 0.0427 -0.5036 1.0115 0.1335 3.0086 0.3581 0.0204 0.23077 1*0243 0.2260 -0.4293 0.9889 0.6652 2.8602 0.8435 0.1952 0.50000 1.0206 0*4964 -0.1333 0*8796 1*2928 2.2667 0.9371 0.4622 0*76923 1.0146 0.7671 0.3817 0.6480 1.4712 1*2341 0,9681 0.7321 0.95309 1.0086 0.9520 0.8590 0.3066 0.8643 0.2769 0.9785 0.9167 1 1 2 u02 01 02 0P2 R TD T 0.04691 -0.1748 0.3526 0.0553 1.0404 0.090 0.089 0.435 0o23077 -0.3513 0.8080 0.5581 2.3141 0.021 0.040 0.845 0.50000 -0*1216 0.8028 1.1889 2.0531 0.008 0.019 0.924 0.76923 0.3581 0.6137 1.3923 1.1569 0.003 0.015 0.952 0.95309 0.8157 0.2948 0.8268 0.2629 0.001 0.013 0.961 PROBLEM 34 TAU = 0.500 ALBEDO FOR SINGLE SCATTERING = 0.30 AO = 1.0000 Al = 1*7321 A2 = 1*0000 F = 0.9330 P = 0.4821 NUMBER OF ITERATIONS FOR M = 0, 1, AND 2 WERE 4, 3, AND 3 o o o 1 1 2 o 0 o'1 *2 1'2 *2 0 01 0,04691 1.0257 0.0408 -0.5052 1.0152 0.1304 3.0125 0.0064 0.0030 0.23077 1.0598 0.2135 -0*4374 1*0093 0.6460 2.8831 0.1539 0.0456 0.50000 1.0600 0O4831 -0*1432 0.9015 1.2733 2.2931 0.4567 0.2259 0.76923 1.0421 0.7579 0*3743 0.6624 1.4604 1.2504 0.6399 0.4624 0.95309 1.0249 0.9464 0.8553 0.3124 0.8609 0.2808 0.7238 0.6369 02 01 012 2 2 R TD T 0,04691 -0.0004 0.0033 0.0040 0.0030 0.131 0.064 0.063 0.23077 -0.0496 0.1312 0.1015 0.3409 0*075 0*083 0.198 0.50000 -0.0502 0.3602 0.5241 0.8568 0.034 0.084 0.452 0.76923 0.2046 0.3757 0.8164 0.6591 0.015 0.079 0.601 0.95309 0.5202 0.2030 0.5400 0,1674 0*008 0.076 0.668 153

PROBLEM 35 TAU = INFINITY ALBEDO FOR SINGLE SCATTERING = 0.30 AO = 1.0000 Al = 1.7321 A2 = 1.0000 F = 0.9330 P = 0.4821 NUMBER OF ITERATIONS FOR M = 0, 1, AND 2 WERE 4, 3, AND 3 o o o 1 1 2 o o 0o *1 2*21 2 2 0 01 0*04691 1.0258 0*0407 -0.5052 1.0152 0.1304 3.0125 0.0000 0.0000 0*23077 1.0620 0.2124 -0.4374 1.0102 0.6448 2.8838 0.0000 0.0000 0.5000"0 1.0670 0.4793 -0.1433 0.9042 1.2696 2.2959 0.0000 0.0000 0.76923 1.0499 0.7538 0.3744 0.6646 1.4578 1.2530 0,0000 0.0000 0.95309 1.0313 0.9427 0.8559 0.3132 0.8602 2.8149 0.0000 0.0000 o 1 1 2 C 02 01 02 02 R TD T 0.04691 0.0000 0*0000 0.0000 0.0000 0.133 0.000 0.000 0.23077 0.0000 0*0000 0.0000 0.0000 0.080 0.000 0.000 0.50000 0.0000 0.0000 0.0000 0.0000 0.042 0.000 0.000 0.76923 0.0000 0.0000 0*0000 0*0000 0.020 0.000 0,000 0.95309 0.0000 0.0000 0,0000 0.0000 0.011 0.000 0.000 PROBLEM 36 TAU = INFINITY ALBEDO FOR SINGLE SCATTERING = 0.90 AO = 1.0000 Al = 1.7346 A2 = 2*2415 F = 0.9336 P = 0.5821 NUMBER OF ITERATIONS FOR M = O, 1. AND 2 WERE 2, 5, AND 5 (NOTE —FOR M = 0 INITIAL FUNCTIONS WERE THOSE OF REFERENCE 4) o o o 1 1 2 o o to *1 t2 t 2 2 o 01 0.0469'1 1.1185 0.0185 -0.5370 1.0516 0.1080 3.1376 0.0000 0.0000 0.23077 1.3550 0.1179 -0.5124 1.0896 0.5911 3.1955 0*0000 0.0000 0.50000 1.4868 0.3121 -0.2079 0.9561 1.2596 2*6406 0.0000 0.0000 0.76923 1.4699 0.5412 0.3771 0.6583 1.5179 1.4695 0.0000 0.0000 0.95309 1.3948 0.7094 0.9349 0.2909 0.9173 0.3331 0.0000 0,0000 o 1 1 2 B02 01 02 02 R TD T 0*04691 0.0000 0*0000 0.0000 0.0000 0.604 0.000 0.000 0.23077 0.0000 0.0000 0.0000 0.0000 0.486 0.000 0,000 0.50000 0,0000 0.0000 0.0000 0.0000 0*376 0.000 0.000 0.76923 0.0000 0.0000 0.0000 0.0000 0.296 0.000 0.000 0.95309 0.0000 0.0000 0.0000 0.0000 0.256 0.000 0.000 154

Table B-l Comparison of Exact and Approximate Integrated Reflectances Isotropic and Rayleigh Scattering o0 Exact. iso* Six-flux, iso. Six-flux,Ray. Two-flux,iso. Richards,iso. and Ray. We= 0.90,1=oo: 0.00000 0.677 0.569 0.569 0.519 0.500 0.04691 0.654 0.560 0.560 0.519 0.489 0.23077 0.581 0.527 0.527 0.519 0.448 0.50000 0.507 0.484 0.485 0.519 0.400 0.76923 0.452 0.440 0.441 0.519 0.361 0.95309 0.422 0.389 0.391 0.519 0.338 1.00000 0.416 0.333 0.337 0.519 0.333 o.= 0.90,rl=l.:0: 0.00000 0.624 0.533 0.534 0.519 0.457 0.04691 0.598 0.520 0.521 0.519 0.446 0.23077 0.504 0.466 0.467 0.496 0.393 0.50000 0.394 0.375 0.376 0.421 0.316 0.76923 0.316 0.303 0.304 0.330 0.257 0.95309 0.279 0.259 0.260 0.293 0.227 1.00000 0.271 0.137 0.138 0.284 0.220.=0.90,r1= —O.50: 0.00000 0.584 0.507 0.509 0.519 0.428 0.04691 0.548 0.492 0.494 0.519 0.413 0.23077 0.417 0.390 0.391 0.420 0.333 0.50000 0.279 0.267 0.267 0.285 0.231 0.76923 0.209 0.199 0.199 0.211 0.171 0.95309 0.179 0.166 0.167 0.182 0.147 1.00000 0.170 0.079 0.079 0.175 0.141 c.=0.90,y1=0.25: 0.00000 0.571 0.488 0.491 0.519 0.405 0.04691 0.498 0.463 0.467 0.505 0.388 0.23077 0.300 0.283 0.283 0.302 0.243 0.50000 0.173 0.1.166 0.175 0.146 0.76923 0.125 0.117 0.117 0.124 0.103 0.95309 0.105 o.o096 o.o96 0.103 0.084 1.00000 0.101 0.043 0.043 0.099 0.081 =o=0.60,-1=0.50: 0.00000 0.345 0.235 0.234 0.225 0.176 0.04691 0.318 0.228 0.227 0.225 0.168 0.23077 0.232 0.190 0.193 0.211 0.135 0.50000 0.154 0.137 0.138 0.164 0.092 0.76923 0.114 0.103 0.104 0.129 o.o069 0.95309 0.097 0.082 0.084 0.112 0.059 1.00000 0.093 0.034 0.035 0.108 0.057.=O.30, 1=.50: 0.00000 0.156 0.091 0.091 o0.089 0.064 0.04691 0.142 0.088 0.088 0.089 0.061 0.23077 0.100 0.075 0.076 o.o89 0.048 0.50000 0.066 0.056 0.057 0.086 0.033 0.76923 0.048 0.042 0.043 0.079 0.025 0.95309 0.041 0.032 0.034 0.074 0.021 1.00000 0.040 0.012 0.014 0.072 0.020 * Extrapolated by Lagrange's formula for = o and 40.0. 155

Table B-2 Comparison of Exact and Approximate Total Transmissions Isotropic and Rayleigh Scattering Exact. iso.* Six-flux iso. Six-flux,Ray. Two-flux,iso. Richards,iso. o. =0.90,rl=1.0: 0.00000 0.183 0.181 0.173 0.000 0.211 0.04691 0.197 0.194 0.185 0.001 0.221 0.23077 0.271 0.265 0.256 0.189 0.269 0.50000 0.398 0.413 0.412 0.420 0.385 0.76923 o.498 0.523 0.523 0.549 0.492 0.95309 0.536 0.567 0.568 o.60g 0.549 1.00000 0.541 0.562 0.564 0.621 0.562 4. =O.90,T1= — 50: 0.00000 0.257 0.240 0.228 0.000 0.278 0.04691 0.286 0.256 0.244 0.025 0.292 0.23077 0.426 0.415 0.412 0.395 0.403 0.50000 0.607 0.611 0.610 0.621 0.582 0.76923 o.693 0.710 0.710 0.724 0.686 0.95309 0.728 0.749 0.749 0.767 0.733 1.00000 0.742 0.742 0.743 0.776 0.742 ==0.90,-1=O.25: 0.00000 0.277 0.277 0.264 0.000 0.322 o.04691 0.365 0.308 0.294 0.137 0.340 0.23077 0.607 0.597 0.594 o.600 0.574 0.50000 0.768 0.768 0.768 0.776 0.748 0.76923 0.827 0.837 0.837 0.844 0.823 0.95309 0.854 0.864 0.864 0.871 0.853 1.00000 0.861 0.859 0.859 0.877 0.859 w.=o.60, 1= —0.50: 0.00000 0.139 0.102 0.095 0.000 0.104 0.04691 0.154 0.109 0.101 0.000 0.108 0.23077 0.283 0.255 0.252 0.061 0.217 0.50000 0.502 0.491 0.491 0.269 0.454 0.76923 0.623 0.618 0.620 o.421 0.588 0.95309 0.675 0.668 0.670 o.496 0.648 1.00000 o.694 0.661 0.664 0.512 0.661'O=0.30,v1-0.50: 0.00000 0.057 0.038 0.035 0.000 0.039 0.04691 0.064 0.040 0.037 0.000 0.041' 0.23077 0.184 0.168 0.167 0.026 0.156 0.50000 0.425 0.417 0.418 0.186 0.399 0.76923 0.566 0.561 0.562 0.335 0.546 o.95309 0.628 0.622 0.622 0.411 0.612 1.00000 0.651 0.626 0.629 0.430 0.626 * Extrapolated by Iagrange's formula for 0o = o and o = 1.0. 156

Table B-3 Comparison of Exact and Approximate Integrated Reflectances Anisotropic Scattering, F = 0.9330, P = 0.4821 Exact* Six-flux Two-flux Exact* Six-flux Two-flux.= 0.90 o=0.90 71 = -~: ~1=1.0: 0.00000 0.628 0.467 0.195 0.567 0.433 0.195 0.04691 0.595 0.457 0.195 0.530 0.421 0.195 0.23077 0.489 0.425 0.195 0.396 0.362 0.143 0.50000 0.377 0.383 0.195 0.239 0.259 o0.0o89 0.76923 0.292 0.338 0.195 0.131 0.196 0.064 0.95309 0.244 0.286 0.195 0.079 0.161 0.054 1.00000 0.233 0.228 0.195 0.067 0.074 0.052 C. =.g90 = 0.50: W,= 0.90 1l= 0.25: 0.00000 0.547 0.417 0.195 0.552 0.405 0.195 0.04691 0.500 0.400 0.186 0.466 0.371 0.155 0.23077 0.333 0.277 o.o98 0.238 0.183 0.056 0.50000 0.163 0.170 0.052 0.097 0.099 0.028 0.76923 0.078 0.120 0.035 0.043 0.068 0.019 0.95309 0.041 0.098 0.029 0.018 0.055 0.015 1.00000 0.030 0.041 0.028 0.013 0.022 0.014. 0.60 o*= 0.60 1: 71= 0.50: 0.00000 0.336 0.165 0.046 0.330 0.162 0.046 0.04691 0.303 0.160 0.046 0.293 0.157 0.046 0.23077 0.206 0.142 0.046 0.180 0.119 0.039 0,50000 0.124 0.120 0.046 0.085 0.077 0.027 0.76923 0.073 0.097 0.046 0.040 0.055 0.020 0.95309 0.048 0.071 0.046 0.021 0.041 0.017 1.00000 0.043 0.041 o.046 o.016 0.012 o0.016 0, = 0~30 ~0,. 0.30 r1 = oo: t1= o.o: 0.00000 0.153 0.057 0.014 0.151 0.056 0,014 0.04691 0.133 0.054 0o014 0.131 0.054 0.014 0.23077 0.080 0.047 0.014 0.075 0.043 0.013 0.50000 0.042 0.039 0.014 0.034 0.029 0.011 0.76923 0.020 0.030 0.014 0.015 0.020 0.008 0.95309 0.011 0.020 0.014 0.008 0.014 0.007 1.00000 0.010 0.008 0.014 0.007 0.003 0.007 * Extrapolated by Lagrange's formula for io = o and 1.0. 157

Table BHA Comparison of Exact and Approximate Total Transmissions Anisotropic Scattering, F = 0.9330, P = 0.4821 11o Exact* Six-flux Two-flux Exact* Six-flux Two-flux. = 0.90. — =0.90 = 1.0: T1=0.50: 0.00000 0.236 0.220 0.000 0.291 0.269 0.000 0.04691 0.259 0.233 0.041 0.331 0.288 0.198 0.23077 0.372 0.339 0.511 0.507 0.515 0.711 0.50000 0.553 0.523 0.730 0.723 0.707 0.853 0.76923 0.692 0.632 0.814 0.833 0.790 0.902 0.95309 0.750 0.673 0.847 0.885 0.821 0.920 1.00000 0.761 0.662 0.853 o.go4 0.811 0.923 o. = 0.90'1 = 0.25: 0.00000 0.297 0.299 0.000 0.04691 0.396 0.354 0.439 0.23077 o.669 o.698 0.842 0.50000 0.853 0.835 0.923 0.76923 0.933 0.889 0.949 0.95309 0.979 0.905 0.959 1.00000 0.991 0.900 0.961 = 0.60 o — =0.30 = 0.50: T1 =0.50: 0.00000 0.143 0.100 0.000 0.041 0.033 0.000 0.04691 0.165 0.106 0.009 0.063 0.035 0.000 0.23077 0.319 0.303 0.386 0.198 0.175 0.210 0.50000 0.565 0.550 0.644 0.452 0.445 0.487 0.76923 0.701 0.671 0.751 0.601 0.587 0.626 0.95309 0.761 0.716 0.794 0.668 0.645 0.686 1.00000 0.783 0.710 0.803 0.692 0.652 0.698 * Extrapolated by Lagrange's formula for to = o and 1.0 158

Table B-5 Comparison of Exact and Approximate Integrated Reflectances and Total Transmissions Anisotropic Scattering, F1 = 0.750, 0w= 0.90, r1 = 1.0 Integrated Reflectances Total Transmissions 1o Exact* Six-flux Two-flux Exact* Six-flux Two-flux P = 0.5821: 0.00000 0.600 0.484 0.402 0.184 0.171 0.000 0.04691 0.563 0.471 0.402 0.203 0.183 0.006 0.23077 0.438 0.411 0.357 0.303 0.289 0.310 0.50000 0.303 0.315 0.261 0.472 0.467 0.560 0.76923 0.213 0.248 0.201 0.605 0.583 0.678 0.95309 0.173 0.209 0.173 0.657 0.629 0.727 1.00000 0.165 0.105 0.168 0.665 0.621 0.738 P = 0.4821: 0.00000 0.598 0.483 0.402 0.203 0.191 0.000 0.04691 0.564 0.470 0.402 0.222 0.203 o.oo6 0.23077 0.445 0.410 0.357 0.318 0.295 0.310 0.50000 0.310 0.314 0.261 0.476 0.469 0.560 0.76923 0.218 0.247 0.201 0.598 0.582 0.678 0.95309 0.175 0.207 0.173 0.645 0.626 0.727.ooo0000 0.166 0.103 0.168 0.652 0.618 0.738 P = 0.3889: 0o.ooooo0000 o.589 0.480 0.402 0.221 0.213 0.000 0.04691 0.559 0.467 0.402 0.238 0.225 o.oo6 0.23077 0.448 o.409 0.357 0.329 0.308 0.310 0.50000 0.315 0.314 0.261 0.480 0.471 0.560 0.76923 0.220 0.246 0.201 0.596 0.580 0.678 0.95309 0.174 0.205 0.173 0.646 0.622 0.727 1.00000 0.164 0.101 0.168 0.656 0.614 0.738 * Extrapolated by Lagrange's formula for,o = o and 1.0 159

APPENDIX C Computer Programs All computer programs were written in MAD (Michigan Algorithm Decoder) compiler language and were compiled and executed on the IBM-704. A description of the language is given in reference (10). Programs written in this compiler language may be compiled and run on either the IBM-704 or 709 and in the near future provision will be made for processing MAD programs on the IBM-7090. The programs are presented in the following pages with brief explanatory notes. 1. Solution of the Integral Equations a) Input Data One general program was written to solve any set of integral equations for an arbitrary value of N and m. In order to specify the problem and its method of solution the following information is required. (The MAD symbols are placed in parenthesis for use in reading the listing of the programs). 1) Run Number A five-character number for identification of output (RUN) N Number of terms in expansion of phase function (NTERM) p Number of points to be used in numerical integration (NPOINT) 160

Ir1 Optical thickness of dispersion (TAU) TOL Iteration will be stopped whenever the maximum absolute value of the change in any function is less than TOL (TOL) w>O, Albedo for single scattering (ALBEDO) 2) a0... aN Coefficients of an angular distribution function, A(O)...A(NTERMS) 3)....'N Set of values of the independent variable for use in numerical integration (See equation I-35 ), U(O)...U(NPOINT) 4) WO ***YW N Weight coefficients for numerical integration (See equation I-35 ), W(O)...W(NPOINT) 5 ) q m() Initial value of Psi Functions, PSI(O,L,I) i- O-., p A 0, *, N 6) +, (;)() Initial values of PHi Functions, PHI(O,L,I) i --- o...- p The preceeding six groups of information are referred to as the "basic data package." They completely define a solution to the integral equations or any trial solution. Thus, the computer starts with a basic data package as an initial approximation; it then performs a certain number of iterations; and, finally, produces as output a new basic data package. The basic data package is preceeded by a header card 161

giving the maximum number of iterations allowable and information as to which of the special options (described next) are desired. b) Special Options The following options are available in the routine: 1) Items 5 and 6 (i.e. the initial values of the functions) of the data above may be deleted in which case the following expressions will be used for the initial functions. (o) rn ha (Al,) e PA (l;) 2) The results of a computation may be punched on cards in the form of a basic data package for use in performing more iterations later or for use in another program if desired. 3) A reversal of direction of the change in /' or ~; may be interpreted as convergence and the iteration stopped when this occurs if desired. In addition the following features are available for speeding up the solution under certain special circumstances. a) if r', So80, calculation of the ~~4p.)will be deleted and F(4; ) will be set to zero b) if any angular distribution coefficient is zero (for example if ak = O) then calculation of the corresponding functions will not be calculated but set to zero 4; (,) 3 0 The number of equations in a set for use in estimating the computer time or cost in I-56 or I-58 is equal to the number of non-zero terms in the phase function. 162

c) Restrictions At the time of compilation of the MAD program an upper limit must be placed on the following so that the computer may properly allocate storage to the various arrays. 1) Maximum value of N, i.e. the maximum number of terms in the phase function (MAXL) 2) Maximum number of consecutive approximations which may be stored in the computer at one time (NMAX). Iteration will proceed computing f*(O), f*(l),....,f*(NM AX) after which f*(O) will be set to f*(NMAX) and the process repeated. 3) Maximum number of points to be used in numerical integration (MAXPT) Large values of NMAX result in a slightly more efficient program, however NMAX may be any value greater than zero. The total storage required including the approximately 5000 locations required by the program and its subroutines is given by STORAGE = 2(NMAX + 1)(MAXL + 1)(MAXPT + 1) + 4(MAXL +1)(MAXPT +1) + 4(MAXPT) + 5000 In practice, the large computer time required for large values of N and/or large values of p is a more serious restriction than the limitation on storage. An additional restriction is that all numbers expressed in floating point form must be of the following range. 37 -37 10 > IxI> 10 163

Thus, if Pmin is the smallest non-zero value of pi, then e must be greater than ]los in order to evaluate which is needed in performing some of the integrations. (Note: This does not apply if 71~ 80 as the 01 are not computed in this case.) The only restriction on the integration procedure is that F -7 P W= 1.o0 *o0, IE wi;r. o. t io i=0 o ) or in other words, the integrals ~ Jx21, ( xdx o~ must be done correctly. d) General Outline of the Programs The general outline of the programs is as follows: 1) A few initializing calculations are performed to give the computer information as to the maximum size of various arrays. 2) The header card is read followed by the basic data package. 3) Iterations are performed until convergence is indicated or until the allowable number of iterations is exceeded. 4) The results are punched and/or printed. A printed copy of the MAD Program follows. 164

$COMPILE MAD. PUNCH OBJECT, PRINT OBJECT 5THI0001 R R GENERALIZED PROGRAM FOR SOLUTION OF LIGHT SCATTERING R INTEGRAL EQUATIONS R BOOLEAN C1l C2. C49 C5. PCH, BL,BJBK,CONVGTHICK,OMIT,SIGN, 1 CONVG 1 INTEGER I,J,KL,MNNPOINTNTERMSMAXITITCNT,INT1,INT2, 1 FIXARG~RUN~NMAX1,TYPE~TIME. 0, I 112 DIMENSION BLOOP(100) VECTOR VALUES NMAX = 3 VECTOR VALUES MAXL = 3 VECTOR VALUES MAXPT = 20 R FOLLOWING DIMENSIONED (NMAX+1)*(MAXL+1)*(MAXPT+1) - 1 DIMENSION PSI(335,PSIDIM),PHI(335,PHIDIM) R FOLLOWING DIMENSIONED (MAXL+1)*(MAXPT+1) + 1 DIMENSION S(83.SDIM),POL(83,POL(DIM)CHPSI(83,CHPSID), CHPHI (83,CHPHID) R FOLLOWING DIMENSIONED MAXPT DIMENSION SS(20),W(20),U(20),BJ(20) R FOLLOWING DIMENSIONED MAXL DIMENSION PSIDEV(3),PHIDEV(3),OMEGA(3), A(3),BL(3),BK(3) 1 OMIT ( 3 ) DIMENSION PSIDIM(3) PHIDIM(3),SDIM(2) POLDIM(2),CHPSID(2), 1 CHPHID(2) VECTOR VALUES PSIDIM = 3 VECTOR VALUES PHIDIM = 3 VECTOR VALUES SDIM = 2 VECTOR VALUES POLDIM = 2 VECTOR VALUES CHPSID = 2 VECTOR VALUES CHPHID = 2 PSIDIM(3) = MAXPT + 1 PSDIMD(2) = MAXL + 1 PSIDIM(1) = PSIDIM(2)*PSIDIM(3) + PSIDIM(3) + 1 PHIDIM(3) = MAXPT + 1 PHIDIM(2) " MAXL + 1 PHIDIM(1) - PHIDIM(2)*PHIDIM(3) + PHIDIM(3) + 1 SDIM(2) = MAXPT + 1 SDIM(1) = SDIM(2) + 1 POLDIM(2) = MAXPT + 1 POLDIM(1) = POLDIM(2) + 1 CHPSID(2) = MAXPT + 1 CHPSID(I) = CHPSID(2) + 1 CHPHID(2) = MAXPT + 1 CHPHID(1) = CHPHID(2) + 1 START EXECUTE ZERO.(BL...BL(MAXL),BJ...BJ(MAXPT),BK...BK(MAXL)) READ FORMAT INBOOL, I,BL...BL(I) READ FORMAT INBOOL, I,BJ...BJ(I) READ FORMAT INBOOL, IoBK...BK(l) READ FORMAT IN1, MAXIT, TYPEPCH,SIGN READ FORMAT INA, RUNMNTERMSNPOINTTAUTOLtALBEDO PRINT FORMAT FORM1, RUN,M,NTERMSNPOINTMAXITTOL,TAUALBEDO THICK = OB WHENEVER TAU.GE. 80., THICK = 1B WHENEVER TYPE *E. 0 PRINT FORMAT FORM9, $READ IN$ OR WHENEVER TYPE.E. 1 PRINT FORMAT FORM9, $CALCULATED FROM LEGENDRE POLYNOMIALS$ END OF CONDITIONAL 165

READ FORMAT IN2, A...A(NTERMS) PRINT FORMAT FORM2 PRINT FORMAT FF1, A...A(NTERMS) EXECUTE ZERO. (OMEGA...OMEGA(NTERMS),OMIT...OMIT(NTERMS)) THROUGH ST6G FOR L = M,1,L.G.NTERMS OMEGA(L) = A ( L)*FACT. (L-M)/FACT. ( L+M)*ALBEDO ST6 WHENEVER OMEGA(L) *E 0O., OMIT(L) = 18 PRINT FORMAT FORM3 PRINT FORMAT FF1, OMEGA...OMEGA(NTERMS) READ FORMAT IN2, U...U(NPOINT) PRINT FORMAT FORM4 PRINT FORMAT FF1, U...U(NPOINT) READ FORMAT IN2, W..*W(NPOINT) PRINT FORMAT FORM5 TEST: -0.5 WEIGHT = -1.0 THROUGH SL6* FOR I = O,1,I.G.NPOINT TEST = TEST + W(I)*U(I) SL6 WEIGHT = WEIGHT + W(I) PRINT FORMAT FF1, W...W(NPOINT),WEIGHT,TEST WHENEVER.ABS.WEIGHT.G.OOQOOO0001 PRINT FORMAT BCD $1SSUM OF WEIGHTS NOT 1-0$ OR WHENEVER *ABS.TEST.G. 0000001 PRTNT FORMAT BCD, $1SUM OF W(I)*U(I) NOT o0.5 OR WHENEVER M-G.NTERMS PRINT FORMAT BCD, S1M IS GREATER THAN NTERMSS OTHERWISE TRANSFER TO AWAY END OF CONDITIONAL EXECUTE ERROR. AWAY CONTINUE THROUGH S1, FOR L = M,1TL.G.NTERMS THROUGH S1i FOR J = O,1,J.G.NPOINT S1 POL(L,J) - LEG.(U(J),ML) PRINT FORMAT FORM6 THROUGH ST3# FOR L =M,1r,L.GNTERMS ST3 PRINT FORMAT FF4,L,M,POL(LO)...POL(LNPOINT) WHENEVER TYPE.E. 0 THROUGH SL1, FOR L = M,1*L.G.NTERMS SL1 READ FORMAT IN2, PSI(OL.O)...PSI(OL,NPOINT)*PSIDEV(L) THROUGH SL2, FOR L = M,1,L.G.NTERMS SL2 READ FORMAT IN2, PHI(QO,LO)...PHI(O,LNPOINT),PHIDEV(L) OR WHENEVER TYPE.E. 1 THROUGH SL7A, FOR L = M,1,LG.NTERMS PSIDEV(L) = O0 PHIDEV(L) = 0. WHENEVER OMIT(L) EXECUTE ZERO.(PSI(0,L,O)...PSI(O,L,NPOINT)) EXECUTE ZERO.(PHI(O0LO)...PHI(O,LNPOINT)) TRANSFER TO SL7A END OF CONDITIONAL THROUGH SL7, FOR J = 0t,1J.G.NPOINT PSI(OLJ) = POL(LJ) WHENEVER U(J).G. 0..AND. TAU/U(J).L.80, PHI(O,LJ) = EXP.(-TAU/U(J))*POL(LJ) OTHERWISE PHI(O,LJ) = 0. SL7 END OF CONDITIONAL 166

SL7A CONTINUE OTHERWISE EXECUTE ERROR. END OF CONDITIONAL PRINT FORMAT FORM7 THROUGH ST4~ FOR L = M*,IgL.GNTERMS ST4 PRINT FORMAT FF2,L,M,PSI (OtL~0)...PSI(OLNPOINT) PRINT FORMAT FORM8 THROUGH ST5. FOR L = M,1*LG.NTERMS ST5 PRINT FORMAT FF3,L,M,PHI(O,LO).e.PHI(OLNPOINT) CONVG = OB ITCNT = 1 PRINT FORMAT BCD,I1BEGINNING OF ITERATIVE SOLUTIONS TIME1 = TIME.(O) LOOP THROUGH S6. FOR N = 1,ltNG*NMAX.OR.CONVG.OR.ITCNT*G.MAXIT EXECUTE ZERO.(PSIDEV(M)...PSIDEV(NTERMS)) THROUGH A1l FOR L = M,l1L.G.NTERMS WHENEVER OMIT(L) EXECUTE ZERO.(PSI(NLO)*..PSI(N~L~NPOINT)) TRANSFER TO Al END OF CONDITIONAL THROUGH A29 FOR J = O0,,J.G.NPOINT WHENEVER U(J) E., 0. PSI(N,L9J) = POL(L,J) TRANSFER TO A5 END OF CONDITIONAL EXECUTE SETCHK. SUMX = 0. THROUGH A3. FOR I = Ot1,IeG.NPOINT XXX = 0O THROUGH A4. FOR K = M,liK.G.NTERMS WHENEVER OMIT(K), TRANSFER TO A4 TS1 = (-1.0).P.(K+L)*F.(O) WHENEVER C1~S(KI) = TS1 XXX = XXX + TS1 A4 CONTINUE WHENEVER C2$SS(I) = XXX A3 SUMX = SUMX + W(I)*XXX*POL(LI) PSI(N,L,J) = POL(L,J) + O.5*U(J)*SUMX A5 TS1 = PSI(NLJ) - PS.I(N-1.LJ) EXECUTE CHECK. A2 WHENEVER.ABS.TS1.G.PSIDEV(L), PSIDEV(L) =.ABS. TS1 Al CONTINUE Re..........eee-e —-,e-e-eOOOeeO 004e.eo.* e 4* 00 ***-.e*..0.... EXECUTE ZERO.(PHIDEV(M).-.PHIDEV(NTERMS)) WHENEVER THICK THROUGH SL109 FOR L = M,1,L.G.NTERMS SLIO EXECUTE ZERO (PHI(N,L,O)..PHI(NL,NPOINT)) TRANSFER TO ON END OF CONDITIONAL THROUGH B1, FOR L = M,1,L.G.NTERMS WHENEVER OMIT(L) EXECUTE ZERO.(PHI(NL,O)...PHI(NLNPOINT)) TRANSFER TO BI END OF CONDITIONAL THROUGH B29 FOR J = O.i,J*G.NPOINT WHENEVER U(J).E. 0.

PHI(NLJ) = O0 TRANSFER TO B5 OR WHFNEVER TAU/U(J).G. 80. EXPON = O. OTHERWISE EXPON = EXP.(-TAU/U(J)) END OF CONDITIONAL EXECUTE SETCHK, SUMY = O. THROUGH B3, FOR I = O,1,I.G.NPOINT YYY = 0. THROUGH 84, FOR K = Mt1,K.G.NTERMS WHENEVER OMIT(K), TRANSFER TO B4 TS1 = G.(O) WHENEVER C4, S(KtI) = TS1 YYY = YYY + TS1 B4 CONTINUE WHENEVER C5, SS(I) = YYY 83 SUMY = SUMY + W(I)*YYY*POL(LI) PHI(NL,J) = EXPON*POL(LJ) + 0.5*U(J)*SUMY B5 TS1 = PHI(NL,J) - PHI(N-1,L,J) EXECUTE CHECK. B2 WHENEVER.ABS.TS1 *G. PHIDEV(L), PHIDEV(L) - oABS. TS1 B1 CONTINUE Re.... o o^eoo oo @0000000C le040. Oo4.00.0000004....00oo@0000000 ON PRINT FORMAT HEAD1.ITCNT PRINT FORMAT FF1, U...U(NPOINT) CONVG 1= 8 CONVGi = OB THROUGH S7A9 FOR L = M,1,L.G.NTERMS WHENEVER OMIT(L),TRANSFER TO S7A THROUGH S11 FOR J = O,1,J.G.NPOINT TS1 = PSI(N,LJ) - PSI(N-1,LJ) WHENEVER L.NE.O, TRANSFER TO 511 WHENEVER TS1*CHPSI(L.J).L.0. *AND. ITCNT.G,1 *AND. SIGN, 1 CONVG1 = 1B Sl CHPSI(L,J) = TS1 PRINT FORMAT FF2,LMCHPSI(L,O)...CHPSI(LNPOINT),PSIDEV(L) S7A WHENEVER PSIDEV(L).G.TOL, CONVG = OB THROUGH S7B, FOR L = M,1,L.G.NTERMS WHENEVER OMIT(L), TRANSFER TO S7B THROUGH S12, FOR J = 0,1,J.G.NPOINT TS1 = PHI(NL,J)-PHI (N-1, L,J) WHENEVER L.NE.0, TRANSFER TO S12 WHENEVER TS1*CHPHI(L,J).L.0. *AND. ITCNT.G.1 *AND. SIGN, 1 CONVG1 = 1B S12 CHPHI(L,J) T TS1 PRINT FORMAT FF3,LM,CHPHI(LO)...,CHPHI(LNPOINT),PHIDEV(L) S78 WHENEVER PHIDEV(L).G.TOL, CONVG = 08B CONVG = CONVG.OR. CONVG1 S6 ITCNT = ITCNT + 1 NMAX1 = N-1 WHENEVER CONVG PRINT FORMAT BCD,S4SUCCESSFUL SOLUTIONS OR WHENEVER ITCNT.G. MAXIT PRINT FORMAT BCDS4DESIRED ACCURACY COULD NOT BE OBTAINED 1 IN SPECIFIED NUMBER OF ITERATIONS$ OTHERWISE 168

THROUGH ALPHA5, FOR L = M,1,L.G.NTERMS THROUGH ALPHA5, FOR I = O1,1IGG.NPOINT PSI(O,L,I) = PSI(NMAX1,L,I) ALPHA5 PHI(O,L,I) = PHI(NMAX1,L,I) TRANSFER TO LOOP END OF CONDITIONAL TIME2 = TIME.(O) PRINT FORMAT FTIME (TIME2-TIME1)*3.6 PRINT FORMAT BCD, $OVALUES OF FUNCTIONS AFTER LAST ITERATION$ N = NMAX1 THROUGH S13, FOR L = M,1,L.G.NTERMS S13 PRINT FORMAT FF2,L,M,PSII(NLO)...PSI(N,L,NPOINT),PSIDEV(L) THROUGH S14, FOR L = M,1,L.G.NTERMS S14 PRINT FORMAT FF3,L,M,PHI(N,L,O)...PHI(N,L,NPOINT),PHIDEV(L) WHENEVER PCH PRINT FORMAT BCD, $ORESULTS PUNCHED$ INA(4) = RUN IN2(3) = RUN PUNCH FORMAT INA,RUN,MNTERMS,NPOINTTAU,TOL,ALBEDO PUNCH FORMAT IN29 A..oA(NTERMS) PUNCH FORMAT IN2, U..oU(NPOINT) PUNCH FORMAT IN2, W..oW(NPOINT) THROUGH SL3, FOR L = M,1tLG.NTERMS SL3 PUNCH FORMAT IN2, PSI(NMAX1,L,O)..*PSI(NMAXI,L,NPOINT), 1 PSyDEV(L) THROUGH SL4, FOR L: MqIL.G.NTERMS SL4 PUNCH FORMAT IN2, PHI(NMAX1,L,O)...PHI(NMAX1,L,NPOINT), 1 PHIDEV(L) END OF CONDITIONAL TRANSFER TO START R R INTERNAL FUNCTIONS R INTERNAL FUNCTION ENTRY TO F. FUNCTION RETURN OMEGA(K)*(PSI(N-1,KJ)*PSI(N-1,K,I) 1 - PHI (N-,K,J)*PHI (N-1,KI ) )/(U( I )+U(J) R ENTRY TO G. WHENEVER I.NE.J FUNCTION RETURN OMEGA(K)*(PHI(N-1,K,J)*PSI(N-1,K.I) 1 - PSI (N-1,K,J)*PHI (N-1,KI ) ) /(U(J)-U( I ) OTHERWISE WHENEVER I.E.O I1 0 1 I2 = 2 OR WHENEVER I.E.NPOINT I0 = NPOINT - 2 I1 = NPOINT - 1 I2 = NPOINT OTHERWISE 10 = I-1 I1 - I I2 = I+1 END OF CONDITIONAL X = U(I) XO = U( IO) 169

X1 = U(I1) X2 U (12) FO a pHI(N-1,KIO) F1 = pHI(N-1*K9I1) F2 = PHI(N-19KI2) PHPR a (F2-FO)/(X2-XO) + (2.*X-XO-X2)* 1 ((x2-X1)*FO - (X2-XO)*F1 + (X1-XO)*F2)/ 2 ( (X-XO)*(X2-XO)*(X2-X ) FO = PSI(N-lKIO) Fl = PSI(N-lKI1) F2 = PSI(N-1,KI2) PSPR m (F2-FO)/JX2-XO) + (2.*X-XO-X2)* 1 ((X2-Xl)*FO - (X2-XO)*F1 + (Xl-XO)*F2)/ 2 ((X1-XO)*(X2-XO)*(X2-Xl)) FUNCTION RETURN OMEGA(K) *(PSI(N-1, K,J)*PHPR 1 - PHI(N-1,K,J)*PSPR) END OF CONDITIONAL END OF FUNCTION R INTERNAL FUNCTION (FIXARG) INTEGER FIXARG ENTRY TO FACT. TS1 = 1. THROUGH ALPHA2, FOR INTl = 1,1INTl.G.FIXARG ALPHA2 TS1 = TS 1 TS*INT1 FUNCTION RETURN TS1 END OF FUNCTION R INTERNAL FUNCTION ENTRY TO CHECK. WHENEVER BJ(J)*AND.BL(L) PRINT FORMAT NGFO,S$J$,JU(J) THROUGH SL5S FOR K = M1,,K.G.NTERMS SL5 WHENEVER BK(K), PRINT FORMAT NGEO, $KSK. 1 S(KO)...S(K,NPOINT) PRINT FORMAT GEO, J,SS...SS(NPOINT),TS1, 1 PSI(N,L,J),PHI(N,L,J) END OF CONDITIONAL FUNCTION RETURN R ENTRY TO SETCHK. WHENEVER BL(L)*AND.BJ(J) C1 = lB C2 = 1B C4 = lB C5 = 1B OTHERWI SE C1 = OB C2 = OB C4 = OB C5 = OB END OF CONDITIONAL FUNCTION RETURN END OF FUNCTION R R FORMAT SPECIFICATIONS AND PRE-SET CONSTANTS R VECTOR VALUES INA = $C5,3I5$,S93F100$,$5,522,$,$ 5H$,S $, 170

1 $*S VECTOR VALUES IN1 = $415,5F10.5*$ VECTOR VALUES IN2 = S1P5E14$9S.7,52t,$$ 5H$,s $9$*$ VECTOR VALUES INBOOL - $I5,S5,60I1/(7011)*$ VECTOR VALUES HEADI = $37HOCHANGE IN FUNCTION DURING ITERATIO 1N I3*$ VECTOR VALUES GIO = $1H019,11110/(S10011110)*$ VECTOR VALUES GFO = $1HO,I5,S4,10FlI.6/(SlOlOF11.6)*$ VECTOR VALUES GEO = $1HO, IS4,P10E11.3/(S1O,1PlOEll.3)*$ VECTOR VALUES BCD = $20C6*$ VECTOR VALUES NGFO = $1H09C1,2H =,I3,3,,10F11.6/ 1 (S1O,1OF11.6)*$ VECTOR VALUES NGEO = $1HO.C1,2H =,I3,S3,1Pl0E11.3/ 1 (SiO,lPlO E11.3)*S VECTOR VALUES FORM1 = $40HlSOLUTION OF INTEGRAL EQUATIONS FOR 1 RUN C5//6H M = I5S,SS,18HNUMBER OF TERMS = I5,S5, 219HNUMBER OF POINTS = I5/32H MAXIMUM NUMBER OF ITERATIONS 3= I5,S5,33HDESIRED DEGREE OF CONVERGENCE = 1PE10.1/ 4 7H TAU = 1PE11.3,S5,32H ALBEDO FOR SINGLE SCATTERING z 5 1PE11.3*$ VECTOR VALUES FTIME = $19HOTIME REQUIRED WAS F6.1, 1 8H SECONDS *$ VECTOR VALUES FORM2 = $30HOPHASE FUNCTION —A.*.A(NTERMS) *$ VECTOR VALUES FORM3 = 1 $47HoMODIFIED PHASE FUNCTION —OMEGA...*OMEGA(NTERMS) *S VECTOR VALUES FORM4 = $28HOVALUES OF MU —U..*U(NPOINT) *S VECTOR VALUES FORM5 = $30HOWEIGHT FACTORS —W...W(NPOINT) *S VECTOR VALUES FORM6 = 1 $47HoLEGENDRE POLYNOMIALS —P(L,M) FOR U*-*U(NPOINT) *S VECTOR VALUES FORM7 = 1$47HOINITIAL FUNCTIONS —PSI(L*M) FOR U..*U(NPOINT) *$ VECTOR VALUES FORM8 = 1$47HOINITIAL FUNCTIONS —PHI(L9M) FOR U..*U(NPOINT) *$ VECTOR VALUES FORM9 = $19HOINITIAL FUNCTIONS 8C6*$ VECTOR VALUES FF1 = $lHO,591lP9El2.4/(SlO,1P9El2.4)*S VECTOR VALUES FF2 = S5HOPSI( I11 1Hi Ilt 2H) 1P 9E12.4/ 1 (Slo,1P 9E12.4)*$ VECTOR VALUES FF3 = $5HOPHI( II, 1H, II, 2H) iP 9E12.4/ 1 (Slo,1P 9E12.4)*$ VECTOR VALUES FF4: $5HO P( It 1HH Ii, 2H) 1P 9E12.4/ 1 (Slog1P 9E12.4)*$ END OF PROGRAM 171

2. Subroutine to Calculate Associated Legendre Polynomials A subroutine was written to calculate Pl(x) The program should be self explanatory. 172

*COMPILE MAD,PUNCH OBJECT LEG 001 EXTERNAL FUNCTION (X,M,N) INTEGER N,MK,FIXARGtINT1,TWOLOGo ENTRY TO LEG. WHENEVER X.Go1..OR. X *L.O. *OR. M.G. N PRINT FORMAT BCD,$4BAD ARGUMENT FOR LEGENDRE POLYNOMIAL SU 1BROUTINE$ PRINT FORMAT FORMAXMN ERROR RETURN END OF CONDITIONAL U: (1-X)/2e T = FACTe(N+M)/(FACTs(M)*FACT. (N —M)) SUM = T THROUGH S1* FOR K = 1,1.K-G.N-M T = -T*(N+M+K)*(N-M-K+1) /(K*(M+K)) WHENEVER TWOLOG.(T) + TWOLOG,(U).*L -115, 1 TRANSFER TO OUT T = T*U S SUM = SUM + T OUT WHENEVER M.E.O, FUNCTION RETURN SUM WHENEVER X.E.1., FUNCTION RETURN O0 V = SQRT.(1.-X*X)/2, FUNCTION RETURN SUM*V.P.M INTERNAL FUNCTION (FIXARG) ENTRY TO FACT. TS1 = 1. THROUGH BETA2, FOR INT1 = 1,1,INT1.G.FIXARG BETA2 TS1 = TS1*INT1 FUNCTION RETURN TS1 END OF FUNCTION VECTOR VALUES BCD = S20C6*S VECTOR VALUES FORMA = SlHO 1PE15.6,2115*S END OF FUNCTION *ASSEMBLE,PUNCH OBJECT 2LOG 001 ORG 0 PGM PZE SIZE PZE BCD 1TWOLOG PZE START REL ORG 0 START CAL 1,4 STA *+1 CLA ** SSP SUB TWO ARS 27 TRA 2,4 TWO OCT 200000000000 SIZE SYN * END 173

3. Program to Calculate Integrated Reflectance and Transmission The program to calculate the integrated reflectance and transmission follows the same general outline as the program to solve the integral equations. It accepts as input a basic data package. The results for R(0o) and T(0o) are calculated using A-4.2 and A-4.4. 174

*COMPILE MAD9 PUNCH OBJECT,PRINT OBJECT REF 001 R R INTEGRATED REFLECTANCE PROGRAM R INTEGER I,J,K,LM,N,NPOINT,NTERMS,INT1iINT2,FIXARGRUN BOOLEAN PCH VECTOR VALUES NMAX = 0 VECTOR VALUES MAXL = 6 VECTOR VALUES MAXPT = 50 R FOLLOWING DIMENSIONED (NMAX+1)*(MAXL+I)*(MAXPT+1)- 1 DIMENSION PSI(356,PSIDIM), PHI(356,PHIDIM) R FOLLOWING DIMENSIONED MAXPT DIMENSION W(50),U(50),T(50),R(50) R FOLLOWING DIMENSIONED MAXL DIMENSION OMEGA(6),A(6),PSIDEV(6),PHIDEV(6) DIMENSION PSIDIM(3),PHIDIM(3) VECTOR VALUES PSIDIM = 3 VECTOR VALUES PHIDIM = 3 PSIDIM(3) = MAXPT + 1 PSIDIM(2) = MAXL + 1 PSIDIM(1) = PSIDIM(2)*PSIDIM(3) + PSIDIM(3) + 1 PHIDIM(3) = MAXPT + 1 PHIDIM(2) = MAXL + 1 PHIDIM(1) = PHIDIM(2)*PHIDIM(3) + PHIDIM(3) + 1 START READ FORMAT IN1, PCH READ FORMAT IN1, RUN,M,NTERMS,NPOINTTAU,TOLgALBEDO WHENEVER M.NE. 0 PRINT FORMAT BCD, $1M NOT ZEROS EXECUTE ERROR, END OF CONDITIONAL READ FORMAT IN2, A...A(NTERMS) EXECUTE ZERO.(OMEGA...OMEGA(NTERMS)) THROUGH ST6, FOR L = M,1,L.G.NTERMS ST6 OMEGA(L) = A(L)*FACT.(L-M)/FACT.(L+M)*ALBEDO READ FORMAT IN2, U...U(NPOINT) READ FORMAT IN2, W.*.W(NPOINT) THROUGH SL1, FOR L = M,1,Lo.GNTERMS SL1 READ FORMAT IN2, PSI(Q,L,0).*.PSI(OLNPOINT),PSIDEV(L) THROUGH SL2, FOR L = M,1*L.GNTERMS SL2 READ FORMAT IN2, PHI(0,L,0)..*PHI(0,LNPOINT),PHIDEV(L) PRINT FORMAT BCD, $lINPUT DATA$ PRINT FORMAT GIO, NPOINT, M, NTERMS. RUN PRINT FORMAT GFO, NTERMS, OMEGA...OMEGA(NTERMS) PRINT FORMAT GFO, O, U...U(NPOINT) PRINT FORMAT GFO,O,W.. W(NPOINT) PRINT FORMAT GEOOTAUTOL THROUGH ST3, FOR L = M,1,L.GNTERMS ST3 PRINT FORMAT GEO,LPSI(OL,O) 0..PSI(O,LNPOINT),PSIDEV(L) THROUGH ST4, FOR L = M.1,L.G.NTERMS ST4 PRINT FORMAT GEO,L,I-**(0,L,Q)...PHIOL(OLNPOINT)gPHIDEV(L) THROUGH A2, FOR J = O,1,J.G.NPOINT SUMX = O. SUMY = 0O THROUGH A3, FOR I = O,1,I*G*NPOINT XXX = O. YYY = 0. THROUGH A4, FOR K = M,1,K.G*NTERMS XXX = XXX + (-1.0).P.K * F.(O) A4 YYY = YYY + G.(O) 175

SUMX = SUMX + W(I)*U (I)*XXX A3 SUMY = SUMY + W(I)*U(I)*YYY R(J) = O.5*SUMX A2 T(J) = 05*SUMY PRINT FORMAT BCD, SOREFLECTION AND TRANSMISSION COEFFICIENTS$ PRINT FORMAT GEO,1,R...R(NPOINT) PRINT FORMAT GEO,2,T...T(NPOINT) WHENEVER PCH PUNCH FORMAT IN1, RUN, M, NTERMS, NPOINT, TAU, TOL, ALBEDO PUNCH FORMAT IN2, R...R(NPOINT) PUNCH FORMAT IN2, T...T(NPOINT) END OF CONDITIONAL TRANSFER TO START R R INTERNAL FUNCTIONS R INTERNAL FUNCTION INTEGFR IOI1,I2 ENTRY TO F. FUNCTION RETURN OMEGA(K)*(PSI(N-1,K,J)*PSI(N-1,KI) 1 _HI - PHI (N-1,KJ)*PHI (N-1,K, I) )/(U( I )+U(J)) R ENTRY TO G. WHENEVER I.NE.J FUNCTION RETURN OMEGA(K)*(PHI(N-l,KJ)*PSI(N-1,K,I) 1 - PSI (N-1,K,J)*PHI (N-1,K,I ) ) / (U(J)-U( I ) OTHERWISE WHENEVER I*E.O 0 = O0 I1 = 1 I2 = 2 OR WHENEVER I.E.NPOINT IO = NPOINT - 2 I1 = NPOINT - 1 I2 = NPOINT OTHERWISE IO 1 I-1 I1 = I 12 = I+1 END OF CONDITIONAL X = U(I) XO u(IO0) X1 = U( 11) Xl - U(I1) X2 = U(12) FO = PHI(N-1,KIO) F1 = PHI(N-1,K,I1) F2 = PHI(N-1,K#I2) PHPR = (F2-FO)/(X2-XO) + (2.*X-XO-X2)* 1 ((X2-X1)*FO - (X2-XO)*Fl + (X1-XO)*F2)/ 2 ((Xl-XO)*(X2-XO)*(X2-X1)) FO = pSI(N-1,K,IO) F1 = PSI(N-1,K,I1) F2 = PSI(N-1,K,I2) PSPR = (F2-FO)/(X2-XO) + (2.*X-XO-X2)* 1 ((X2-X1)*FO - (X2-XO)*F1 + (Xl-XO)*F2)/ 2 ((X1-XO)*(X2-X)*(X2-X1) ) FUNCTION RETURN OMEGA(K)*(PSI(N-1,K,J)*PHPR 1 - pHI(N-1,K,J)*PSPR) 176

END OF CONDITIONAL END OF FUNCTION R INTERNAL FUNCTION (FIXARG) INTEGER FIXARG ENTRY TO FACT. TS1 = 1. THROUGH ALPHA2g FOR INT1 = 1,1,INT1o.GFIXARG ALPHA2 TS1 = TSl*INT1 FUNCTION RETURN TS1 END OF FUNCTION R R FORMAT SPECIFICATIONS AND PRE-SET CONSTANTS R VECTOR VALUES IN1 = $415,5F10.5*S VECTOR VALUES IN2 = S1P5E14.7,S2*S VECTOR VALUES GIO = SHOI9,11IlO/(SlO911I10)*$ VECTOR VALUES GFO = S1HOI5tS4,lOFl.6/(SlO,1OF11.6)*S VECTOR VALUES GEO = SlHOI5,S4,1lPlOEll3/(SlOlPIOE11.3)*$ VECTOR VALUES BCD = $20C6*$ VECTOR VALUES N = 1 END OF PROGRAM 177

REFERENCES 1. S. Chandrasekhar, Radiative Transfer, Oxford, Clarendon Press, 1950. 2. C. M. Chu and S. W. Churchill, "Multiple Scattering by Randomly Distributed Obstacles - Methods of Solution", IRE Trans. on Antennas and Propagation, AP-4,142 (1956). 3. S. W. Churchill, J. H. Chin, G. C. Clark, B. K. Larkin and J. A. Leacock, "The Transmission of Thermal Radiation through Real Atmospheres" —AFSWP-1035, Contract No. Nonr-1224(17), Project No. NRo87-063, The University of Michigan (April 1957). 4. Symposium on "Atmospheric Transmission of Thermal Radiation" held at AFSWP Headquarters, Washington, D. C. (April 16, 1957). 5. S. W. Churchill, C. M. Chu, J. A. Leacock, and J. Chen, "The Effect of Anisotropic Scattering on Radiative Transfer", DASA-1184, Contract No. Nonr-1224(17), Project No. NR o87-o63, The University of Michigan (March 1960). 6. Internuclear Corporation - report to be issued (1961). 7. A. N. Lowan, N. Davids, and A. Levenson, "Tables of Zeros of the Legendre Polynomials of Order 1-16 and the Weight Coefficients for Gauss' Mechanical Quadrature Formula", Bull. Amer. Math. Soc. 48, 739 (1942). 8. S. Chandrasekhar, D. Elbert, and A. Franklin, "The X- and Y-Functions for Isotropic Scattering", Ap. J., 115, 244 (1952). 9. C. M. Chu and S. W. Churchill, "Representation of the Angular Distribution of Radiation Scattered by a Spherical Particle", J. Opt. Soc. Am. 45, 958 (1955). 10. B. Arden, B. Galler, and R. Graham, "Manual for the Michigan Algorithm Decoder" The University of Michigan, June 1960. (Available from The University of Michigan Computing Center, North University Building, Ann Arbor, Michigan). 11. C. M. Sliepcevich, S. W. Churchill, G. C. Clark and C. M. Chu, "Attenuation of Thermal Radiation by a Dispersion of Oil Particles". AFSWP 749, Part II, Contract No. DA-18-108-CML-4695, ORD CP3-601, Project No. 4-12-01-005, The University of Michigan (May 1954). 178

12. C. M. Chu and S. W. Churchill, "Numerical Solution of Problems in Multiple Scattering of Electromagnetic Radiation", J. Phys. Chem. 59, 855 (1955). 13. P. I. Richards, "Multiple Isotropic Scattering", Physical Review, oo, 517 (1955). 14. B. Davison, "Neutron Transport Theory", Oxford Press, 1955. 15. C. C. Grosjean, "A High Accuracy Approximation for Solving Multiple Scattering Problems in Infinitive Homogeneous Media", Il Nuovo Cimento, 3,1261, 1955. 16. S. Stein, D. E. Johnson and A. W. Starr, "Theory of Antenna Performance in Scatter-Type Reception", Hermes Electronics Co. Report, AF CRC-TR-59-191 September, 1959. 17. Devenport and W. Root, "Random Processes in Automatic Control", McGraw Hill (1956). 18. R. Miltra, "On the Solution of Eigen value Equation of the Weiner-Hopf Type in Finite and Infinite Ranges", I.R.E. Convention Record, Part 4, 170, (1959). 19. Morse and Fishbach,'Methods of Theoretical Physics" McGraw Hill (1953), Chapter 8. 179

DISTRIBUTION LIST ADDRESSEE ARMY NO OF CYS Chief of Research and Development, DA, Washington 25, D. C., 1 ATTN: Atomics Division The Quartermaster General, DA, Wash 25, D. C., ATTN: 1 R and D Division Commanding General, Army Medical Service School, Brooke Army 1 Medical Center, Ft. Sam Houston, Texas Commanding General, U. S. Army Chemical Corps, Research and 2 Development Command, Washington 25, D. C. Director, National Aeronautics and Space Administration, 2 1520 H St. N. W. Washington 25, D. C. Assistant Chief of Staff Intelligence, DA, Washington 25, D. C. 1 Commanding General, Aberdeen Proving Ground, Aberdeen, Md., 2 ATTN: Director, BRL Commanding Officer, Engineer Research and Development Lab., 1 Ft. Belvoir, Virginia, ATTN: Chief Tech Support Branch Commanding Officer, Picatinny Arsenal, Dover, New Jersey, 1 ATTN: ORDBB-TK Commanding Officer, Chemical Warfare Lab., Army Chemical Center, Md., 1 ATTN: Tech Library Commanding Officer, U. S. Army Signal Research and Development Lab., 1 Ft. Monmouth, New Jersey, ATTN: Technical Documents, Evans Area Director, Operations Research Office, Johns Hopkins University, 6935 1 Arlington Rd. Bethesda 14, Md. NAVY Chief of Naval Operations, DN, Washington 25, D. C., ATTN: OP-75 1 and OP-03EG Chief, Bureau of Naval Weapons, DN, Washington 25, D. C. 3 Chief, Bureau of Ships, DN, Washington 25, D. C., ATTN: Code 423 2 180

ADDRESSEE NAVY NO OF CYS Chief, Bureau of Yards and Docks, DN, Washington 25, D. C. 1 ATTN: D-440 Chief of Naval Research, DN, Washington 25, D. C., ATTN: Code 811 1 Superintendent, U. S. Naval Postgraduate School, Monterey, California 1 Commander, U. S. Naval Ordnance Lab., White Oak, Silver Spring, Md., 1 ATTN: EE Division and R Division Director, U. S. Naval Research Lab., Wash 25, D. C., ATTN: Code 2029 1 Commander, New York Naval Shipyard, Brooklyn 1, New York, ATTN: 1 Director, The Material Lab. Commanding Officer and Director, U. S. Naval Radiological Defense 3 Lab., San Francisco 24, Calif., ATTN: Tech Information Division Commanding Officer, U. S. Naval Development Center, Johnsville, Pa. 1 AIR FORCE HQ USAF (AFTAC) Wash 25, D. C. 1 Director of Research and Development, DCS/D, Hq USAF, Wash 25, D. C., 1 ATTN: Guidance and Weapons Division Air Force Intelligence Center, HQ USAF, ACS/I (AFCIN-3V1) Wash 25, 2 D. C. Commander, Air Research and Development Command, Andrews AFB, Wash 25, 3 D. C., ATTN: RDRWA Director, Air University Library, Maxwell AFB, Alabama 1 Commandant, School of Aviation Medicine, USAF Aerospace Medical 1 Center (ATC), Brooks AFB, Texas, ATTN: Col Gerritt L. Hekhuis Commander, Wright Air Development Center, Wright-Patterson AFB, Ohio, 1 ATTN: WCOSI AFCCDD, Laurence G. Hanscom Fld, Bedford Mass, ATTN: CRQST-2 1 AFSWC, Kirtland AFB N lex, ATTN: Tech Information Office 3 181

ADDRESSEE OTHERS NO OF CYS Director of Defense and Engineering, Washington 25, D. C., 1 ATTN: Tech Library Commander, Field Command, Sandia Base, Albuquerque, New Mexico 6 Chief, Defense Atomic Support Agency, Washington 25, D. C. 5 Chief, Defense Atomic Support Agency, Washington 25, D. C., 2 ATTN: Major Vickery Commander, ASTIA, Arlington Hall Station, Arlington 12, Virginia, 30 ATTN: TIPDR Los Alamos Scientific Laboratory, P. O. Box 1663, Los Alamos, New 1 Mexico, ATTN: Report Librarian Director, Lincoln Laboratory, Massachusetts Institute of Technology, 1 P. 0. Box 73, Lexington 73, Mass., ATTN: Publications (For Prof. G. C. Williams) Los Alamos Scientific Laboratory, P. O. Box 1663, Los Alamos, New 1 Mexico, ATTN: Report Librarian (For Dr.Alvin C. Graves) University of California, 405 Hilgard Street, Los Angeles 24, 1 California University of California, Lawrence Radiation Laboratory, P. 0. 1 Box 808, Livermore, California University of California, Lawrence Radiation Laboratory,,Berkeley 4, 1 California University of Rochester, River Campus, River Road, Rochester, New 1 York, ATTN: Dr. H. Stewart Director, Oak Ridge National Laboratory, P. 0. Box X, Oak Ridge, 1 Tennessee Director, National Bureau of Standards, Washington 25, D. C., ATTN: 1 Library 182

ADDRESSEE OTHERS NO OF CYS U. S. Atomic Energy Commission, Germantown, Maryland, ATTN: Library 1 Director, The RAND Corporation, 1700 Main Street, Santa Monica, 1 California, ATTN: Dr. Diermingian Director, Scripts Institution of Oceanography, 3602 La Jolla 1 Shore Drive, La Jolla, California 183

UNIVERSITY OF MICHIGAN III III9151lill ll02827ilii I4II 3 9015 02827 4929