016234-1-F THE UNIVERSITY OF MICHIGAN COLLEGE OF ENGINEERING DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING Radiation Laboratory SCATTERING AND ABSORPTION OF ELECTROMAGNETIC RADIATION BY CLOUDS OF DISC SHAPED AEROSOLS 16234-1-F = RL-2290 Final Report on Grant No. 04-78-B01-1 By Herschel Weil and C. M. Chu December 1978 for National Oceanic and Atmospheric Administration U.S. Department of Commerce Boulder, CO 80302 Ann Arbor, Michigan

Scattering and Absorption of Electromagnetic Radiation by Clouds of Disc Shaped Aerosols H. Weil and C. M. Chu

1. Introduction The purpose of this research has been to develop theoretical and numerical procedures which can contribute to the interpretation and prediction of lidar experimental observations of the scattering and absorption of polarized light by snow and ice crystals in cloud formations. Our results can also be used for the study of absorption and scattering of natural unpolarized or partially polarized radiation by any disc-like particles of lossy dielectric material. They can serve as inputs to studies of the influence on cloud formation of radiation absorption, of weather and climate as it is influenced by the greenhouse effect, and of various effects of air pollution. The scattering and absorption of electromagnetic radiation has been treated purely analytically only for spheres and infinite circular cylinders. For these cases the well known Mie theory is used. For all other cases recourse must be made to numerical (or experimental) methods or to approximate theories such as Rayleigh-Gans or geometric optics which are valid only in the extremes of very long or very short wavelength respectively relative to the body dimensions. This report describes a numerical technique we have developed to efficiently compute the scattering and absorption by thin circular discs of lossy dielectric (complex refractive index material). Sample numerical results are given for the case of infrared radiation incident on a disc of ice, the wavelengths being comparable to and somewhat less than the disc radius and about 10 times the disc thickness. Such a disc could be thought of as simulating a thin flat plate hexagonal ice crystal, particularly when the crystal is symmetrically rimed. Such crystals fall flat side horizontal in natural clouds (Zikmunda and Vali (1972)) and since their orientation in the horizontal plane is random a circle represents an ensemble average shape as well as approximating the individual crystal. The work we have done on this grant extends earlier work described in Chu and Weil (1976) and Weil and Chu (1976). In those papers a method,

-2 - tailored for its efficient applicability to thin discs, was developed to compute the current density distribution induced by an incident electromagnetic wave of arbitrary polarization and from this to compute the scattered fields and absorbed power. While the method could in principle apply for any directions of incidence relative to the discs, the formulas were in fact worked out and programmed only for incident propagation directions perpendicular or parallel to the flat surfaces of the discs. This of course severely limited the practical situations where the results could be applied. On the present grant, the method was pursued to generate the detailed formulas and associated computing programs for arbitrary directions of incidence. The new program, of course, automatically handles the broadside and edge-on cases when these are called for. It does this in an improved fashion based on some reformulation of the equations used previously. The present program, incorporating as it does all angles of incidence can now be used with appropriate number densities and probability distributions for crystal dimensions to model the effects of clouds of disc scatterers. However, the only experimental data on size distributions for naturally occurring plate crystals is that summarized in Auer and Veal (1970) which leads to the empirical relation between thickness 6 and radius a: 6 = 2.02 (2a)0 49 m. Multiple scattering and absorption (multiple photon-particle interactions) may be treated by radiation transfer computations. Some radiation transfer codes depend intrinsically on the use of Mie theory and so are limited to exploring the effects of spherical particles. The program we have developed for disc particles can however be used as an external subprogram or function input to those radiation codes which are not intrinsically limited to use of Mie theory.

-3 - 2. The Integral Equation and Reduction to Algebraic Equations The integral equation for current density J induced in a homogeneous dielectric object of refractive index n is i -1 J E = --- [L-J - ] (1) - n -- 1 where L is the self-adjoint integral operator defined by L-A = -j- dv'J(R')G(R,R') D -ff ds'J(R').n' (R')VG(R,R'), (2) S D and S represent the volume and surface of the object and G = exp(j IR - R' I)/(4T'IR - R' I). (3) All lengths are normalized by multiplication by k0 = 27T/X where A is the free space wavelength. The time dependence of the fields is assumed to be of the form ej t so that the refractive index n has a negative imaginary part except for ideal nondissipative, nonconductive dielectrics. Eq. (1) is reduced to a set of linear algebraic equations via the method of moments. To do so J(R) is expanded in a series of basis functions W.(R). J(R) = -jW ajWj (R).(4) The basis functions should satisfy the requirement V-'i = 0 in D W-n # 0 on S (5) to ensure the similar properties (which exist in a homogeneous medium) for J(R) itself. Substitution of (4) into (1) yields i- (6) E o= cti L ((W. - 1))W]. (6) i: - Using the inner product

-4 - < > = fdv' f(R')^g(R') (7) with (4) then yields an algebraic matrix equation to solve for the ai <W i> = i i[Zij - W /(n2 - 1)] (8) i ij where Z = <W. L W>, W.. <W,Wj>. (9) ij i'-i j 'ij -.] ( ' 2 The symmetric matrix [Z.. - W./(n - 1)] is termed an impedance 1j 1j 2 -1 matrix. All material effects are confined to the factor (n - 1)1 while shape effects are included in Z.. and W... To specify the basis functions W. we use a cylindrical coordinate system (p,{,z) with z along the symmetry axis of the disc. The disc fills a region given by -6/2 < z < 6/2, 0< p < a. Then, bearing (5) in mind, we choose the Wi from the following sets of functions. To most i concisely express the induced current when E1 is perpendicular to z we use the TE mode W(nm) Vx{zR(p)Z (z)o (p)} (10) ^E n m where Z (z) = z; m (<) = sin m) or 0 (%) = -cosmr. (11) n m m The form of Z (z) is chosen to permit easy integration over z when evaluating the impedance expressions. In accord with the thin disc 2 assumption only terms of order 6 or greater will be kept in the impedance, hence only n = 0 and 1 functions will be needed. The following odd TE modes resulting from the choice of sin mp in (13) result TEO ODD: W 0,m) = p(mR(p)/p)sin mq + {R' (p)cos m{ (1, m) TE1 ODD: W = p(mR(p)/p)z sin mc + 4R'(p)z cos m4, (12) etc. Similarly, TM modes defined by

-5 -W(n,m) V x (n,m) (13) most concisely express the currents induced when the incident magnetic field H. l z. The odd TM nodes are -- TMO ODD: W0) = z[(R' (p)/p) + R (p) - (m2R(p)/p2)] sinm4 TM1 ODD: m) = PR'(p)sinm( + (mR(p)/p)cosm4 + z[-(R' (p)/p) - R' (p) + (m2R(p)/p )]sinmn. (14) etc. For even modes m4 is to be replaced by m% + (7rr/2) in (11) and (14). When neither E nor H are perpendicular to 2, a combination of TE and TM mode W functions will be required. With these W functions there is no coupling between the different values of m; that is, the impedance terms corresponding to use of W(n,m) and W(nm') are zero when m ~ m'. Likewise there is no coupling between even and odd modes. We achieve further decoupling by using our assumption that the discs are thin and dropping all terms of order 63 or higher relative to the lower order terms when carrying out the volume integrations called for in determining the matrix elements. As a result the impedance matrix has many zero elements and can be partitioned as shown in Eqs. 3.3 - 3.6 of section 3 which greatly facilitates the numerical solution. In our earlier work involving solely broadside and edge-on incidence only TEO modes were excited and hence only matrix terms involving these modes were not zero [the z(TEO,TEO) terms of the present section 3]. We choose a set of radial functions R(p): {Ri(p)} so as to facilitate the integration. Each R. (p) is chosen to be zero over all but a fractional part of the interval 0 < p < a. The subintervals are taken small enough to permit the integrals to be evaluated by approximate means with adequate accuracy, yet an excessively large number of intervals is avoided with the choice of R.(p) given below. We are then able

-6 - to treat the singularity region to desired accuracy by reducing the problem to the singularity integration for a statics problem plus correction terms, all of which are evaluated analytically. We do end up, however, with rather lengthy expressions for the matrix elements. One needs second derivatives of R(p) so we choose symmetric functions made of segments of second order polynomials. Specifically, we divide the interval 0 < p < a into N contiguous segments of length a/N. Then the functions R.(p) are the spline functions illustrated schematically in Fig. 1 and defined for i = 1,2,...,N by (15). (a/N) (i-l)<p<(a/N)i, (1/2) (N/a) [p-(a/N) (i-l)]2 (a/N) i<p< (a/N) (i+l), -1+2 (N/a) [p-(a/N) (i-l)] -(1/2) (N/a) [p-(i-l) (a/N)] (15) (a/N) (i+l)<p<(a/N) (i+2) -1+2(a/N) [(a/N) (i+3)-p] -(1/2) (N/a) [(a/N) (i+3)-p]2 (a/N) (i+2)<p<(a/N) (i+3), (1/2) (N/a) [(a/N) (i+3)-p] provided p<a. For p>a, Ri(p)=0. This treatment of the radially dependent part of the problem is typical of techniques used in finite element methods so our numerical method may be thought of as a hybrid momentfinite element method. Note that there is an overlapping of the R!s in the integrations for Z.. and W... Ri (p) Ri-l R R1*+ Ri+2.5 i-I i i*l i*2 i3 P radial distance in units of a N Fig. 1. Illustration of the function defined in Eq. (15).

-7 - Numerical evaluation of the integrals is facilitated by use of average values over short intervals for R (p), and functions involving Ri(p), and its derivatives. For this purpose each interval (a/N)(i-l) < p < (a/N)(i+3) is subdivided into eight equal parts denoted by p=l,...,8 each extending from aip to.ip and centered on yi (see Eqs. 2.5 - 2.7 of section 3). We define averages over these intervals Bip 16) Aip = <R(p)> = (a/2N) f Ri(p)dp (16) ip and similarly B.p C = <Rp) C = )> ip 1 p ip 1 p Dip = <R(p)> E. = <R.(p)/p> - (17) 1p l p 1p 1 p (17) Table II in section 3 lists the explicit formulas for Ap... Ep. ip p The W.. and excitation terms <.j-E > each involve only a single volume integral and can be evaluated analytically with our choice of polynomial functions for the R.(p) by straightforward integration. The Zij each involve double volume and surface integrals with singular integrands. They are carried out by partly numerical and partly analytical methods in which the two z integrations are done analytically 2 to 0(6 ), the singularities are handled analytically and the two angular integrations reduced to a single integral approximated by a finite sum. The result of the analytic treatment of the singularities plus the use of the averages defined above reduces the radial integrals to finite sums. The results of this somewhat lengthy analysis leads to the formulas which are given in section 3 and where were programmed and used for computation. Once the a. and hence J are determined, classical EM theory permits one to compute the reradiated and the absorbed power by straightforward integrations over the disc volume which can be carried out analytically.

-8 - Specifically, the far-zone electric field in the direction R (s,s ) may be expressed as -jR-R -Far 4R - R4 s where N= kk3 // dv' ej s J(RR) (19) Then by substituting Eq. 4 for J one can express N as -2 N = -jwGok2 kPk(s,) (20) where p* ^ s*.The I and Is functions are given by Eqs. (4.5) of Section 3 for the present choice of basis functions Wk. They are used also in evaluating the excitations Wj,E of Eq. (8). The incident fields are taken to be of unit magnitude with elliptical polarization. They are expressed as E= (p + s)et (22) where p and s are complex and satisfy II2 + ql2 = 1. Then Eqs. (23) - (27) which follow and define the various cross-sections which characterize the scattering and absorption are evaluated. The final formulas which result are given by Eqs. (5.1) - (5.7) of Section 3.

-9 - a(O, = radiated power per solid angle in direction (G,2) 3) B incident power per area. the absorption cross-section = power dissipated in the object (24) A incident power per area the total scattering cross-section, TS = f B(8,)d (25) 47r the total cross-section, also called extinction cross section = aTS A (26) and the assymetry factor aASY = faB(')cos sdQ (27) where E is measured relative to the incident direction. aB and aA are given in terms of the a. by Eqs. 5.1 and 5.2 of section 3. 3. Computational Formulas The explicit formulas to be programmed are given in this section. 3.1 Ordering the Coefficients and Partitioning the Matrix The procedures involve solving the linear algebraic matrix equation (2.8) which we rewrite as z a = I (1.1) for unknown vector a, then computing several functions of the elements ai of a. The matrix z may be decomposed in the form z = Z - W/(n2 - 1) (1.2) where the elements of Z and W depend explicitly on two real parameters a 2 and 6 while n is a given complex number. I, the excitation vector depends on a and 6 plus four other real parameters 90, p0, p and s. Each

-10 - row of (1.1) (corresponding to the element Ik of I) is associated with a number of parameters; integers i, j=l...N, m=l...M and "mode labels" TMO, TEO, TM1, even and odd. The association with index k is given by Table I which covers the case m = 1. The k numbering continues from k = 6N+1 through 12N following the same ordering as the m=l case, but with m=2. This pattern continues until m=M. The maximum k value is K=6NM. Table I Index numbers, k, for the excitations Ik and coefficients ak sub-. matrix mode:see Eq. k label m i or j (1.4) 1 TM0 ODD 1 1 I N '. N N+l TEO ODD 1. - 2N. N 2N+1 TM1 ODD. 1 II 3N N 3N+l TMO EVEN 1 I 4N N 4N+l TEO EVEN. 1 5N N 5N+1 TM1 EVEN ' 1 TI 6N. N I..

-11 - Although the overall array of 6NM equations may be large, the solution is reduced analytically to a number of subproblems involving much smaller arrays as follows. First z may be shown to reduce to a diagonal array of submatrices so that equation subsets correspondong to "odd" and "even" and to the different m values are uncoupled and may be solved separately. Thus z(m=l, ODD) 0 0 - - - z= 0 z(m=l, EVEN) 0 0 0 z(m=2, ODD) - - - (1.3) Furthermore, although the excitations are different for each m in the odd and even cases, the z's are not; z(m=mo, ODD) = z(m=m EVEN) = z(m ). Each of the submatrices z(m) may be further partitioned as follows: z(m) = 0 II(m) (1.4) where the z(m) and z(m) are associated with mode combinations as indiI II cated by the following notation (and by the fifth column of Table I). z z(TMO, TMO) (1.5) z(TEO, TEO) z(TEO, TM1) zI z(TMl, TEO) z(TMl, TM1) (1.6) -~ In (1.5) and (1.6) the dependance on m is "to be understood". Each of the matrices on the right hand sides is specified by symmetric N x N arrays of complex elements. Thus the largest set of equations which cannot be decoupled has 2N complex unknowns and is a symmetric complex matrix.

-12 - 3.2 The Matrix Elements Z.. 1j The individual elements of the submatrices in (1.5) and (1.6) will be expressed in terms of a number of auxiliary functions of the given parameters a and 6. These functions are given in Tables II and III and formulas (2.1) - (2.23). In addition to the subscript index numbers previously introduced (i,j=l...N, m=l...M) there will also be used p,q=l,...N and ==1,...L where N, M and L are all integers whose values must be determined by numerical experimentation but are fixed in a given computer run. To distinguish integer values i and j from the symbol for /T1 we use j = /-1. The auxiliary functions are R(r,r') = R = (r2 + r12 - 2rr' cosu + 62)1/2 (2.1) R(a,a') = [2a2(1 - cosu) + 62]1/2 (2.2) D(r,r') = D = R (r2 + r2 - 2rr' cosu)1/2 (2.3) 6=0 j 0 if i / j i if i = j (The Kronecker 6) (2.4) oip = (a/N)[i - 1 + (p-l)/2] (2.5) B = (a/N)(i - 1 + p/2) (2.6) iP yip = (a/N)[i - (5/4) + (p/2)] = (1/2)(ai +. ip) Ap -lpi -'p (2.7)

-13 - = (T/)Q-(1/2)] (2. 8) F p;j,q) ~i~p I jq + ai~p ftjq) - F(f3,cp Itjq) - F (ca,fp jq) (2. 9) F(r,r')= r'.Zn(r-r' cosu + R) + r kn(r'-r cosu + R) ~ 2(inu/6 T 12tan 1[r 6sinu/[(r-r' cosu)R + D 2] + r2 tan- [r6, sinu/[ (r'-r cosu)R + D 2]]} + (D 2 cos- 2rr' sin 2u){TR - D - 26 k 6R D 1 -(sin 2u/6 2{r' 3Zn [(r-r' cosu + R) /(r-r' cosu + D)] + r3 kn[(r'-r cosu + R)/(r'-r cosu + D)] (2.10) G(jjq) =GU3.j ) - G(a. ) (2.11) G(r) = k~n[r-a cosu + R(a,r)] + (a/6)sinu tan- {a6 s~;inu/ [D (a, r) + (r-a cosu) R(a, r)] + kn{[c3 + R(a,r)]/D(a,r)}/6 (2.12) H(i,,P;j,rq) =HU3.,6. ) + H(cxa r t ip jq ip jq H H(~..cx..- ) - H (ao 4. (2.13) i~p jq jip jq H~~' r'Zn [(r-r' cosu + D) /(r-r' cosu + R)] *H(0,0) = - U 6sinu

-14 - + rkn [(r-r' cosu + D) /(r'-r cosu + R)] - (6/sinu)tan ' [(R6sinu)/(rr'siu + 6 2 cosu)]}/62 K(U)= (2/6).Zn[(R(a,a) + 6)/D(a,a)] 2 - (2/6 )[R(a.,a) - D(a,a)]-j S(i,p;j,q)= S(~.,p".J) + S(cy,.,ctjq S(ip ojq) - ip rjq) S(r,r')= rzn(r'-r cosu + D) + r'9kn(r-r' cosu + D) Table II (2.14) (2. 15) (2.-16) (2. 17) Mean Values of the Radial Functions

-15 - The following combinations of the functions in Table II occur frequently here and in Section 4. 2 F. (m) = F. = (C. + D. - m E. )(N/A) (2.18 iq iq iq iq iq ) Gq(m) ip Hcq (m) lp jq _ 2 =G = m AA. + B B. ip ip jq ip 3q = H = m(A. B. + B. A. ) lp ip jq lp ]q (2.19) (2.20) (2.21) Pi (m) = mAiq + Biq, Miq(m) = mAiq - Bi iq i-q i- q iq iq Table III Summation Limits p0 and q0 i P0' q0 < N - 3 N - 2 N - 1 8 6 4 N 2 Further auxiliary functions are Po ij = l p=l q -jD(Y.,Y fe ip J [(1 + jD(yi.'Yjq) ) q=l -p jq 2 2.q F(i,p;j,q) - j a /(2N) ][Gq cos mu cos u ip + Hjq sin mu sin u]} ip u=uZ and 0 q0 -jD(y.,pY ) ij = Z {e ip Jq [(1 + jD(Yip Yj )) ij p=l q=l - q - 2 2 jq F(i,p;j,q) - 3 a /(2N) ] [Gp sin mu sin u + cos mu cos u] + Hjq cos mu cos u]} ip u=uS (2.22) (2.23)

-16 - Note that..ij is obtained from.ij by interchanging Hjq and Gj. In terms of all the above auxiliary functions the impedance submatrices are TMO, TMO 2 L Z.. (TMO, TMO) = 7r6 /(2L) Y cos m u 13 Z& 9= 1 PO q0 p=l q=l -jD (y ip. q) e 'p F. F. [(1 + iq jq jD(yip,yjq)) (F(i,p;j,q) - 2H(i,p;j,q)) - j a / (2N) 2]} (2.24) TEO, TEO Z.. (TEO, TEO) = 76 /(2L) i] L {..i - e aa - "i cos mu 2 + + (1/2)/2) + m K(u) [(1/2)6 + 6 + (1/2)1 I[(1/2)6 + N N-2 N 6-1 N-1 + (1/2) 6 _2] u=u (2.25) TEO, TM1 Zj. (TEO, TM1) = [(1/2)6N + 6N1 -N N-1 Tr6 /(2L) L x { ij. + m cos mu Z=l 1] q 0 -jD + (1/2)6N-2] 2 {e q=l (a,.q) jq [ (1 + ]jD(a, yjq ))G(j,q) - j a/(2N)]} -mNejD(aa) cos mU-K(u)[(1/2)6N + 6N-1 (63 - 6-2) (6N N 2 U=UP TM1, TM1 *2 L Z.. (TM1, TM1) = 7r6 /(2L) 2 {&ij Z=l1 + (1/2)6N12] (2.26)

-17 - PO ^q -cos mu[ ~ C [e + jD(yiplYjq)) p=l q=l p 3 2 2 S (i,p;j,q) - a /(2N) ] FipFq PO -D(a,Yip) -N [ [e (1 + jD(a,yip))G(i,p) p=l -j a/(2N)] (6S - 6N-2)Fip q0 -jD(a,yjq) -N q [e 1 + jD(a,yjq ))G(j,q) q=l i i -j a/(2N)] (6N - N-2) Fjq N2e-j(a K(U) - ) - 6 2)] }u N N- 2 N N-2 U'uU (2.27) 3.3 The Matrix Elements W.. 13 The W.. will be expressed in terms of the integrals I )(a) = 1 xndx 1 x dx with explicit expressions as follows provided a7O. In this 0 (x+a) section a is not the radius of the disc which is denoted by a. m=l Il = nl+a 11 a 1 ln -n-l I(nIl) n=2,3,4,5 m=2 - 1 21 a(l+a) 1 l+a I = - ---- + n - 22 l+a a I 1 _ (n-l)a I n=3,4,5 2n (l+a) (n-2) n-2 I2(n-1)' n=3,4,5

-18 - m=3 1 1 1 I31 [ 2 - -2 ] 31 2 2 2 (l+a) a _ (4+3a)a 3 l +a 33 2 2 a 2 (l+a) 1 (n-l)a n-2 3n 2(n-3) n-3 3(n-1)' m=4 1 1 1 41 3~ 3 ~ 1 3 (l+a)3 a 2 3 3a+(a/2)a +(11/6)a 11 l+a _I_____ --- —- - - + 9n - 44 (+a) 3 6 a (l+a) 1 1 (n-l)a 4n (l+a)3 n- 4 n-4 I(nl) n=2,3,5 (3.4) These integrals will be used-to express a number of functions f1' *f4' gl9..g3' h1, h2 and W.i+ as shown in Table IV. Table IV Structure of W.. 13 i W. W W W Wi1i Wii+l i,i+2 ii+3 N f 0 0 0 N-l f +f N-1 2 1 0 0 N-2 f +f +f3 g +g2 h1 0 < N-3 f f +f+f g +g2g3 h+h2 W 1 2 3 4 g1+g2+g3 h1+h2,i+3......... I............. The f,,g. and h. functions are different for each submatrix. In their I 1 formulas we have not written the argument a explicitly in each I (a). Note that a takes on different values in different formulas. Two

-19 - special cases corresponding to a=O occur. These are for fl when i=l and are treated separately in footnotes. TMO - TMO 2 2 f = 76 (N/a)2 [(3/2)+a+I -m (I13+I24 13 1324 + (m /4) I35 a=i-l,iil 2 f2 = (N/a)2 [a-(l/2)+Il- 2I12+ 113 + m (Ill+ 212- I13- I22+ 3I23- 24 11 12 13 21 22 23 24 + (m /4)(I31+ 4I32+ 233-434+ I35)]a=i f3 = r6(N/a)2[(3/2)+ a+ 113 + m (2Il- I1+ 2122 24) + (m /4)(4I31- 4133+ I35)]a=i+l 2 f4 = r (N/a) [-(1/2)+ a+ Ill- 2I12+ I13 m(- + m2 ( -Ill+ 212 - I13+ I2 1 3 323I24) + (m /4)(I31- 4132+ 6133- 434+ I35)]a=i+2 gl = r6 (N/a) [-(1/2)- a+ Ill-I12 -(im2/2) (Il+ 211- 213+ I2+ 313- 21) /)( 11 2I12 2I13+ I22+ 3I23- 224) + (m /4) (I33+ 2134- I35) a=i g = 76 (N/a) [(1/2)+ a- I12+ I13 *When a=i-l,i=l so that a=O f = (7T6/2) (N/) 2[4-2m2+(m4/4)]

-20 - (n2 /21 - 21 - 21 + 31 +31 -21 + (m2/2) (3Il+ 2112- 2113- 2121 322 23 24) + (m4/4)(213+ 4I - 3. - 21 I 31 32 33 34 23) a=i g3 = 6(N/a) [-(1/2)- a+ I12- 13 + (m2/2)(-Il- 2I1+ 2I1+ 21- I2 3I+ 2124 11 12 13 21 22 23 24 + (m4/4)(2I31- 413+ I33+ 2134 I35)a=i+2 h1 = '6(N/a) [-(3/2)- a- 113 + im (-Il+ I23 I22+ I24) + (m /4)(2I33- 35)]a=i+ h = r6(N/a) [(1/2)- a- I+ 212- 13 2 11 12 13 + m2(-212+ I+ 212- 31 + 24) 12 13 22 23 24 + (m4/4)(I- 433+ 434 I35)a=i+2 31 33 34 35 a-i+2 W. = r6(N/a)2[-(1/2)+ a- I12+ 113 1,i+31 -(m/2) (I - 212+ 21+ I22 323+ 2124) 11 12 13 24 + (m4/4) (I33- 234+ I35)]a=i+2 TEO - TEO and TM1 - TM1 These two cases have identical W. f= T [(1/4)+(a/3)a/3)+(m /4)1 ] 15*When a=i-, i=l so that a=0 *When a=i-1, i=1 so that a=0 f = (/4) [1 + (m2) f = (Tr6/4)[1 + (i /4)] 1

-2 1 - f2 = Tr6 [(I/12) + (a/3) + (m2/4)(Iii+ 4112+ 2I13 RIi4+1I5 H~ f3 = 'rr6[(1/4) + (a/3) + (M2/4) (4I - 4113 115a i+ f4= 'ir6[(1/12) + (a/3) + Cm2/4( 41 + 61 - 41 + I H /4) 1i 12 13 14 15 ~a-i+2 = 7r6 [(1/12) + (a/6) + (M2/4)(113+ 2114 I15aH 92= 'ir6[-(1/12) - (a/6) + (M2/4) (2Iii+ 2I12 2I1- 2114+ 115 Ha+ 93= 'ir6[(1/12) + (a/6) + (M2/4) (21Ii- 4112+ 11 + 2I14 H1~-+ h = 1T6[-,(1/4)-(a/3) 1 + C /4(213~ 115 a=i+1 h2 = 'ir6[-(1/12) - (a/3) 22 + (m2/4)I ~11- 4I13-I1 ~+ Wi,i1+3 = rr6 [- (1/12) - (a/6) + (m2/4) (I13 21 + Is)~ 13 14 1 ~+ TE~ TM1 =1 7r6m/4 f 2= rr6m (3/4 )

-22 - f3 f4 4 = -7rmm(3/4) = -Trc6m/4 = T r6m/2 = 0 - -7rc6m/2 - Tr6m/4 = -7T6m/4 w.. = 0 1, 1+3 In this case (TEO-TM1) the results of applying Table IV may be summarized as in Table V. This can give an easy check on a bit of the programming. Table V W.. for TEO-TM1 1,j Case w i N N-l N-2 <N-3 Wi. 1/4 1 1/4 0 i,i Wi O 1/2 1/2 0 W +2 0 0 1/4 0 i,i+3 2 i,i+3 4. The Excitations Ik(0,10 ) The index k is determined by the scheme in Table I. It will facilitate writing the subsequent formulas to introduce the notation m = 6N(m-l) (4.1)

-23 - In the following J (x) indicates the Bessel function of the first kind of order m and argument x. We will use the following auxiliary functions IIk - IIIk to define the Ik(0,0). For excitations we use 6=e0, q=q0 but the same Ik functions, evaluated for other 0, values, are needed in Section 5. Po IIk(0) = F(m) Jm(ykp sin (4.2) p=l kp P0 (II )kp n(e) - Mk Pm)JJ +(y sin))} IIIk) = p- P ( m)-(Ykp sine) + Mk(m) Jm+l(kP sin)} p=l p-1 kp m k mil k Here Fk (m), Pkp(m) and M p(m) are given by (2.18) and (2.21); the constant n is defined as n = 6(-)m-1 a/(2N) (4.3) Each Ik is decomposed into two parts with factors p and s, 0 <Ip[< 1, 0 <j|s| 1, |PI2 + Is12 = 1, and p and s either real or imaginary. Then Ik(,) = pIP(0,) + sI(0,() (4.4) (The symbols p and s stand for the German words parallel and senkrecht meaning parallel and perpendicular.) The IP and Is are specified as k k follows: 1 + m < k < N + m (TMO ODD) IP(0,) = 2j sinO sinm4 rIk(0) k I1k Ik (e,) = 0 k

-24 - N + 1 + m < k < 2N + m (TEO ODD) 1k(e, ) = Ik (',() = 2N + 1 + m < k If(e,~) = I (e = 3N + 1 + m < k cosO sinme IIIk(e) cosmk IIIk,(e) < 3N + m (TM1 ODD) cosO sinm) IIII(6) k' = k-(n+m) k' = k-(2N + m) cosmp IIIk(e) < 4N + m (TMO EVEN) --- I (,, f) = 2j1 sine cosmP IIkz(6) k' = k-(3N + m) s (e, ) k4 + 4N + 1 + m < = 0 (4.5) k S 5N + m (TEO EVEN) - i (e, ) 5N + 1 + m < = cose cosm I IIk,(8) = -sinmp IIIk,(9) k' = k-(4N + m) k < 6N + m (TM1 EVEN)... i (e, ) Ik (e,~) 1k9<> = cose cosme I IIk,(e) k' = k-(5N + m) = -sinme IiIk,(e) 5. Cross-Sections Once the ak are computed for a given set of data (80, ~0, p, s; a, 6, n), one or more of the quantities given below will be desired. The formulas use notation |AI I (A) and A* for absolute value, imaginary part and complex conjugate of complex quantity A.

-25 - Bistatic or differential cross-section M 3N+m B (6 ) = (4rk)-2 {I I [ a kIP ( 0L (_ ' k s s m=l k=l+m 6N+m + akIk (9s', ] I (5.1) k=3N+l+m M 3N+m + I I [ I a I ^ akIk (,rs ) m=l k=l+mk s 6N+m + k kI (s' s ) 1 } k=3N+l+m Eq. (5.1) could be expressed more compactly by not grouping the summations over k into subsums of 3N terms, but this grouping is necessary in (5.2) and in fact corresponds to the data storage procedure used in our computing code. To determine the integrated or total scattering cross-section the integration over cs can be performed analytically leading to the result M 3N+m gB(es) = 7 (4k0)-2 [I [ akk (e,T/2m) 2 m=l k=l+m kk s 3N+m 6N+m + I kIk ( 0) +! I (e 0) k=l+m k=3N+l+m k k s 6N+m * + N+ kIk (e,T/2m) 2]. (5.2) k= 3N+l+m The integration over 8s must be done numerically leading to the following results: Total Scattering Cross-Section iT &TS = aB(Os)sines d (5. 3) 0 B 5

-26 - Absorption Cross-Section 2,2 -22K K a = -I (n )In -_I ko il j- c. W.i - A m 0'_Ij1131 Total or Extinction Cross-Section UT = TS + A (5.6) Assymmetry Factor IF UASY = B(S) sines cos(Os-e0)dO (5.7) The forward scatter theorem* provides a relationship between aT and the component of the far-zone scattered electric field which is parallel to the incident electric vector (assuming plane incident polarization) and which is propagating directly forward; that is, in the same direction as the incident wave so that O =O90, %=QO' This relationship provides a check on aT as determined by (5.6) although one would expect its accuracy to be less for approximate numerical solutions such as ours. This is true partly because the theorem hangs the entire determination of aT on the value of the field in one direction; there is no smoothing of errors by integration as there is in (5.3) which provides one term in (5.6). A second reason is that the theorem is based on the interactions between excitations of the various elements of the scatterer. It will not apply for an approximate theory which does not take these into account sufficiently accurately. For example it does not hold for lowest order Rayleigh theory (Van de Hulst, 1957) and would not hold for our analysis if we had not adequately taken account of interelement interactions in the impedance evaluations. In any case we have programmed and evaluated the following *Formula (110) Section 13.5 of Born and Wolf (1959 or 1975). We will apply the theorem for elliptical incident polarization and so use a slight generalization of the derivation and formula in Born and Wolf since they considered only plane incident polarization.

-27 - formula which is derivable directly from the forward scatter theorem. We use the notation of Section 4. Extinction Cross-Section 2 3NM 3NM (T = k0 I{P*( + q* + ~kk ~ 0 ) ( 5.8 ) 4. Computational Code and Results A skeletal description of the FORTRAN programs for numerically evaluating the formulas detailed in Section 3 follows. All equation numbers refer to the numbering in that section. Programs ZMAIN and WMAIN govern computation and storage of the symmetric complex matrices Z and W respectively. In each case the matrices are packed for storage into vectors using column major order using separate vectors corresponding to each of the 3N x 3N submatrices for each m value. The matrices were stored since the slowest and most expensive part of the program lies in computing the Z elements, yet they are needed for each computation involving the same disc dimensions relative to wavelength regardless of different incident wave amplitudes, directions and polarizations or different disc materials. The W matrix elements, although inexpensive to compute are similarly reusable. Each of the submatrices is calculated by separate subroutines: ZTMO, ZTM1, ZEOM1, ZTM1 for Eqs. (2.24) - (2.27) corresponding to the different mode combinations for Z elements and WTMO, WTEO, WEOM1, WTM1 for Eqs. (3.6) - (3.9) for the W elements. The excitations, Eqs. (4.5) are computed as needed for specified incident waves by program IMAIN using function subprograms for the sums

-28 - of Eq. (4.2). The excitation values are not stored since they are very inexpensive to compute. Program AMAIN uses the stored Z and W elements plus the excita-.m tions computed by IMAIN to solve Eq. (1.1) for the coefficients ak with the help of subfunctions CSPCO and CSPCL available on the Michigan Terminal System (MTS). CSPCO computes the condition number and a 2 decomposition of matrix z = Z - W / (n - 1) into the form UDU* where U is a product of elementary unit lower triangular and permutation matrices and D is a block diagonal matrix with blocks of order 1 or 2. CSPSL then solves for the ak by back substitution. Finally, program XSECTS governs computation and printing of the various cross-sections, Eqs. (5.1) -(5.8). The cross sections, and all auxiliary functions are computed in function subroutines. The integrations involved are done by calling subroutine QCRP, an adaptive quadrature scheme available in the MTS system. Apart from debugging the program we have had time and money on the grant to obtain only a few illustrative numerical results. Tables l(a)-(c) list normalized*cross-sections for a disc of normalized dimensions a = 2.95, 6 = 0.1 irradiated by plane polarized radiation of wavelength X = 3.0 pm and composed of ice. This wavelength corresponds to the center of an infrared absorption band for ice and the refractive index is n = 1.130 - iO.2273. To exhibit how the scattering was influenced by the absorptivity we also list aB (Os's) results for n = 1.130, I n = 0. m Results are presented for three incidence angles, 0 =0~ (broadside), 9 =45~ and 0 0 G0 = 90~ (edge-on), each incident wave propagating in the x-z (q=0~) plane with electric vector parallel to the y axis. For each incidence *Length units are kO

-29 - angle three scattering patterns are presented; i.e., a (0,8 ) as a BD s s function of 0 for fixed values of 0; s = 0~, 90~ and 180~. When s s s = 90~ the forward scatter directions correspond to p = 180~ and o s the backscatter directions to =0~. We also list the values of 5 aTS and aT; the latter computed both by Eq. (5.6) and Eq. (5.8), and the assymetry factors. We note the agree.aent. in the a values obtained by integration and T by the forward scatter theorem. In light of the discussion just prior to (5.8) it may therefore be inferred that the computed induced current is quite accurate. The listed data clearly show various flat disc scattering characteristics which are radically different from those of particle scattering. Thus for broadside incidence the scattering is symmetric in the forward and backward hemispheres corresponding to the 0~ < 0 < 90~ and 90~ < 0 < 180~ results while in the plane of the s s disc in the direction of E the scattering drops almost to zero. The resemblance to Rayleigh scattering from a sphere, simple dipole radiation is strong. On the other hand there is an appreciable forward to backward assymmetry for edge-on incidence; a forward bulge and a depression in the direct backscatter direction. There is also only a mild dip in the scattering pattern in the plane of the disc in the direction of the incident E vector; i.e., for E = s = 90~. Thus 0 0 there is no vestige of the independence to particle orientation which is inherent for spherical scatterers. The pattern for broadside incidence can be thought of as due to a superposition of parallel thin dipole currents in the y direction induced by the incoming electric field and all excited and radiating almost in phase because 6 << 1. The effects of the loop currents generated by the incident magnetic field is negligible since oppositely

-30 - directed current elements on the top and bottom surfaces can be thought of as giving mutually cancelling radiations because they too are radiating in phase. For edge-on incidence the pertinent phase differences for excitation and radiation are not small since the disc diameter is about a wavelength. Hence the y directed dipoles due to incident E are not excited in phase so that forward to backward symmetry should not be expected. Similarly the oppositely directed current elements in the current loops induced by the incident magnetic field need not produce mutually cancelling fields. In this way one can explain the "fill-in" of the pattern in the plane of the disc (relative to the almost zero result purely y directed current elements would yield) as being due to the circumferential current components.

Table 1. Normalized Cross Sections for disc, a=2.95, 6=0.1 with incident electric polarization parallel to flat disc surface (a) Broadside direction of incidence, 0=0~ BISTATIC CROSS SECTION aB ( 8ss ) B> S n = 1.130 - iO.2273 s n = 1.130 0 degs. degs. s 0 90O 180~ 00 90~ -i I 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 0.169E-01 0.159E-01 0.133E-01 0.998E-02 0.689E-02 0.451E-02 0.293E-02 0.201E-02 0.154E-02 0.140E-02 0.154E-02 0.201E-02 0.293E-02 0.451E-02 0.689E-02 0.998E-02 0.133E-01 0.159E-01 0.169E-01 0. 169E-01 0.154E-01 0.115E-01 0.707E-02 0.363E-02 0.156E-02 0.561E-03 0.163E-03 0.295E-04 0.362E-14 0.295E-04 0.163E-03 0.561E-03 0.156E-02 0.363E-02 0.707E-02 0.115E-01 0.154E-01 0.169E-01 0.169E-01 0.159E-01 0. 133E-01 0.998E-02 0.689E-02 0.451E-02 0.293E-02 0.201E-02 0.154E-02 0.140E-02 0.154E-02 0.201E-02 0.293E-03 0.451E-02 0.689E-02 0.998E-02 0.133E-01 0.159E-01 0.169E-01 0.453E-02 0.426E-02 0.357E-02 0.269E-02 0.187E-02 0.124E-02 0.816E-03 0.564E-03 0.437E-03 0.398E-03 0.437E-03 0.564E-03 0. 816E-03 0.124E-02 0.187E-02 0.269E-02 0.357E-02 0.426E-02 0.453E-02 0.453E-02 0.411E-02 0.307E-02 0.189E-02 0.978E-03 0.423E-03 0.153E-03 0.477E-04 0.818E-05 0.102E-14 0.818E-05 0.447E-04 0.153E-03 0.423E-03 0.978E-03 0.190E-02 0.307E-02 0.411E-02 0.452E-02 ABSORPTION CROSS SECTION 0.145E+01 1.61 DB 0.0 TOTAL SCATTERING CROSS SECTION 0.457E-01 -13.4 DB TOTAL CROSS SECTION EQ. (5.6) 0.149E+01 1.64 DB TOTAL CROSS SECTION EQ. (5.8) 0.151E+01 1.80 DB 0.124E-01 0.124E-01 0.183E-01 0.768E-09 -19.1 DB -19.1 DB -17.4 DB ASSYMMETRY FACTOR 0.743E-09

Table l(b). 45~ direction of incidence, e = 45~ O BISTATIC CROSS SECTION a(O,4s ) B S S n = 1.130 - i0.2273 s 0 s n = 1.130 s degs. 0 90O 180~ 00 90~ 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 0.560E-02 0.745E-02 0.901E-02 0.996E-02 0.102E-01 0.991E-02 0.932E-02 0.872E-02 0.829E-02 0.814E-02 0.829E-02 0.872E-02 0.932E-02 0.991E-02 0.102E-01 0.996E-02 0.901E-02 0.745E-02 0.560E-02 0.560E-02 0.516E-02 0.411E-02 0.299E-02 0. 216E-02 0.168E-02 0.143E-02 0.130E-02 0.123E-02 0.121E-02 0.123E-02 0.130E-02 0.143E-02 0.168E-02 0.216E-02 0.299E-02 0.411E-02 0.516E-02 0.560E-02 0.560E-02 0.384E-02 0.244E-02 0.148E-02 0.884E-03 0.546E-03 0.362E-03 0.266E-03 0.220E-03 0.206E-03 0.220E-03 0.266E-03 0.362E-03 0.546E-03 0.884E-03 0.148E-02 0.244E-02 0.384E-02 0.560E-02 0.153E-02 0.201E-02 0.241E-02 0.265E-02 0.272E-02 0.264E-02 0.249E-02 0.233E-02 0.222E-02 0.218E-02 0.222E-02 0.233E-02 0.249E-02 0.264E-02 0.272E-02 0.265E-02 0.241E-02 0.201E-02 0.153E-02 0.153E-02 0.140E-02 0.111E-02 0.785E-03 0.550E-03 0.417E-03 0.351E-03 0.317E-03 0.299E-03 0.293E-03 0.299E-03 0.317E-03 0.351E-03 0.417E-03 0.550E-03 0.785E-03 0.111E-02 0.140E-02 0.152E-02 ABSORPTION CROSS SECTION 0.114E+01 0.55 DB 0.0 TOTAL SCATTERING CROSS SECTION 0.330E-01 -14.8 DB TOTAL CROSS SECTION EQ. (5.6) 0.117E+01 0.67 DB TOTAL CROSS SECTION EQ. (5.8) 0.117E+01 0.67 DB 0. 878E-02 0. 878E-02 0. 723E-02 0.428E-02 -20.6 DB -20.6 DB -21.4 DB AS SYMMETRY FACTOR 0.162E-01

Table 1(c). Edge-on direction of incidence, 0 = 90~ 0 BISTATIC CROSS SECTION aB(9,6 ) n = 1.130 - iO.2273 0 s degs. pS n = 1.130 _s 0 90~ 180 0~ 90~ 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 0. 140E-02 0. 252E-02 0.407E-02 0. 583E-02 0. 746E-02 0. 870E-02 0. 946E-02 0. 983E-02 0. 996E-02 0. 999E-02 0. 996E-02 0. 982E-02 0. 946E-02 0. 870E-02 0. 746E-02 0. 583E-02 0. 407E-02 0. 252E-02 0. 140E-02 0. 140E-02 0. 135E-02 0. 128E-02 0. 129E-02 0. 126E-02 0.129E-02 0. 129E-02 0. 126E-02 0.122E-02 0.121E-02 0. 122E-02 0. 126E-02 0.129E-02 0.129E-02 0. 126E-02 0. 124E-02 0. 128E-02 0.135E-02 0. 140E-02 0.140E-02 0. 729E-03 0. 398E-03 0. 258E-02 0.212E-03 0.207E-03 0. 218E-03 0. 233E-03 0. 242E-03 0.246E-03 0.242E-03 0. 232E-02 0. 218E-03 0. 207E-03 0.211E-03 0.258E-03 0. 398E-03 0. 729E-03 0. 140E-02 0. 398E-03 0. 694E-03 0. 110E-02 0. 157E-02 0. 200E-02 0. 233E-02 0.253E002 0.263E-02 0.267E-02 0.268E-02 0.267E-02 0.263E-02 0.252E-02 0.233E-02 0.200E-02 0.157E-02 0.110E-02 0.694E-03 0.398E-03 0.398E-03 0. 380E-03 0. 343E-03 0. 316E-03 0. 310E-03 0. 313E-03 0. 312E-03 0. 303E-03 0.293E-03 0. 289E-03 0.293E-03 0.303E-03 0.312E-03 0.313E-03 0. 310E-03 0. 316E-03 0.343E-03 0.380E-03 0. 398E-03 ABSORPTION CROSS SECTION 0.114E+01 0.55 DB 0.0 TOTAL SCATTERING CROSS SECTION 0.231E-01 -16.4 DB TOTAL CROSS SECTION EQ. (5.6) 0.116E+01 0.65 DB 0. 612E-02 0. 612E-02 0.516E-02 0.4792-02 -22.1 DB -22.1 DB -22.9 DB TOTAL CROSS SECTION 0.116E+01 ASSYMMETRY FACTOR 0. 182E-01 0.65 DB

-34 - References 1. Auer, August H. Jr., and Donald L. Veal, (1970), "The Dimensions of Ice Crystals in Natural Clouds", J. Atmospheric Sciences, 27, pp. 919-926. 2. Born, Max and Emil Wolf, (1959), (1975), "Principles of Optics", Pergamon Press. 3. Chu, C-M., and Herschel Weil, (1976), "Integral Equation Method for Scattering and Absorption of Electromagnetic Radiation by Thin Lossy Dielectric Discs", J. Computational Physics, 22, pp. 111-124. 4. Van de Hulst, (1957), "Light Scattering by Small Particles", Wiley, p. 66. 5. Weil, Herschel and C-M. Chu, (1976), "Scattering and absorption of electromagnetic radiation by thin dielectric discs", Applied Optics, 15, pp. 1832-1836. 6. Zikmunda J. and Gabor Vali, (1972), "Fall Patterns and Fall Velocities of Rimed Ice Crystals", J. Atmospheric Sciences, 29, pp. 1334-1347.