AN ASYMPTOTIC SOLUTION FOR THE TWO-FREQUENCY MUTUAL COHERENCE FUNCTIONS OF A RANDOM SLAB by Hsiao-fei Maa A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1985 Doctoral Committee: Professor C. M. Chu, Co-Chairman Professor V. V. Liepa, Co-Chairman Professor E. N. Leith Professor V. C. Liu Professor C. T. Tai RL-795 = RL-795

RL 795 ABSTRACT AN ASYMPTOTIC SOLUTION FOR THE TWO-FREQUENCY MUTUAL COHERENCE FUNCTIONS OF A RANDOM SLAB by Hsiao-fei Maa Co-Chairmen: C. M. Chu, V. V. Liepa The parabolic equation method has been mainly used in dealing vwith scattering problems associated with communications. The backscattering effects are often ignored in these applications. However, in remote sensing problems the calculation of the backscattered field intensity is essential, and neglecting the effect posed restrictions on the application of that method. In this thesis a modified parabolic equation method is employed to formulate the two frequency mutual coherence functions of a statistically homogeneous body of a continuous random medium. The backscattering effect is included in the formulation. An asymptotic solution is found and the coherence functions of a flat slab for both the transmitted and the backscattered fields are solved. The scattering effects of the slab upon a rectangular pulse are calculated and plotted. The results are found to be physically meaningful, the energy conservation law is observed by the solutions, and they represent a realistic picture of the scattering mechanism; moreover, this formulation also provides a means to evaluate the backscattering of a random medium, thus is suitable for both communication and remote sensing applications.

To My Grandmother ii

ACKNOWLEDGMENTS I would like to thank the two co-chairmen of my doctoral committee, Professor Chu and Professor Liepa. for their help and guidance. I am thankful to Professor C. T. Tai, Professor E. N. Leith and Professor V. C. Liu for their useful suggestions. as well as the time they devoted serving on my dissertation committee. A special 'Thank You' goes to Dr. Leo L. Beranek of Boston, Massachusetts, for his invaluable help and encouragement. iii

TABLE OF CONTENTS DEDICATION...................................................... ii ACKNOWLEDGMENTS............................................ iii LIST OF TABLES........................................... vi LIST OF FIGURES........................................... vii LIST OF SYMBOLS............................................. viii CHAPTER I. INTRODUCTION............................................. 1 1.1 Random Media and Associated Problems 1.2 The Coherence Functions 1.3 Continuous Random Media Problems 1.4 Parabolic Equation Method II. FORMULATION........................................ 6 2.1 Equations Governing the Electric Field 2.2 Fluctuating Quantities 2.3 The Closure Relations 2.4 The Coupled Equations 2.5 Boundary Conditions 2.6 Summary of Chapter II III. THE ASYMPTOTIC SOLUTIONS.............................. 19 3.1 Basic Approach 3.2 The Green's Function 3.3 Green's Function Properties 3.4 Integration of the Equations 3.5 Stationary Phase Method IV. SOLUTIONS FOR THE SLAB.............................. 3 4 4.1 Laplace Transform Method 4.2 The Matrix (sT - B) 1 4.3 The General Solution 4.4 Solution for a Single Layer 4.5 Conventional PEM Result 4.6 Summary of Chapter IV V. A CALCULATED EXAMPLE.................................. 56 5.1 Basic Formulae 5.2 Fourier Series of the Pulses 5.3 The Scattered Pulses 5.4 Discussion and Summary iv

SUM M ARY.................................................. 74 APPENDIX.................................................... 75 REFERENCES................................................... 85 v

LIST OF TABLES 5-1. Parameters Used in the Example.................................59 5-2. Energy Relations of the Pulses.................................. 62 5-3. Accuracy of Eqs. (5-5) as Energy Level Estimators....................63 vi

LIST OF FIGURES 2-1. Structure of the Slab........................................ 18 4-1. Physical Layout of the Problem.......................................... 55 5-1. Scattering of a Rectangular Pulse (Data set No. 1)........................... 65 5-2. Scattering of a Rectangular Pulse (Data set No.2)................................ 66 5-3. Scattering of a Rectangular Pulse (Data set No.3)........................... 67 5-4. Scattering of a Rectangular Pulse (Data set No.4)........................... 68 5-5. Scattering of a Rectangular Pulse (Data set No. 5)............................... 69 5-6. Scattering of a Rectangular Pulse (Data set No.6)........................... 70 5-7. Transmitted Pulse: comparison (Data set No. 1).................... 71 5-8. Transmitted Pulse: comparison (Data set No.2)........................... 72 5-9. Transmitted Pulse: comparison (Data set No.4).......................... 73 vii

LIST OF SYMBOLS r (u,v) coherence function of the slab for forward scattered field r (u,v) coherence function of the slab for backward scattered field p horizontal position vector, p = xx + yy 4 (kl,k2) coherence function for the complex envelope of the field A(p) auto-correlation of the permittivity fluctuations in the horizontal plane A,A(0) A(p) at p = 0 E(r) total scalar electric field, E= E +EE (z, p) forward propagating electric field, 0 z L E (0,p) incident field at z= 0 E (z, p) backward propagating electric field, 0< z < L E (L, ) incident field at z = L G Green's function for V+ G Green's function for V H Green's function for W H Green's function for W k free space wave number, k=L thickness of the random slab r position vector, r = xx + yy + zz T (L) transmitting matrix of the random slab viii

t.. element of matrix T 13 u dimensionless variable, u=(k + k2)A(0) v dimensionless variable, v = (k -k2)A(O) V correlation function of E + V <E E > V correlation function of E, V = <E E > ++ + - - W cross-correlation, W = <E E > W cross-correlation, W - <E E > ix

CHAPTER I. INTRODUCTION 1.1 Random Media and Associated Problems By) definition, a random medium is a medium whose characteristics vary in a manner that they can not be decided exactly for a given position at a given time. The medium of interest could be either discrete or continuous. An example for the discrete case is a medium containing randomly distributed moving particles such as rain drops or hail, whereas one for a continuous model would be a body of air whose density due to temperature and wind velocity fluctuations., is a continuous random function of both space and time. When electromagnetic waves propagate in such media, fluctuations will be induced. Due to the uncertainty of the fluctuations, the only practical way to deal with the quantities of interest is to find their average values, and hence probabilistic methods are appropriate. Of all the quantities of interest, the average field <E> and the average power density represented by the quantity <EE > are undoubtedly the most important. A single particle or irregularity in space will cause scattering in all directions. When a cluster of particles or irregularities exist, multiple scattering generally occurs to make the problem more complex. Various approximate models have been developed to include multiple scattering effects to some degree. When the typical size of the scatterers is very large when compared with the wavelength of the field, the scattered energy will concentrate within a small angle relative to the propagating direction of the wave. In such situations, the so-called 'parabolic equation method' is widely used to calculate the forward propagating field quantities. This name stems from that under the 'small angle' assumption, the differential equations for the field quantities are of parabolic type. In the conventional parabolic equation method primarily used in communication 1

2 problems, the backscattering is often assumed to be negligible to simplify the formulation. Jedrzejewski proposed a formulation to include backscattering for the coherence functions for a single-frequency plane wave case[l]. In the present work, the two-frequency mutual coherence functions of the forward- and back-scattered fields is formulated and the scattering of a rectangular pulse by a flat slab of a continuous random medium calculated. 1.2 The Coherence Functions An electric field propagating in the z direction can be written as ikz E(k,r) = U(k,r) e where U(k.r) is the complex envelope, k is the fiee space wave number of the field, and the time dependence e is assumed. The two frequency mutual coherence function of the complex envelope is then defined as (kk2,rl,r2) = <U(kl,rl)U (k2r2)>, and for a statistically homogeneous medium 4(kl,k2,rl,r2) = 4(klk2, Irl-r21)At k =k2 k and r _ =_, the real part of F represents the average power spectrum of the field. If we consider a body of random medium as a system, the system response due to a delta input

3 (kl k2) = 4(k, k2 output plane) k input plane u if 2~1'(k k2 'iinput plane>=1 medium. The correlation function B(t,t2) is defined by ([2]) B(tlt2) = < Uout(t ) Uout(t2 > co co iltl+iw2t2, — 00 —00 = f fd 1 d2< in ( )Uin (W2)>FI u(l( 2)e1 1 2 2 where u.n(() = dtU. (t)e is the complex spectrum of the input field. in 2 dt inn From its definition, B(t,t2) is proportional to the power density of the output signal Pout(t) when tl=t2t. 1.3 Continuous Random Media Problems We start with the scalar wave equation 2 2 (-) [V2 + k e(r)]E(r) = 0, where e(r) is the relative permittivity of the medium. The medium is continuous and is varying slowly compared with the field. Depolarization effects of scattering are neglected. The medium permittivity can be written as E(r) = e + e(r) - 1 + e(r) where (r)< 1 is a random variable with a zero mean. It represents the fluctuating part of the permittivity. There are two general approaches to solve the problem. The first one is to solve Eq.(1-1) for E and then obtain <E> and <EE > from the expression for E. The second approach is to derive equations for <E>, <EE > or whatever quantities of interest from

4 (1-1) and then try to solve those equations. Because of the complexity of multiple scattering mechanism, assumptions and approximations have to be made for both approaches in order to arrive at meaningful solutions for specific situations. (Such as the magnitude of fluctuation, the size of scatterer, etc.) The first approach results in an expression of a complicated series. The common practice is to truncate the series after one or two terms. Since the convergence is very slow with even moderate fluctuations, this approach is suitable for very weak fluctuations and short propagation distances only. The difficulties encountered in the second approach stem from the fact that c(r) and Er) are mutually correlated, and therefore new variables like <e(r)E(r)> are created when equations such as (1-1) are averaged. Assumptions and approximations are usually made to circumvent this so-called 'closure problem', and despite those approximations, higher order scattering terms can be included in the formulation to make this approach suitable for larger permittivity fluctuations and longer propagation distances. 1.4 Parabolic Equation Method For a partial differential equation N2u N au Lu = Z Z a.ij(x,t)a + b. i(xt)- + c(x,t)u - = 0.lljJ 1 i=3 i. if matrix (a..) is positive definite, i.e., if for any real vector Q = (... ) ~ a..itit.>O, the operator L is of parabolic type. In the process of deducing equations for <E> or <EE >, this type of equations will be encountered if the scattering is assumed to be in small angles to the incident wave. According to scattering theory, the 'small angle' assumption is always true if the size of the permittivity irregularities is much larger than a wavelength. In this work a high frequency incident wave is assumed to justify the small angle approximation. The following is a simple example of the deduction of a typical parabolic equation

5 from Eq.(1-1). Let E(r) = U(r)e. Substituting it into (1-1) gives 82 8a22 a2 2 + i2ka)U(r) + a2 + ay + k 2(r)]U(r) = 0. (1-2) The 'small angle' approximation implies a-2U(r) << i2k zU(r) (1-3) and when it is combined with Eq.(1-2), a parabolic equation + + i2k- 2+ a2+ a + i + kE(r)]U(r) = 0 results. In the past. the parabolic equation method was mainly used to investigate the forward scattering problems associated with communication systems. The majority of those works either ignored the multiple scattering effects or assumed that the scattering occurs only in the forward direction. But models that include backscattering are needed for remote sensing problems. Jedrzejewski introduced a backscattering model and gave some numerical results of the mutual coherence functions for a single frequency case in his Ph.D thesis[l]. In the present work the two-frequency mutual coherence functions of the forward and backward scattered waves are formulated and analytical asymptotic solutions are presented. As a result, it is possible to calculate both the transmission and the backscattering of a random medium, as shown by an example in Chapter V. The physical structure for which the solutions are given is a flat slab of a continuous random medium with normal incident waves on both sides, the slab is located in free space, with thickness L in z-direction and is infinite in x and y directions. In the formulation, the forward and the backward scattered fields are mutually coupled, as a result, multiple scattering effects have been included.

CHAPTER II. FORMULATION 2.1 Equations Governing the Electric Field To account for the effects of backscattering caused by medium permittivity fluctuation, it is convenient to explicitly introduce a backscattered field E, such as that; in Jedrzejewski's work ([1]): (2-1) E(r,k) = E (r,k) + E (r,k). A second relation of E and E is chosen as a-E(r,k) = ik[E (r,k) - E (r,k)], (2-2) where k is the wave number in free space. The typical size of the disturbance is assumed to be much larger than a wave length so that the scattering is primarily in the forward and backward direction. Therefore E+ and E can be writter as + ikz (2-3a) E (r,k) = U(r,k)e and -ikz (2-3b) -ikz E (r,k) = U (r,k)e The field E(r,k) inside the random medium satisfies equation a2 2 [V2 + + k (l+e)]E(r,k) = 0, (2-4) where V 2 = + a2 a= 2 ay2 r r 6

7 e = <E > 1 r r e << 1 is a real quantity. The equations for E+ and E- are readily derived from (2-1), (2-2) and (2-4): ikLa 2 2 22V - [i2k +V2+k(2+)]E(r,k) = -(V +k E)E(r,k) (2-5a) a L 2 -2 2-5 +) [-i2k^a +V2+k2(2+E)]E (r,k) = -(V2+k2e)E (rk), (2-5b) whereas the complex conjugates E and E obey a 2 2 _( 2 2" -* [-i2kz +V2+k (2+)]E+*(r,k) = -(V +k )E (r,k), (2-5c) and [i2ka +V2 +2 E*2 2 [i2k - +V +k2(2+k)]E (r,k) = -(V2+k2 )E+(r,k). (2-5d) z + (2-5d 2.2 Fluctuating Quantities As an intermediate step essential to the task of this chapter, equations for the quantities E (r kl)E (rbk2), E rkl)E (rbk2) -a - -(r kl)E+ * - - * E (rak)E (rb'k2) and E (rakl)E (rbk2) are needed. By averaging those equations, the equations for the coherence functions can be derived. Take two points in a horizontal plane r = Pa + z, -a and Lb = Pb + zz. Multiply Eq.(2-5a) for point r and k by E +,k2), subtract the product of 2-a 1note E (r,k ) and Eq.(2-5c) for r and k, and notice that — a 1 b 2

8 a 2Ebr) = V2bE(ra) = 0 The resulted equation is. 2 v2 aL - i a b + + i +k+*+ +k E) = aZ 2.k a2b 6 a2 2bla2. 2 v2 i a - +* - b + -* j* - +..* (-E E -E E )+ -(k E E E -k E E) 2kla 2b k la 2b 2 2. a la 2b 2 b laE2br (2-6a) where ea = 6 (z,Q a) Eb = E(zPb ) El E(k,~ la 1. a Eb = E (k2,,b p = xx + yy. SimilarlyN v2 2 [a +ia - b - E-E* 2 2 V + *-V +* i + - (-~E Eb j-E12) 1(keE1 1 2 a2 a2b 2 1aL. 2 2 [a - a + b-~ k)E az k 2 laE2b2. 2 ia- - b + +* j (- (TE Eb + ~El E b+ 1(k eE E 2k1 l b kl2 2 la la,f(kl 6ElaE2b-k 26bElaE2b) = aE2bk2Eb laE2b) (2-6b) f(k9 6ElE2b.k 26bElE2b) = E~b+k~+ +*~b (2-6c) and

9 2 2 V V - i a bl + i * i- + + *)+i( [+ 2(k1 2 1k2 1k laE2b+ E2 la2b +k 2blaE2b V2 V2 i a + +* -* + +* - 2 k 1a2b k 2la 2b 2 (k a la2b 2 b la 2b Equations (2-6) are a set of mutually coupled partial differential equations with +- +; - -+-_' + -* - + random coefficients a and eb for quantities E E E E E E and E E The coherence functions are defined as the time-average of those quantities: V+ +(kz )> V (k,k,z,Z,Pb)= <E (k,Z,a)E (k2,z,pb)> V (kl'k2'za,pb) = <E (klZpa)E (k2,z,p)> - + * W (klkZ,,a ) = <E (kl,Z,a)E (kz,)> (k1,k2,ZpaPb) = <E (klZa)E (k2,z,pb)>. Of the above four quantities, only V and V (the auto-correlation functions of the forward and backward scattered fields) are of interest. The two cross-correlation functions, W and W-, will be treated as intermediate variables and be eliminated eventually. 2.3 The Closure Relations In the averaging process of Eqs.(2-6), variables other than the four correlation functions defined in the previous section, such as <e E + E >, are produced. It is a la 2b essential to express these quantities in term of the correlation functions so that a set of partial differential equations for the correlation functions can be derived. In an isotropic, homogeneous medium the permittivity correlation is sharply peaked in the direction of propagation and can be approximated by the expression ([3])

10 <e(r'')E(r")> = 6(z'-z")A(P'-~"), where A('-p") is the transverse correlation of the medium permittivity. Furthermore, since the medium is statistically homogenous, A(p'-p") = A("-') = A(p) P = P -"!I A(p-+o) = 0. In this work e is assumed to have a zero mean Gaussian distribution. Under these assumptions, the Furutsu-Novidov relationship can be applied to give ([5]) - 6E2 L 2y <fa E- > = Jdz' dp' { 6 (z-z')A( -p' ) [< E; >+<E - >] } a - a 6(r') 2- 1%(r') (2-7a) and E +* +;. L 3 5* + E2 <E aE E> = fdz' f d'{6(z-z')A(p - <') [< E2 >+<E, ->]}, 0 - a 2ry 1%E(r') (2-7b) where = e (z,'a) E3 = E(klZ,,o) E2 =E (k2,z,p ) Pa takes the value of either p~ or p, and 6E(k,z,^) is a functional derivative. 6E(z', P ')

11 In Eqs.(2-7) either the upper or the lower row of signs. but never a mixed combination, should be used for the entire equation. A discussion about the concept of functional derivatives and the derivation of that of a linear functional can be found in [4]. In our particular case, the expressions are 6E+(k,z,p) = 0 (z'>z) - = 0 (z'<z) 6E(z' r') and 6E (k,z,~) ie(Z, -) z+Z 6E (k,z,p) - \ z94z+ 6E(Z',P') ik = -[E (k,z,E')+E (k,z,p')]~i(p-p'). (2-8) 4 By substituting Eq.(2-8) into (2-7), the closure relations are derived. These are _+ +* i r A+ + <eE E2> = { [klA(Pa)-kA(p) ]- -kv a)W+kA(p (2-9a) a 12432 4L 1 a ay 2 and + + + < E- E2> = 4[klA(p )-k2A(pa )]W -k2A(Pa)V-+kA(pa)Vf}, (2-9b) where Pa = \E oa Pa = lpa -l Pa = p Or pa = or

12 The quantities V, V -, W and W are defined in Section 2.2. 2.4 The Coupled Equations The means of Eqs.(2-6") can now be taken and the closure relations (2-9) substituted. Observe that in (2-9). either Pa = 0 or p a =0. Thus, z 2 k1 2 Vb -)-i(k -k )+ k 1 2 k2 +k2 1 2 + ~ A(O)IV+ = k k - A(pb)(V ++V +++ -) - 4 ab 2 2A0~ i 2 L2k2b 2k1 a k1 ]w (2-10a) [. + az 2 k1 v2 k2+k2 - F)+i (k1-k2) 1 2A0)]IV = k2 12 4 4 ab (V++V- +W++W ) k2 - 2~- - 2k - 21V ZAO4 W+ 2k b 2 + -(0)]IW (2-10b). 2 v2 k2-k2 az - + b-)-_i(k +k2)+ 2 A() ]W+ = a 2k k21 2 4 i2 = [-V2 2 k + -A0]+[ji 2 - + () IV [~'v A(O 1 (2-10c) and [a+ i a az 2 k1 2 + b )+i(k +k )1 2 k2k2 1 2A(0) ]W- = 4

13 k2 2 i 2 1 2 - V A(O)V -[ - + ()]v, (2-10d) 2ka 2k b ' where V = <E (k1,Z,p )E (k 2Z') V = <E (klZpa)E (k2'zpb)> + + -* W = <E (k,za)E (k2,z, pb)> 11zaL a 2 & W = <E (k1 z )Ea (k)E zPb Eqs.(2-10) can be further simplified by substituting in the expressions 2 = (a+Pb) = xx+yy and Q = ea-Pb = xx+yy. Thus V2 12 2 Thus V2 V -V +V'V+V2 a 4 2 1 2 2 Vb = -V -V *V+V, b 4 2 2 where V' = + a ax'2 ay,2 2 2 ax ay = 2 +ay2 The physical structure of the random medium in this work is assumed to be a flat slab with plane wave normal incidences. It is also assumed to be statistically homogeneous in the x-y plane. Hence V2 (2-11) V = 0O which implies

14 2 = 2 = 2 a b that is, V, V, W+ and W depend only on the horizontal distance p = pa-_. Thus Eqs.(2-10) can be re-written as +kl-k2 2 az 2k k 1 2 2+-A(k2 4 )2a(0-i(k -k2)IV+ k k k2 -2 - [ -V2 + V2 4 2 1 k2 - A(0)]W = = S (2-12a) k -k az 2kk 2 2 - 4 A(0)+i(kI-k2)V = 41 k k k2 k =- + - + - 2 1 -i -2 2 k A(p)(V +V +W +W ) - [-V - -A(O)]W++[V + A(O)IW = 4 2k 4 2k2 4 1 2 = s (2-12b) k +k a kl+k22 [ - -i v 1 2 + -Ak2)W = + (0)-i(kl+k2)]W = 2i 2 = [2k 2 -2 2 2 + i - 1 4 2k 4 - = R (2-12c) and k +kv az 2k k 1 2 12 2 - 4 (O)+i(kl+k2) ]W = k2 iV2 - 10 + i 2 2k1 2k2 k2 + -A(0)]V =

15 = R. (2-12d) The above equations are a set of mutually coupled partial differential equations. The left side of each equation is a parabolic operator with constant coefficients whereas the right hand side will be regarded as the source term for the equation. 2.5 Boundary Conditions Assume a flat random slab in free space with a symmetric pair of incident waves E- and E, as illustrated by Fig.2-1. The field E travels in the positive z direction and illuminates the slab from the underside (z=0), whereas E propagates along the negative z direction and hits the upperside of the slab (z=L). The boundary values of these waves E (z=0) and E (z=L) are known quantities. From the definitions of V, V, W and W, we have + + +* V = <E E > V = <E E > +W W = <E E > W = <E E >. When E (L) = 0 or <E E >l=L = V (L) = 0, z-L W (L) = W (L) = 0, since there is only one wave propagating in the positive direction at z=L, because no scattering exists outside the slab. Similarly, when E (0)=0 or <EE+ > =V (0)=0, = > Z=0=V(0)=0, W (0) = W (0) = 0. Thus the boundary conditions at z=0 and z=L are, if V and V being considered as primary variables and W and W secondary or induced variables,

16 (at z=0) W+ () IV+()=O = W (O) IV+(0)=0 (2-13a) = w = (2-13a) and (at z=L) W (L) I-()=O = W (L)IV-() = 0 (2-13b) As for the boundary condition at P=4j -A -4 a, since the electromagnetic wave is propagating inside the random slab, the degree of coherence decays with increased horizontal distances between the two points. At, o, the medium correlation falls to zero, and therefore the boundary conditions inside the slab are V (p4) =V (p-*)=W (pc )=W (p >)=O, (2-14) and the derivatives of these quantities are VV (p-o)==VV (p-o')=VW (p-o+)=VW (p-o+)=0. (2-15) Equations (2-12) together with the boundary conditions(2-13) through (2-15) define a single valued solutions if the boundary values V+ (z=0) and V (z=L) are specified. 2.6 Summary of Chapter II 1. Medium structure: the random medium occupies the volume defined by O < z < L, -B< x,y <. Medium is continuous and is fluctuating slowly when compared with the incident field quantities. 2. Statistical properties of the medium: e = c(l+e) = 1+e <e> = e _ 1 <e> = 0 -,(z'p )^(z" " )> = 6(z'-z")A(p) with p = lc'-p"| A(p) is a Gaussian function with a zero mean. 3. Incident wave: plane wave, normal incidence.

17 4. Scattering by the medium: the size of irregularities is much larger than the wavelength, so that scattering is concentrated in the forward and backward directions. 5. To derive the closed set of equations, the Furutsu-Novikov relation is used.

18 Free Space incidence at z=L output at z=L L forward scattering Random Medium 0 backward scattering incidence at z= incidence at z=O output at z=O Free Space Figure 2-1. Structure of the Slab

CHAPTER III. THE ASYMPTOTIC SOLUTIONS 3.1 Basic Approach In this and the following chapter a set of approximate solutions for Eqs.(2-12) are obtained by using a combination of Green's function techniques, asymptotic integration methods and matrix manipulations. Since the quantities of interest are <E E > and <E E > or V and V W and W- have been eliminated in the process. The solution is given in the form of a 'scattering matrix' with V (z= 0) and V (z=L) as input and V (z=L) and V (z 0) as the output, similar to the scattering matrix of a two-port network V (O,p) - V+(O,) V (L,o) V(L,p) = S(kk-L,_Fl. By setting V (L,p) = 0, this solution yields the correlation functions of backward and foreward scattered fields at z=0 and z=L,respectively, as functions of the correlation function of the incident field at z= 0. 3.2 The Green's Functions Define a set of Green's functions G + G-, H + and H- for Eqs.(2-12) 2 2 k -k k +k [- -i V +i(kl-k2)- A()]G(z,z',,') 2k[k- 2 4 = 6(z-z')6(s-e') (3-la) 2 2 k -k k +k 1 2. 2 1 2 [- + V -i(k -k) >- A0)]G (zz'*,,') = 2kk 12 12 19

20 = 6(z-z')6(p-p') k +k k-k2 [k- +j1 2V2+i(k +k ). '2A(O)]H+(zz'tprp) = az 2k1k2 12 4 - 6(z-z')6(p-p') (3-1b) (3-ic) and.k1+k 2 2 [a V-i(k +k ). az 2k k 1 2 1 2 k2k2 + 1 2A(O)]H-(Zz1z'Pp) = = 6(z-z')6(p-p,). (3-id) Also define the Fouri~er transform pair f~s =(2r).f t()exp~i[K (X-X')+Kc (y-y')]}I dKc t(K) = ff(p)exp{-i[K c(X-X')+K (y-y'))} dp. Transforming Eqs.(3-1), we get a k -k2 [-z +i21k2 - 6(z-zt ) a k -k2 8az 2kk 21k2 K2+i (k -k )1 2 k2k2 1 2A(O)]G+(Z,z',,) = 4 (3-2a) K2_-i(k -k )+ 1 2 k2k2 4 ~ - 6(z-z' ) (3-2b).k1+k2 [L -i 1 az 2k k1 K2 +i (k +k ) - 1 2 k2k2 1 2A(O):Ifi+(tzr) = 4 - 6(z-z' ) (3-2c) and

21 k +k k -k [L +i K -i(kl+k2)+ 4-A(O]H (z,',K) = az k 12 4 2k1k = 6(z-z'), (3-2d) 2 2 2 where K = 2K = K +K x y The solutions to Eqs.(3-2) are easily found, and for z<z' k -k k +k ^= 2 2 1 2 G -exp{[-i K -i(kl-k2) 4 ()](z-z')} 2k k 2k1k2 2 2 k +k k -k A= 12 2 1 2 H = -exp{[i K -i(kl+k2)+ 4-A(0)](z-z')} 2k k 11 G = H = 0, and for z>z' 2 2 kl -k k +k = exp[ K +i(k-k2)- - A(0)](z-z')} 2k k 12 2 2 k +k k -k 12 2 1 2 H = exp{[-i K+i(k+k2)- - A ](z-z) 2k k1 A+ ^+ G = H = 0. The physical meaning of the solutions is that G and H are considered as 'forward propagating' quantities, therefore their values at a point z' are affected by a source located at z only if the source lies in the path of the ray, that is, only if z'>z. Similarly, G and H are 'backward propagating' quantities and their values at z' are affected by a source at z only if z> z'. It is also easy to obtain the Green's functions from their Fourier transforms. For z<z',

22d ik 1k2 G (zfz',pp,p) =1 1 2 4r(k -k )(z-z') 12 2-~1k2 ________ *exp{[-i(k - k2 ) 4 A(O)](z-z')+i 12} (-a 1 42(k - k2)(z-z') ik 1k2 G(z,z',p,p')= 2ir(k +k )(z-z') 1 2 k2k2 kk 2.1rf~(k +-k )2 i 12 }0 (-z(3-3b) expLl\1 22(k +-k2)(z-z') ik1 k2 G (z,z',p,p') =1 2ir(k1-k2 )(z-z') k2k2 kk 2 1 2__ __ ___2.exp{[i(k -k2 )- }01z-' (3-3c) 1 22(k - k2)(z-z') H+ (z,z',e~p') 12 zr~~,)=0 3.3 Green'2Function.Pr)pertie Thefolowig pope 22e 2r sflfreautntechrnefntos 1. From the inverse Fourier transform formula,

23 f(p) = 1 ftOK)expji[K (X-X')+Kc (y-y')]}dKc (2vr)2 xy and the expressions of G's and H's, the limiting values of the Green's functions are lrn G (zFz',p,p') = 5p, 34a +-Z lrn H (z,,z',p,p') = -6(p-p') (3-4b) lrn G (z,z',p,p') = 5(p-p') (3-4c) Z-+Z' + lrn H (z,z',p,p') = 6(p-p'). (3-4d) Z-+Z' + 2. Since all G's and H's are functions of the quantity (z-z'), aF =~ F. (3-5) az z 3.In their non-zero regions and for z~tz', V'2G+(z',zrp',,) 2k k k2k2 =- 1 2 [-a+i(k _k 1 12A(O)]G+(z',zp',p) (3-6a) kk az 1 2 4 -V'2 G (z',z~p',2p)= 2k k k2+k2 = i 12[- - -i(k -k )+ A (O)]G_(z',z,p',p) (3-6b) kk az 1 2 4 V12H(,ZP,) 2k k k-_k = - 1 2 [-a +i(k +k 1 2A(O)]H+(z',z,p',p) (3-6c) k +k az 1 2 4 1 2

24 and V' H (z',z,p',p) = 2k k k2 -k2 =-i [- a k-ikl+k2)+ i -i(k +k)+ 4 A(O)]H (z',z,p',~). (3-6d) k1+k2 3.4 Integration of the Equations With the help of the Green's functions derived in Section 3.1, the set of differential equations (2-12) for the coherence functions are easily solved to yield a set of integral solutions. From Eqs.(2-12) and (3-1), with the primed and unprimed coordinates interchanged, we have (V+G+ k -k (V G )+i 1 2V*(G+V'V+-V+V'G+) = 2klk2 1 2 ++ + = G S +V 6(z'-z)6(p'-p) (3-7a) k -k az'(V G )-i V'-(G V'V -VV'G) = 2klk2 = G S +V 6(z'-z)6(p'-p) (3-7b) a k +k,(W H )-i V'*(H V'W -W V'H ) = az' 2klk ++ + = H R +W 6(z'-z)6(p'-p) (3-7c) and a -- kl+k2 az'(W H )+i V'-*(H V'W -WV'H ) = 2k=k = H R +W 6(z'-z)6(p'-p). (3-7d)

25 Next define two volumes V and V2 as 1 = 0 0 < z' < z; -~ < x',y' < } V = { z < z' < L; - X << co Note that the horizontal plane z=z' is excluded from both volumes. Also, all the z-values are bounded by 0 < z < L and 0 < z' < L. Integrating Eq.(3-7a) over V1, using the boundary values Eqs.(2-14) and (2-15), the limiting values Eq.(3-4a), in view of the fact that z&z' in volume V we get V+ (zp) =-Jdp'V'(0,p')G'(O,z,p',p) - Z co -fdz' fdp' S+(z',rp')G (z',z,,p p). (3-8a) O -co Similarly, V (z,p) = fdp'V (L,p')G (L,z,p',p) - L - fdz' fdp' S-(z',p')G ( zzp',p) (3-8b) z -C W+(z,p)= -fdp'W+(O,P')H+(Oz,2',P ) - -co fdz' faa' R (z',2')H (Z',Z P, ') (3-8c) O -ap W-(z,p) = fdp'W(L,p')H (Lzp'r2) - L co - fdz' dap' R-(z',p')H-(z',z,t',p), (3-8d) Z -cm

26 where in the process the second Green's Theorem has been used. In the integral solutions Eqs.(3-8), all the source terms S, S - R and R contain V V and V W - and can be treated as the following example: -00 I = fd'G+(z',z, P',p)V W(, )= co oo = fdp'V'.(G+V'W+-W+V'G+) + fdp'W+V'2G = -00 -00 = Jdp'W+V' 2G -00 By substituting Eq.(3-6a) into the above expression, we get 2k k k2+k2 I=-i [- +i(k -k)- 4A(O)]. k-k az 1 2 4 -k2 00 jdp'W (z',p')G (z'2z,p'',), and thus Eqs.(3-8) can be reduced to the general form X(z,p)=+ f dp'X(z',' )F(z',z,p'p)-dz' jdp'f(z' p')F(z' z,p',P), -oo -oo where X and F designate one of the correlation functions and the corresponding Green's function, respectively. From Eqs.(3-8a) and (3-6a), v (z,) = -f d'V+(O,~')G+(O,z,~',~ ) - -co

27 - 12 E+ + - + (ttp IFdz' f dp'A(p')[ +V +W +W ]G(z,2')+ 4 0 +-:A(0) dz'CA d'rk2W +k2W J-Z~~2 1]G (z',ztp',P)0 - k2k2 1 1 2 -+i(k -k )- A(0)]. k1-k2 az 1 2 4 122 0 - from Eqs.(3-8b) and (3-6b), V (z,p)= fdp'V (Lfp')G (L..z,p',p) + -co k k L 0 +f2 dz' fdp'A(p')[V +V +W +W]IG(z',ztp',p) - 4 z - 1 L 2-,-,, — A(0)fz Wd'[ +k WI z,~'p fdz ' 1 2WGz,,) z - k2k2 k k azk +k4 1 2 L *fdz' fdp'[k2w -k1WIFG (z',Z,p',P), 3-b z - from Eqs.(3-8c) and ~(3-6c), W~ (z, ) =fdp'W+ (0,P')H+ (0,zp',rp)-co 1 Z +V 0 -c

28 2 2 0 — k and from Eqs.(3-8d) and (3-6d), CO 1 L 2+ 2 — - A(0)fdz Jdp '[k1V -k2V] H(z',z, ', p) + 2 2 L c *Jdz' dp [kV +k V ]H (z'zr^P' P). (3-9d) z -o In Eqs.(3-9), all V's and W's are evaluated at ( z', '), unless explicitly noted otherwise. 3.5 Stationary Phase Method The well-known Kelven's stationary phase method is widely used to find the asymptotic values of integrations with the form b 1(z, ) = fdxP(x)exp[i7f(x)]; z > 0. a It states that if — OO i. f(z) and O(z) are analytic where z = x,+iy, ii. f(z) is real for real values of z, and iii. in the integration region f(x) has a finite number of 'stationary' points where iii. in the integration region f(x) has a finite number of 'stationary' points where f(x) = 0,

29 then b I(7) = fdx?(x)exp[iYf(x)] = a N = in 271f,,(an) l(an)exp[i7f(a )~i4]}+0(1/y). n=l ( n 4 The a (n= 1,2,...,N) are the stationary points, and the sign in the exponential of each term is the same as that of the corresponding second order derivative of f(x) at that point. In the set of Green's functions G + G H + and H. the x and y dependences are of the form (see Eqs.3-3) exp{ia [ (x-x')2+(y-y)2] } -ZZ where a is a constant. The stationary phase method can be applied to carry out the x and y integrations in Eqs.(3-9) so that a set of asymptotic solutions are obtained. Notice that the only stationary points of functions (x-x') and (y-y') are x=x' and y=y', respectively, and the second order derivatives of these functions are constants 2>0. From Eq.(3-9a) k2+k2 k +k V (z,) V(0,)exp{[- 1 2A(0) +i(k -k )](z-O)} + 4 12 k k2 z + -A(p)fdz'[V (z',p)+V (z',~)+W (z',p)+W (z',p)] 4 0 k2+k 2 *exp{[- 1 A(O) +i(k -k )](z-z')} - - 1A(O)dz[k2W (z ',)+kW (z )] 0 2 1

30 k2 +k2 *exp{[- 1 A(0)+i(k 1-k2)~-vl 31a from Eq. (3-9b) k2k2 V(z,p) V'(L,p)exp{~[ 142 AO k1 k2)]'(z-L)} + k~k L 4 z k2+k2 12 *exp{ A()-+k-2)(2z) 1A(O)fdzY[k2W (ztp)+k2W (z',p)]. z k2+k2 121 w+ ('P)w +O,*exp{[- 1A(O)-+i(k -k )(-~ 4 1 2](-)}(-lb k2-k2 (zp *O~exp{[- 12A(O) +i(k +k2)](z-zO)}r +31c 221 k2-k2 W*zp -Lpexp{j[-1 A(0)+i46(k1 +k2)](z-z)}, (3-c 4 12 z

31 2 2 k -k ~exp{ [ 4- (0)-i(k+k2) ] (z-z ). (3-10d) A set of mutually coupled differential equations are derived from Eqs.(3-10) by direct differentiation. These are +k2 k k k k a 12 12 + 1 2 -v = A(O)+i(kl-k)+ A(p)]V + - A(p)V + a 4 1 2 4 4 4 4 2 2 kk k kk k 1 2 2 2 A - [ = -A - 7 A (O)]W [ a[(10 A( )- -A )]W ()] (3-lla) 4 4 2 2 kk k+k k 1 2 + 1 2 1 - -V =- -— AA(p)V +[ 14A(01(O)-i(k -k - A(p)]V + az4 4 14 kk k k k k 12 1 + k +[- cop(r)+ be-(e)]W++[- -12)AC) re -A(OlsW (3-llb) 4 4 2 2 k2 k k-k a + 2 + 1 2 + z = -A(0)V - A(O)V -[ — 2A(0)-i(k+k) ]W+ (3-c) k k k2-k |-W~ = -A(O)V - A(0)V [ A(O)-i(kk)]W (3-11d) where all variables V, V- W+ and W- evaluated at (z,p). A comparison between Eqs.(3-ll) and (2-12) reveals that the application of the stationary phase method results in the elimination of all the V terms contained in the sources of the original equations. This result is consistent with the assumed normal incidence and the 'small angle' approximation.

32 Eqs.(3-11) are next written in a more concise matrix form a v v az = B + (3-12) w w Since the solution of Eq.(3-12) will be used to calculate pulse propagation problems where only the values of the coherence functions at _pO are required, it is sufficient to solve Eq.(3-12) for_ p=O. The coefficient matrix B of the equation then has the elements k k -k -k b =b 1 1 2A +i(k -k2) 11 22 4 12 klk2 12 12 -b21 4 kA 2 b = -b (k-k2 k A 14 = 123 2 31= -b = 4 2 b32 b41 4 k2-k2 1 2 33 =-b44 — A+i(kl+k2 b34 b43 = where A = A(0).

33 By introducing new variables u = (k + k2) -A and v = (k -k2)A, the matrix B becomes 2 2 u2 +3v2 6 +iv 2 2 u -V 16 2 (u-v) 16 (u+v)2 16 2 2 u -v 16 u +3v2 16 —V (u+v)2 16 (u-v)2 16 u-v — v 8 u+v 8 v -u(-v -i) 4 U+V 8 U-V u-v - -~v 8 0 (3-13) 0 u(- -i) 4 The set of equations for the two-frequency mutual coherence functions, Eqs.(3-12), with coefficient matrix B as given in Eq.(3-13) is solved in the next chapter.

CHAPTER IV. SOLUTIONS FOR THE SLAB In this chapter Eq.(3-12) is solved for the flat random slab shown in Fig.2-1. The solutions are given in the form of a scattering matrix defined at the begining of Chapter III. As an intermediate step, a 'transmitting matrix' T (L), defined as V (L) V (0) + (L) +(0) (4-1) W (L) W (0) is derived first. Appropriate boundary conditions are applied to Eq.(4-1) to eliminate quantities Wk+ and WV, and the scattering matrix derived. 4.1 Laplace Transform Method For a matrix differential equation a az X(z) = B X(z), where B is a constant matrix, the Laplace transform is sX(s) - X0 = B X(s) or 34

35 - — 2 X(s) = (SI - B) Xo. Hence X(z) = L1 {(s - B)1} X 0 where X0 is the boundary value of X(z) at some plane z=z0. Specifically, X(L) = L1 {(si - B)- }z= X(0), and z=L X(0) = L1 {(SI - B) -} X(L). 4.2 The Matrix (sl - B) — 1 The matrix (sI - B) can be calculated following the standard procedure, although it is quite tedious in the present problem. 4.2.1 Eigen values of (sI - B ) Eigen values of matrix (sI- B) are found from the characteristic equation 2 s 1 2 2 2 2 +[u +v +i5u +3v2)]- + -[u2v2- u +3v )+ +7v2 2 = 0. (1- 16 )+ 8 uv = 0. A A It is easily solved to yield s1 = (a+ip)/A s2 =-(a+i3)/A s3 = (S+i7)/A s4 =-(S+i7)/A

36 where U2 a = -vv/[-l+'() 2+1] v2 u 2 +3v2 S = -v -[1+V( v)+l] 1 U2+3v2 7 = -v/[l+( v )2+8]v 4.2.2 The inverse matrix of (sI- B) -1 The elements of matrix (sI- B) P.. can be expressed as -J, where P.. is a A polynomial of no higher than the third order, and A= sI-B| is a fourth order polynomial. That is, -P. (sI-B)-. = ] 13 A where A = s4 [ + [ v2+i (5u2+3v2) +3 2 S - A 1 1 2 2 + -[u v (1 -A4 2 2 2 2 u +3v 2 u2+7v2 2 16 ) +i 8 u v] 2 2 2 2 = [s - -(a+i) ][ s - -( +i?) J. A A (4-2) The expressions for P..'s are Ui

37 3 u2+3v2 s2 2 u2-v2 2 P (s) = s -( -iv)- +(u - v +i 11 16 A 64 2A2 U2 2 2 - Iu +11v A3 16 u2+7v2 +iv( 64 -1)]I P2 (S) = -P 11(-S) 22 11 1216A E 2 v 2 U s+ (I+, A P21 (S) = P1 (S P()=(u-v)v 2 u+v v UV u-3v 13 8A A 8 A P2 (s) =-P 13(-S) p I=-(u+v)v 2 u-v(y -i UV.1+iu83)] p14 S) 8A [s A A i~ p23 (S) -P 14(-S) U 2uv 2 UV U2-U+2 P31(s 16A {(u-vs A ~ v -i(u-v)]s - A (u-v- u4uv2 P42 (S) = P3 (S u~v2 u-v v2 P32(S) y- " {(u+v)s -[ —i(u+v)]s + 32 16 -jA4 (u+v+i u4+uv+2 P41 (S) = -P 32 (-S) P33 (S) = s uv~(- -i) s2+ v 2-2 ~~(+u-v+ u2+3v 8 )

38 uv 2 2 2 2 UV u2+5v2 3u+v 2 3 [ +iv(l- 64 A )] P 44(s) = -P33(-s) 2 2 P34 () = 2v (s-i A) 34 64A2 A 64A P43(s) -P34(-s). 4.3 The General Solution It is now a straight forward procedure to inverse transform the elements of (sT - -B-) to obtain the general solution of the matrix equation (3-12). Let the matrix L { (sT - B ) } be denoted by T (z), with P (s) P.. t..(z) = L 1{-i —} = Z Res{-,sk}exp(skz), (4-3) A(s) k A where sk's are poles of A(s). 4.3.1 Pole type analysis P.. Generally, -- has four simple poles located at A a+ and + n A A so that the residues are P.(s) Pi (sk Res{-2 -—, Sk} = (4-4) A(S) A' (Sk) where A'(sk) is the derivative of A evaluated at sk. However, higher order poles could exist at a+iq3 = 0 and +ir2 = 0, as can be seen from Eq.(4-2).

39 4.3.2 Special cases 1. a+i3 = 0 From the expressions for a and 3, this can occur only if u = 0. There are then P.. three simple poles located at s = 0 and s = ~+(;+i )/A for every J. These expressions A are 4 2 v 3v 4 s2+ -is+ - S+[s + (- +i-)s+ 2 2 P22 6A64A V V A2 v2 s[s + (l+iv-)] PA 2 2 v v s[s + (l+i v)] 2 8 2 2 V VV P13 A A -i) V2 s [ 2 ( +i3 2 8

40 2 V + V(V - P24 8A A 8 A 2v2 A2 8 p1 p1 p 23 A A 2 V ES- -(- -i)]I p31 16A__ __ _A_ __4 _ 16 v2 A 3 S[s + +i) 2 v[s+ - (- -i)]I 42 - 16A A 4______ 2 3 S - +(1+i~v)]I 32 _ 42 A6 A 41 _ 31 2 S + 1 V2 3 T4 +i~v) p 33 -A 2 S~s + +i)I

41 p44 _ 33 A A p 34_ A 4 V 6 4A 2 2 S[s + +1 )I p43 _ 34 A A6 with v2 8 1]' 7 = - V[41+VI/Q)2] The corresponding t.(z)'s calculated from Eqs.(4-3) and (4-4) are v2 3 v 2 3v ~+ +(.+i-v -)C S 64 8V 64)C(6 ') v t1 (z) = 3.+i-3v 8 (4-5a) t22 (Z) = t11 (-Z) (4-5b) 2 ~ +7 64 3.C)+ 6 t 12 (Z) (4-5c) 1+i-3v 8 t2 (z) = t1 (4-Sd (4-5d)

42 V(V -i) (l-C)+ 8+7 S t 3(z) = (4-5e) 13 t24(z) = t13(-z) (4-5f) t (z) = t (z) (4-5g) t23(z) = t24(z) (4-5h) 2 v 64 t (z) = 1- (1+C) (4-5i) 333 1+i-v t44(z) = t33(z) (4-5j) 8(8- 2)(1-C)+ 1 S t 31(z) = - (4-5k) 31 3 1+i-v t 41(z) = t31(z) (4-51) t 32(z) = t31 (-z) (4-5m) t 42(z) = t32(z) (4-5n) 2 V 64 t34 (z) = - (1-C) (4-5o) 3 1+i-v t43(z) = t34(z). (4-Sp) In the above quantities C(z) = cosh[(Q+in) A]

43 S(z) = sinh[(~+in) z ]. 2. ~-i-jn = 0 From the expressions of ~ and r7, this implies v = 0. The results then are p11 A A2 1 S 1 S 2 U 16As 2 2 u 16As 2 p12_ A As2 p21=A As 2 13 =3 14 = 23 = 24 34 = 43= p31 A p32 A p41 A A.uJ. 1 - A = u l 1 16 s S1 A 51A =.u 1 1 16 s Si 51A

44 P33 _ 1 AA 44 1 A S+i u The corresponding t 'is are 11 Z 2 = 1- z 2 t22(Z) = 1+ 2 t12(Z) = 16 2 t21(Z) = - 16A t13 = 14 = 23 = 24 = 34 = 43= (4-6a) (4-6b) (4-6c) (4-6d) (4-6e) (4-6ft) (4-6g) (4-6h) (4-6i) t 32(Z) 41 Z 42 Z = i j -2[l-exp(i K z)] = - i # [1-exp(iu z)] = i - [1-exp(-i~ z)] 16 A U t33 (z) = exp(i-i z) (4-6j)

45 t4(z) = exp(-iA (4-6k) 3. a+iP = -+ir = 0 This is a trivial case since the condition requires that u=v=0 or k = k2 =0, and solutions for this case are readily found from Eqs.(2-5) t. = 1 (i=j), and = 0 (ij). 13 This solution can also be obtained from from Eqs.(4-6) by allowing u go to zero. 4.3.3 Solutions for the general case (u 0, v 0) P.. Generally, when u - 0 and v t 0, there are four simple poles for each -. From A the expressions for A and P..'s, derived in Section 4.2.2, the t..(z)'s can be obtained by substituting Eq.(4-4) into (4-3), giving 2.U 'Z71V 1 2 64 2 t22(z) = — + S +(1+iv - )C + 22 I 641 1 8 64 2 l+i3 a+if 1+i-v 2 2 2 2 u+9v2 u +9v2 16 +1V 128 1 + S2 ] (4-7a) a i; - t11(z) = t22(-z) (4-7b)

46 1 2 t (Z) = 12.364 3~ ~ 2 i-v a+ij3 2 1 642 2 2 2 2_ 2 u-.u-3v 16 +1 128 - + s2 t 21 (Z) = t 12 (-Z) (4-7c) (4-7d) - (- -i) 88 t23 (Z) = 3 1+i-v 8 t1 (z) = t23 (-Z) n H-i) t24 (Z) = 3 1+i-v 8v (C1-C) v( U + Vs ) 8 a+ig3 1 +it1 (C -C) X ( Us V 1 2 8 ca+ig 1 ~i (4-7e) (4-7f ) (4-7g) t 13 (Z) = t 24 (-Z) 1.u+2v 16 2 U+V-1 t41(Z) = {[Xv. -i(u+v)](-C +C + us1 + 4 3 4 1C2) a+i3 1+i-~v 2u +uv+3v2 U+V+1 8 ++i (4-7h) (4-7i) t32 (Z) = t4 (-Z) (-j (4-7j)

47 1. u+2v 16 24 t (z) = [... +i(u-v)](-C +C us42.3 4 1 2 a+ij S 1+ i-v.u2-2uv+3v2 U-V+L 8 + Vs 2} (4-7/k) t3 (Z) = t (-z) (4-71) 2 V 64 1 s2 t (z) = [-C +c +iu(- + )] (4-7m) 43 +3 1 2 a+iO3 ~ +in t3 (z) = t (-Z) (4-7n) 5 5v2 2 ~ v+i(- -41) t4 (z) = 3[(1+i-v- 64L c + aO us + 1+i-v 2 i IU ( S,)(4-7o) 64 2 ~+i7) t33 (Z) = t44 (-z), (4-7p) where in t..(z)'s ii C (z) = cosh[(a+io3)-) 1 A S (z) = sinh[(az+ig3) ] 1 A C2(z) = cosh[(~+i'?)z] 2 A

48 S2(z) = sinh[(S+i7A) ] U2.V =U2 V a = -2[-l+/() 2+1] 2v2 2 2 = + [ U 3v)2+1] = 8v 1 U2+3V2 = -v/~[1+(/( ) +1] 4.3.4 The matrix solution From the solution given by the Laplace transform method, the equation x(z) = B x(z) has the solution x(z2) = L- {(sI-B) }z - x(z 2. The solutions for a flat slab with thickness L are thus V (z=L) V (z=O) + (z=L) (z=) (4-8a) W (z=L) W (z=O) W (z=L) w (z=0) and

49 (z=O) V (z=L) V (z=0) V (z=L) +(z=O) T(-L) ) | (4-8b) W (z=O) W (z=L) W (z=0) W (z=L) where T(z) = L1 {(sI-B) -1 4.4 Solution for a Single Layer For the convenience of the readers, definitions for quantities V+,,, W+ and Ware again listed here. They are V (z,kl,k2) = <E (z,kl)E (z,k2)> V (z,kl,k2) = <E (z,kl)E (z,k2)> W (z,kl,k2) = <E (z,kl)E (z,k2)> W-(z,kl,k2) = <E-,kE(z,k2)>. For a single layer of random medium with incident waves E (z=0) and E (z=L), the boundary conditions are stated in Eqs.(2-13a) and (2-13b). The solutions of the twofrequency mutual coherence functions of a flat slab of a continuous random medium for the forward and backward scattered fields due to a single incident field are derived next. The physical layout is shown in Fig.(4-1). From Eq.(4-8b): V (z=O) = tl(-L)V+(z=L) + t12(-L)V (z=L) + + t3 (-L)W+(z=L) + t4 (-L)W (z=L) (4-9a)

50 V(z=Q) = t21-L)V+ (z=L) + t22(-L)V (z=L) + 21 +2 + t 2 (-L)W (z=L) + t 2 (-L)W (z=L). (4-9b) Applying boundary condition Eq. (2-13b) to Eqs.('4-9) gives V~ (z=O ) V-(zL)=O =tll(-L)V (z=L) V-(z=O) IV~(zL)o=O t21 (-L)V+(z=L). Therefore V+ (z=L) IV-(zL)=o= V+ (z=O) t1t (-L) + 21 (-)+ V (z=O) V-(zL)O = t 21 (-L)V (z=L) = V (z=O). t1 (-L) Assign t1 (-L) -t22 (L)(4l) =t21 (-L) -t12 (L)(4lb t1 (-L) t22 (L) Then for the transmitted field at z =L V+(z=L) = <E + (Lrk )E + (Lrk2)> =I' (urv)V+(z=O) (4-11a) and for the backscattered field at z =0 V-(z=O) =E (O,k )E - (Olk )> = r -(u,v)V+(z=O), (4-11b) where u= (k 1 +k2)A

51 v = (k -k2)A A = A(p=O). The functions F and F are elements s21 and sl of the scattering matrix S, respectively. The complete solutions for the two-frequency coherence functions are: 1. for u=0 and vtO r (o,v) = r (O,v) = 3 1+i-v 2 2 v 3 v 3v S+ij -+i V >c - < - i> S 64 8 - 64 16 v 2 V (+ic) 64 (1-C) + 16,,. (4-12a) (4-12b) 2 64 +(l+i - )c - 8 64 3v - +i S (16 -i V C = cosh[(Q+in) L] A L S = sinh[((+i7) L] A = v- [-1+ V/(jv)2+1] 2 8? = -v /V2 [1+ V(iv) 2+1] 2. for v=0 (including the case u=v=0) 1 r+ (u,0) = u2L 1 + 16 16A (4-12c)

Po"2 F (u,O0) = u2L 2 1 16A (4-12d) 3. u:;O and v:;O F+(urv) = Q1~v (1+i 3v) (4-12e) f,(u, V) = 1 2 _[-(C -C ) Q(ulv) 64 2 2 a+ij3 2 2 u2 3v2 16 1v128 ~s2j ) (4-12f ) where Q (u, v) = 2 lZ64V 1 a+io3 2. 83 2 VI + -64~ C2+ 16 + iv(12 -1i) ~ +i71 c (u,v) = Loh(~i~i s (u~v)c=sih[(a+io3)-] 1 A S (u,v).= csih[(a+i77)-L] 1 'A C2(u,v) =sih(in ] 2 A ca = v /LlJV/V~)2+1i 22 2

53 = V2 [-+(U 2+3V2+1 = - ^ [-l+/(" 83 )2+1] 28v = - /[1+v( U2+3v22 8v 2 1] It is seen that for a single layer problem, only t22(z) and t12(z) are needed since W + and W- are eliminated by application of the boundary conditions. However, for multi-laver problems, all t..'s will be needed to derive the solutions. It is noted that the above-derived solutions satisfy the energy conservation law, represented by r+(u,O) + r-(u,0) - 1, since when v=(k -k2)A(O)=0, the case corresponds to a continuous monochromatic plane wave input and thus r and F- represent the power densities of the forward and backward scattered fields, respectively. 4.5 Conventional PEM Result For later comparisons, the coherence function of the forward propagating field calculated by conventional parabolic equation method is derived in this section. If the backscattering effects are ignored, the equation for the coherence function of the field is a. 12. i i2z +(k k )V2+ [k)A(O)-2k k2A(e)+k A(0)]+2(k -k2)} *<E+ (z,kl^,P)E (zk2 )> = 0, (4-13) where \PI = _Pil-p21. Note that Eq.(4-13) is the two-frequency equivalence of Eq.(1. 14) in Reference [3]. The plane wave solution of (4-13) atp =0 is

54 + (k -k )2A V (z=L) = V (z=O) exp[- 1 - k2 L + i(kl-k2)L], 8 or, in the u-v plane e(u,v) = exp [- ( — iv)]. (4-14) pern. A 8 4.6 Summary of Chapter IV In this chapter the general solutions linking the quantities V, V-. W, W- at two planes z=0 and z=L are derived by applying Laplace transform method and are presented in matrix form(Eqs.(4-8a) and (4-8b)). The expressions of t..'s for various u and v values are given by Eqs.(4-5) through (4-7). For a single layer and a single incident field the boundary conditions are applied to derive the two-frequency mutual coherence functions of the transmitted and backscattered fields as functions of the incident field. These results can be used to solve problems such as pulse propagation and backscattering. Since the solutions are given in the form of a transmitting matrix T (L), it is particularly convenient to extend the solutions to multi-layer problems. In such case the general solutions are V+ V (z=O) W W(z=0) W W (z=0)

55 Free Space transmitted pulse L forward scattering Random Medium i backward scattering incident pulse incident pulse backscattered pulse Free Space Figure 4-1. Physical Layout of the Problem

CHAPTER V. A CALCULATED EXAMPLE In this chapter the coherence functions derived in Chapter IV are employed to calculated the scattering effects of a random slab on a string of rectangular pulses. The resultant pulses are plotted together with the incident pulse. For computation and plotting, The University of Michigan Amdahl 5860 computer was employed. 5.1 Basic Formulae The quantity P(t! = <EE > is proportional to the magnitude of the Poynting vector of a field E(t) and thus is directly linked to the strength of a detected signal. Its value is given by * *co P(t) = <E(t)E (t)> = fdcwlfdc2 <E( ()E (2) )> e 1 2 -00 -00 where E(W) = - f dt E(t) ei. -CO Therefore P(t)'s of the scattered fields are P (z=L,t) = <E (L,t)E (L,t)> <'+ () ' +* i(W-W)t = rdw Sdfo E (Wl E (W2)>e( 1 2t = d m -i -o ) (5-la) = jdwldw2 V L(z=) e- (Wl-00 -00 = jfdlfdw2 rV(z=O) e( 2 -00 -00 and 56

57 P (z=O,t) = <E (O,t)E (O,t)> co co. (5-lb) = fdufd w r (z=0) e, 1 (5-b) -00 -c0 where V (z=0) is the two-frequency coherence function of the incident field at z= 0. There are two properties common to all coherence function expressions in Eqs.(4-12) and (4-14) that are particularly useful in the simplification of the calculations. These are F(u,v) = r(-u,v) (5-2a) and, (5-2b) r(u,v) = r (u,-v) 5.2 Fourier Series of the Pulses For a continuous string of rectangular pulses e- -i2rNt (nM < t < nM+l, n=O, ~1, ~2,...) E (z=O) = - 0 (otherwise) where M > 1 is the period of the string, and N > 1 is the carrier frequency, the frequency domain expression of the incident. field is (notice that the choice of the duration of a pulse as unity simplifies the calculation considerably without losing any generality since that is just a scaling in the time domain) sin(M -N)r E (=) -- exp[i(~ -N)T]6(c- 2n) ( M M n s ~-oo n T chn IM I e 6(c,-2,T(n +N)). -o n The coherence function of the incident field is then

58 (z=O) = <E: (z=0O,&1) ~t*zO,cA)2) - z ex~i2 5.12- +N)'6c -27T(- +N)) 2 ex~lM ~lTM )(2 M n - ( )26(w( N)6&2-2r(- +N)) + tr~ n n.n-J 1 coS inljT Sin. r ir {- e 65(ow-2tr(-! ni 2r( +) 7r~n n-1 + e1M6(wl2r(f +N))6c -2ir(n +N))} +.n.n-2 1 00 XIIRT SinflT0lT i2- n-2 + {e M68(12,(n +N))6 )-27r( — +N)) + v -w n n-2 + ei2 6(wl-2ir(~ +N)) -2tr(n +N))} + - (5-3) The Fourier series of the scattered power densities, F C, can now be derived by substituting Eq.(5-3) into (5-1), and applying relations (5-2). Thus, + 1 coS rl)2+ P S (t) =2:Z ( ) r'((i+ i;)2A, o) + ir n=-= n

59. n. n-m 2.1-2t sin-r sin — IT 2 m i- -mlr MA + - Re{ Z e M Z *. - r ((2+ )A, M)}, (5-4) IT m=l -c n n-m where A = 2rA(0) is a dimensionless quantity, X= is the carrier wavelength, and 1 - F (u,v) and r (u,v) are given by Eqs.(4-12) in Chapter IV. 5.3 The Scattered Pulses The coherence functions r (u,v) and F (u,v) derived in Chapter IV are substituted into Eq.(5-4) to calculate the Fourier coefficients of the forward and backward scattered pulses. Because of the limited computer time available and the large amount of -3 calculations required, the accuracy of each Fourier coefficient had to be set for only 10; the carrier frequency N was chosen as 10, and the period M has been chosen as 4. Those choices pose limits on both the thickness of the slab L and the horizontal correlation, A(0), of the medium. The parameters A and L for all the data sets are listed in Table 5-1. Table 5-1. Parameters Used in the Example Data Set No. 1 2 3 4 5 6 A(0)/X 0.02 0.02 0.01 0.10 0.002 0.01 L/X 2.5 5.5 11.0 2.5 10.0 10.0 The results are plotted in Figures (5-1) through (5-6). The incident pulse reconstructed from its Fourier series, i.e., Eq.(5-4) with r= 1, instead of the original waveform, is plotted in each case. The purpose is to show the influence of the accuracy and the limited number of harmonics used in the calculation, as well as the Gibb's

60 phenomenon. The computer program used in the calculating and plotting of the pulses is listed in the Appendix. In all the figures in this chapter, the incidence starts at t= 0 and ends at t= 1. 5.3.1 Transmitted pulses L i L The transmitted pulses are delayed by r=-= --- *- L which is to be expected and C IN A. can be seen on Figs.(5-1) through (5-6). After a relatively steep rise due to the effect of the arrival of the primary waveform, there is a 'turning point' from where the magnitude increases much more gently, which can be explained as the multiple scattering effect. The I L increase turns into an abrupt drop in the peak value at t= 1+N-.- when the back edge of the incident pulse passes the plane z=L. This drop is followed by a gentle decrease in magnitude. and again, this can be explained as the influence of the multiple scattering effects. These phenomena are most obvious in Fig.(5-4), representing the result of a severe disturbance (X = 0.1). 5.3.2 The returned or backscattered pulses The pulses plotted in Fig.(5-1) represent the scattering effects of a thin random slab. The backscattered pulse starts to rise at the moment the incidence hits the slab. It reaches the turning point at t 2r when the first order backscattered field due to the front edge of the incident pulse generated by scatterers located near the far edge of the slab (z = L) arrived at the plane z =0. From then on the pulse continues to rise, but at a much slower pace, particularly noticeable in Fig.(5-4), an exhibition of higher order scattering effects. The severe fluctuations of the returned pulse show that for large values of A(0) a higher accuracy than that used in the computer program (10 ) and more harmonics should be used in the computation since the exponentials in r are fluctuating very violently. The returned pulse reaches its maximun value when the back edge of the incident pulse enters the slab. The magnitude then falls gradually to zero.

61 Plots in Figs.(5-2) and (5-3) represent the scattering effects of relatively thick slabs. When L is large enough so that the incident pulse ends before the first order backscattered wave caused by the scatterers near z =L arrived at z=0, there is no flattening in the backscattered pulse's rise; the magnitude increases smoothly until the entire incident pulse passes the plane z= 0, then the gradual fall begins. Since in the computation the period M was chosen sufficiently large so that before the arrival of an incident pulse the response of the previous pulses has fallen to practically zero, the results can be considered to be the response of a single pulse instead of the periodic string of pulses. 5.3.3 The influence of A(O) The obvious effects of an increased A(0) are the increased magnitude of and the energy contained in the returned pulse whereas the same parameters will show some decreases for the transmitted pulse. In addition, the shapes of the pulses are also influenced by increased A(0). Plots in Fig.(5-5) represent a small A(O) case in which the returned pulse shows three distinctive regions: from t=O to t=1, a smooth rise; from t=1 (when the incident pulse ends) to t=2 (when the first order backscattered wave due to the front edge of the incident pulse caused by scatterers located near z=L arrives at z=0), a graduate drop; and a steeper drop till t=3 (when the first order backscattered wave due to the back edge of the incident pulse caused by the same scatterers arrives at the plane z=O0). For t>3, the returned pulse is practically zero. The transmitted pulse shape is almost the same as that of the incidence. Plots in Fig.(5-6) represent a slab of the identical thickness but with a more severe disturbance, A() =0.01. The returned pulse now is almost triangular in shape and the magnitude does not drop to zero until t=3.5. Note that the transmitted pulse is noticeably different from the incident pulse. This behavior can be explained as that for a very weak disturbance (small A(O)), the

62 first order scattering effect is predominant in the returned pulse and the zero'th order effect predominant in the transmitted pulse. As the disturbance of the medium increases in magnitude (larger A(0) values), the higher order or multiple scattering effects become more and more important. 5.3.4 Energy conservation property The energy contained in each pulse is calculated and printed under the plotted pulse in each figure and also listed in Table 5-2. As the numbers show, the total energy contained in the transmitted and the backscattered pulses is always equal to the incident energy. Since the medium is assumed to be lossless(E is real), the solutions obtained obey the energy conservation law. Table 5-2. Energy Relations of the Pulses Data Set No. 1 2 3 4 5 6 Transmitted Energy.663.476.477.289.827.501 Returned Energy.327.512.514.701.164.491 Total.990.988.991.990.991.992 Incident Energy.990.990.990.990.990.990 Relative Error.00% -.20% 10%. 00%.10%.20% The transmitted and backscattered energies are functions of A(0), L, the incident pulse period M and the carrier frequency N. A very good estimation suitable for cases that N> > 1 and M> > 1(which are almost always true in communications) is to use the monochromatic plane wave power relations e = e. (5-5a) N 2 1+(ir -) A(0)L C and

63 N 2 (t N)2 A(0)L e = - e., (5-5b) N 2 l+( r -) A(0)L where e. is the total incident energy et is the transmitted energy e is the returned energy c is the velocity of lighe in free space. The two coefficients in Eqs.(5-5) are the transmitting and the backscattering coherence functions for the special case v=0 as given by Eqs.(4-12c) and (4-12d), respectively. When Eqs.(5-5) are used to estimate the energy levels for the examples shown in Figs.(5-1) through (5-6), the relative error for the transmitted energy ranges from 0.10% to 0.98%, whereas the error for returned energy is between 0.52% and 1.91%, as shown in Table 5-3 (the numbers printed under the pulses are calculated by integrating the power magnitudes against time, thus representing the 'actual' energy contained in the pulses). It is seen that at least for parameters used in examples computed, Eqs.(5-5) give quite good approximation to the energy levels contained in the pulses. Table 5-3. Accuracy of Eqs. (5-5) as Energy Level Estimators Transmitted Energy Backscattered Energy Data Set No. 'Actual' Estimated Error 'Actual' Estimated Error 1.663.66957.98%.327.33042 1.04% 2.476.47947.72%.512.52053 1.64% 3.477.47947.52%.514.52053 1.25% 4.289.28840 -.21%.701.71160 1.49% 5.827.83515.98%.164.16485 0.52% 6.501.50328.45%.491.49672 1.15%

64 5.3.5 Results of the Traditional PEM In Figs.(5-7) through (5-9) the transmitted pulses calculated by the traditional parabolic equation method are plotted along with those calculated by the modified parabolic equation method formulated in this thesis. The influence of the random slab, which could A(O)L be roughly measured by the quantity —, is of increasing order in those plots. It is seen A(O)L that for small values of — 2-, the pulse shapes in both cases are quite similar to the A(0)L incident pulse. As the value of; — increases, the difference in maximum magnitudes becomes larger and larger and the shapes are more and more dissimilar. Notice that the magnitude calculated by the conventional PEM remain the same, due to the assumption that the backscattering is negligible. It can also be seen that the conventional PEM predicts no remaining response after the primary waveform passed through the slab, whereas the modified PEM results in a 'spread' effect, especially for large values of A()L as presented by Fig.(5-9). Such spread could affect the detectability of the pulses, and consequently, the error code rate of a digital system. 5.4 Discussion and Summary The scattering effects of a random slab on a string of periodic rectangular pulses are formulated and calculated using the results of Chapter IV. The solutions are compared with those obtained using the conventional PEM. The solutions of the modified PEM have been shown to obey the energy conservation law, and all the results are physically explainable. For weak fluctuations and thin slab sizes, the results agree with those obtained by the conventional parabolic equation method. For more severe disturbance and thicker random slab, the modified PEM solutions are still meaningful whereas the solutions by the conventional PEM are doubtful at best, suggestng that the modified parabolic equation method is more suitable for those situations.

65 Returned Pulse: Maximum Magnitude = 3.59E-01 Total Returned Energy = 3.27E-01 Transmitted Pulse: Maximum Magnitude = 6.68E-01 Total Transmitted Energy = 6.63E-01 Reconstructed Incidence: Magnitude = 1.1E+00 Total Incident Energy = 9.9E-01 PARAMETERS A(0) / Wavelength = 0.020, L / Wavelength = 2.5 Carrier Frequency = 10.0, Period = 4.0 Number of Harmonics Calculated: 200 Figure 5-1. Scattering of a Rectangular Pulse (Data set No. 1)

66 Returned Pulse: Maximum Magnitude = 5.11E-01 Total Returned Energy = 5.12E-01 Transmitted Pulse: Maximum Magnitude = 4.77E-01 Total Transmitted Energy = 4.76E-01 Reconstructed Incidence: Magnitude = 1.1E+00 Total Incident Energy = 9.9E-01 PARAMETERS A(0) / Wavelength = 0.020, L / Wavelength = 5.5 Carrier Frequency = 10.0, Period = 4.0 Number of Harmonics Calculated: 200 Figure 5-2. Scattering of a Rectangular Pulse (Data set No.2)

67 Returned Pulse: Maximum Magnitude = 3.40E-01 Total Returned Energy = 5.14E-01 Transmitted Pulse: Maximum Magnitude = 4.47E-01 Total Transmitted Energy = 4.77E-01 Reconstructed Incidence: Magnitude = 9.9E-01 Total Incident Energy = 9.9E-01 PARAMETERS A(0) / Wavelength = 0.010, L / Wavelength = 11.0 Carrier Frequency = 10.0, Period = 4.0 Number of Harmonics Calculated: 200 Figure 5-3. Scattering of a Rectangular Pulse (Data set No.3)

68 Returned Pulse: Maximum Magnitude = 8.72E-01 Total Returned Energy = 7.01E-01 Transmitted Pulse: Maximum Magnitude = 3.12E-01 Total Transmitted Energy = 2.89E-01 Reconstructed Incidence: Magnitude = 1.1E+00 Total Incident Energy = 9.9E-01 PARAMETERS A(O) / Wavelength = 0.100, L / Wavelength = 2.5 Carrier Frequency = 10.0, Period = 4.0 Number of Harmonics Calculated: 200 Figure 5-4. Scattering of a Rectangular Pulse (Data set No.4)

69 Returned Pulse: Maximum Magnitude Total Returned Energy = 1.64E-01 = 9.50E-02 Transmitted Pulse: Maximum Magnitude = 8.29E-01 Total Transmitted Energy = 8.27E-01 Reconstructed Incidence: Magnitude = 9.9E-01 Total Incident Energy = 9.9E-01 PARAMETERS A(0) / Wavelength = 0.002, L / Wavelength = 10.0 Carrier Frequency = 10.0, Period = 4.0 Number of Harmonics Calculated: 200 Figure 5-5. Scattering of a Rectangular Pulse (Data set No.5)

70 Returned Pulse: Maximum Magnitude = 3.39E-01 Total Returned Energy = 4.91E-01 Transmitted Pulse: Maximum Magnitude = 4.77E-01 Total Transmitted Energy = 5.01E-01 Reconstructed Incidence: Magnitude = 9.9E-01 Total Incident Energy = 9.9E-01 PARAMETERS A(0) / Wavelength = 0.010, L / Wavelength - 10.0 Carrier Frequency = 10.0, Period = 4.0 Number of Harmonics Calculated: 200 Figure 5-6. Scattering of a Rectangular Pulse (Data set No.6)

71 By Standard PEM: Maximum Magnitude = 9.97E-01 By Modified PEM: Maximum Magnitude = 6.68E-01 J' -' Reconstructed Incidence: Magnitude = 1.1E+00 PARAMETERS A(0) / Wavelength = 0.020, L / Wavelength = 2.5 Carrier Frequency = 10.0, Period = 4.0 Number of Harmonics Calculated: 200 Figure 5-7. Transmitted Pulse: comparison (Data set No.1)

72 By Standard PEM: Maximum Magnitude = 9.97E-01 By Modified PEM: Maximum Magnitude = 4.77E-01 Reconstructed Incidence: Magnitude = 1.1E+00 PARAMETERS A(0) / Wavelength = 0.020, L / Wavelength = 5.5 Carrier Frequency = 10.0, Period = 4.0 Number of Harmonics Calculated: 200 Figure 5-8. Transmitted Pulse: comparison (Data set No.2)

73 By Standard PEM: Maximum Magnitude = 9.97E-01 By Modified PEM: Maximum Magnitude = 3.12E-01 Reconstructed Incidence: Magnitude = 1.1E+00 PARAMETERS A(0) / Wavelength = 0.100, L / Wavelength = 2.5 Carrier Frequency = 10.0, Period = 4.0 Number of Harmonics Calculated: 200 Figure 5-9. Transmitted Pulse: comparison (Data set No.4)

SUMMARY The two-frequency coherence functions of a continuous random medium are formulated using functional analysis technique. The asymptotic solution for a plane wave incidence is derived and presented in the form of a 'transmitting matrix'. For a single flat slab of random medium in free space, the two-frequency mutual coherence functions of the forward- and backscattered fields are solved by applying the appropriate boundary conditions. The solutions are shown to be physically meaningful and compare favorably with the standard parabolic equation method results. Although the order of approximation is hard to decide, from the comparisons of the computed results, the modified PEM solutions seem to present a more realistic picture of the scattering mechanism, and thus are suitable in studies of severe scattering situations, be they longer propagation distances or larger fluctuations in the medium. Since this formulation views the forward scattering and the backscattering as a mutually coupled process, the results not only do appear to be an improvement over the conventional PEM, but also provide a means to calculate the backscattered power density. 74

APPENDIX 75

76 Listing of the Computer Program Used in Chapter V C PROGRAM PULSE C C March 28, 1985 C by Hsiao-fei Maa C C C C This program is designed to calculate the scattering C effects of a random slab upon a rectangular pulse C predicted by the Modified Parabolic Equation Method. C C The program generates two plots for each data set. C First, the transmitted and backscattered pulses are C plotted along with the incident pulse reconstructed C from its Fourier series expansion, then on the second C plot the transmitted pulses calculated from both C the standard- and the modified parabolic equation C methods are plotted for comparison. C C Required input data are: C C AK = A(O) / Carrier Wavelength C LK = L / Carrier Wavelength C M = Period of the rectangular pulse string C N = Carrier frequency C NHAR = Number of harmonics to be eveluated C TMIN,TMAX define the time interval for which C the results are calculated and plotted. C TOL = Tolerance for all the harmonics C C Cares must be taken that: C C Matrix SI(N1,N1) be dimensioned no C less than 2M by 2M C C TMIN be less than zero in order to show C the front edges of the pulses. C C TMAX be chosen carefully. If it is set too C large, resulted pulses may be too narrow on the C plot to show any details, whereas the back edges of C the pulses may not be plotted completely if TMAX is C not large enough. C C M be large enough that the influence of the C last pulse died out before the arrival of the C next pulse so that the results can be considered the

77 C approximation of transcient responses to a single C pulse incidence. C C Each data set should be look like: C C 81234567 890123456789012345678901234567890 C &DATA AK=.001,LK=5.,NHAR=150,TMAX=2.5,M=5,N=20, C TOL=1.D-3,TMIN=-.3 &END C C Running command: C C $R Object+*PLOTSYS 5=Datafile 9=Plotfile C C Plot command: C C $R *CCQUEUE PAR=Plotfile C C Plots can be previewed on a graphic terminal C by the following command: C C $R NEW:PLOTSEE 0=Plotfile C C The CPU time required on the Amdahl 5860 computer C is approximately 60 to 100 seconds for each data C set, depending on the parameters used. C C IMPLICIT REAL*8(A-H,O-Z,$) COMPLEX*16 GI,GP,GN,CBO,CB1,CB2,CBI,CBP,CBN,CT,CD1, & CD2,CD3,CD4,CD5,CPI(201),CPP(201),CPN(201),CPS(201) REAL*8 L,SI(40,40),PS(201),PP(201),PN(201), & T(201),PIN(201) REAL*4 AK,LK,LDA,TMAX,TMIN,RM,RMN COMMON /GRAPH/PSMAX,PIMAX,PPMAX,PNMAX,EI,EP,EN, & LK,AK,TMIN,TMAX,NT,JPLOT,M,N,NHAR/CORR/CD1, & CD2,CD3,CD4,CD5,D1,D2,ANGLE,V2,LDA NAMELIST/DATA/AK,LK,M,N,NHAR,TMIN,TMAX,TOL C NT=201 JPLOT=0 PI=DARCOS(-1.DO) ANGLE=65536*PI GI=DCMPLX(1.DO,0.DO) C i00 READ(5,DATA,END=200) C A=2.D0*PI*AK L=2.DO*PI*LK LDA=LK/AK RM=1./FLOAT(M) RMN=1./FLOAT(M*N) PM=PI *RM MM1=M-1 M2=2*M

78 M2Ml=M2 -1 DT= (TMAX-TMIN) /DFLOAT (NT-i) DF=RMN*A A2=2.*A DF2=2.*DF C DO 2. I=1,M2MJl SI (I,M)=PM*DSIN(PM*I) DO 1 J=1,M2Ml IF(J.NE.M) SI(I,J)=DSIN(PM*I)*DSIN(PM*J) 2. CONTINUE C C Zero'th harmonic term. C X=.25*L*A BI=.5*PM**2 BP=BI/(l.+X)N BN=BP*X C 2 J=J+l BO=O.DO Bl=O.DO B2=0.DO C DO 3 JJ=1,MM~11 J=J+1 Xl=X*(l.+RMN*J)**2 X2=X*Gl.-RMN*J)**2 5=. 5*SI(JJ, JJ)/FLOAT(J*J) BO=BO+2.* Bl=Bl+S*(l./(l.+Xl)+l./(l.+X2)) B2=B2+S*(Xl/(l.+Xl)+X2/(l.+X2)) 3 CONTINUE BI=BI+BO BP=BP+BI. BN=BN+B2 IF(BO/BI.GT.TOL) GOTO 2 C DO 4 I=1,NT T( I)=TMIN+DT*FLOAT( I-1) CPI (I )=DCMPLX(BI, 0.DO) CPP( I)=DCMPLX(BP, 0.DO) CPN( I)=DCMPLX(BN, 0.DO) CPS(I)=CPI(I) 4 CONTINUE C K.K=O DO 11 K=1,NHAR KK=KK+l IF(KK.EQ.M2+1) KKJl JMK=M2-K 5 IF(JMK.LE.O) JMK=JMK+M2 IF(JMK.LE.O) GOTO 5

79 JPK=K 6 IF(JPK.GT.M2) JPK=JPK-M2 IF(JPK.GT.M2) GOTO 6 C CBP=DCMPLX(O. DOF 0.DO) CBN=CBP CBI=CBP V=K*DF U=A2-V UJl=U U2=U C C Quantities associated with Vr for GAMA. C V2=(.125*V)**2 D=DSQRT( 16. *V2+1 ) Dl=DSQRT(-.5+.5*D) D2=-DSQRT(.5+.5*D) CD1=DCMPLX(l.DO,.375*V) CD2=CD1-V2 CD3=DCMPLX(40.*V2,V*(6.*V2-1L.)) CD4=DCMPLX(l.DO,.125*V)/16. CD5=4.*V2*CD1 C, C Coefficient of the K'th harmonic. C J=-1 JC=-1 IF(KK.EQ.M.OR.KK.EQ.M2) GOTO 7 CALL GAMA(GP,GNrU,V) CBI=SI (KK,M)/FLOAT(K) *GI CBP=GP*CBI CBN=GN*CBI C 7 J=J+l JC=JC+1 IF(JC.EQ.M2) JC=0 CBO=DCMPLX(0.DO, 0.DO) CB1=DCMPLX(0.D0, 0.DO) CB2=DCMPLX (O. DOr 0. DO) C DO 9 JJ=1,MMf11 J=J+1 JC=JC+1 JMK=JMK+l IF(JMK.EQ.M2+1) JMK=l JPK=JPK+l IF(JPK.EQ.M2+1) JPK=l Ul=Ul+DF2 U2=U2-DF2 IF((JMK.EQ.M.OR.JMK.EQ.M2).AND.J.NE.K) & GOTO 8 IF(J.EQ.K) S=SI(JCM)/FLOAT(J) IF(J.N~E.K) S=SI(JCJMK)/FLOAT(J*(J-K))

80 CALL GAMA(GP,GN,UlV) CBO=CBO+S*GI CBi=CB1+S*GP CB2=CB2+S*GN 8 IF(JPK.EQ.M.OR.JPK.EQ.M2) GOTO 9 S=SI(JC,JPK)/FLOAT(J*(J+K)) CALL GAMA(GP,GNU2,V) CBO=CBO+S*GI CBi=CBI+S*GP CB2=CB2+S*GN 9 CONTINUE C CBI=CBI+CBO CBP=CBP+CB1 CBN=CBN+CB2 JMK=JMK+1 IF(JMK.EQ.M2+1) JMK=i JPK=JPK+i IF(JPK.EQ.M2+1) JPK=i Ui=Ul+DF2 U2=U2-DF2 C IF(CDABS(CBO/CBI).GT.TOL.OR.J.LT.K+2) GOTO 7 C DO 10 I=INT CT=CDEXP(DCMPLX(0.D0,(1.-2.*T(I))*PM*K)) CPI(I)=CPI(I)+CBI*CT CPP(I)=CPP(I)+CBP*CT CPN(I)=CPN(I)+CBN*CT CT=CT*CDEXP(LDA*DCMPLX(-8.*V2,V)) CPS(I)=CPS(I)+CBI*CT 10 CONTINUE 11 CONTINUE C CBI=DCMPLX(0.DO,0.DO) CBP=CBI CBN=CBI PIMAX=O.DO PPMAX=0.DO PNMAX=0.DO PSMAX=.DO DO 12 I=iNT PIN(I)=CDABS(CPI(I)+DCONJG(CPI(I))) PP(I)=CDABS(CPP(I)+DCONJG(CPP(I))) PN(I)=CDABS(CPN(I)+DCONJG(CPN(I))) PS(I)=CDABS(CPS(I)+DCONJG(CPS(I))) CBI=CBI+CPI(I) CBP=CBP+CPP(I) CBN=CBN+CPN(I) PIMAX=DMAX(PIMAXPIN(I)) PPMAX=DMAX(PPMAX,PP(I)) PNMAX=DMAX(PNMAX,PN(I)) PSMAX=DMAX(PSMAXPS(I)) 12 CONTINUE

81 C C Calculate total energy contained in each pulse C CBI=CBI-.5*(CPI(1)+CPI(NT)) CBP=CBP-.5*(CPP(1)+CPP(NT)) CBN=CBN-.5*(CPN(1)+CPN(NT)) EI=DT/PI**2*CDABS(CBI+DCONJG(CBI)) EP=DT/PI**2*CDABS(CBP+DCONJG(CBP)) EN=DT/PI**2*CDABS(CBN+DCONJG(CBN)) C C Normalize the pulses before plotting. C DO 13 I=1,NT PIN(I)=PIN(I)/PIMAX PP(I)=PP(I)/PPMAX PN(I)=PN(I)/PNMAX PS(I)=PS(I)/PSMAX 13 CONTINUE C C Calculate the peak values. C PIMAX=PIMAX/PI**2 PPMAX=PPMAX/PI**2 PNMAX=PNMAX/PI**2 PSMAX=PSMAX/PI**2 C C Plot the calculated pulses. C CALL PLOT(T,PIN,PP,PN,PS) JPLOT=JPLOT+1 IF(JPLOT.EQ.4) JPLOT=0 IF(JPLOT.EQ.0) CALL PLTEND GOTO 100 200 CALL PLTEND STOP END C C FUNCTION DMAX(D1,D2) C C Pick out the variable with C larger absolute value. C IMPLICIT REAL*8(A-H,O-Z,$) Cl=DABS(D1) C2=DABS(D2) C=C1 IF(C.LT.C2) C=C2 DMAX=C RETURN END C C SUBROUTINE PLOT(T,PIN,PP,PN,PS)

82 C C Plot the pulses using the University C: of Michigan Plotting System. REAL*8 T(NT),PP(NT),PN(NT),PIN(NT),PS(NT), & PIMAX,PNMAX,PPMAX,PSMAXEIEPEN REAL*4 LK,TMAX,TMIN LOGICAL TL(1O)/'PARAMETERS'! COMMON /GRAPH/PSMAX,PIMAXPPMAXPNMAX,EIEPENr & LK,AK,TMIN,TMAXNTJPLOTM,NNHAR C CALL PALPHA('SANSERIF.1 ',O) XORG=2. YORG=2.5+FLOAT(JPLOT)*7.8 DX=(TMAX-TMIN)/5. DY l. C XT=4.5-.5*PSMLEN(TLlO,.15) YT=YORG-1.5 CALL PNUMBR(2.7,YT-.5,.lNHAR,O., & '"Number of Harmonics Calculated ",JI3*1,O) CALL PNUMBR(2.7,YT-.35,r.J,FLOAT(N),O., & "'Carrier Frequency = ",JF3.1*",O) CALL PNU`MBR(-O.,-O.r.1,FLOAT(M),O., & '", Period = ",JF2.1*,O) CALL PNUMBR(2.7,YT-.2,.l,AK,O., & "'A(O) / Wavelength = ",JF2.3*',O) CALL PNUMBR(-O.,-O..1 LKO., & '", L / Wavelength = ",JF2.1*",O) CALL PSYM(XTYT,.15,TL,O.,10,0) C CALL PNUMBR(2.7,YORG-.45,.lEIO., & "'Total Incident Energy =",WE1.1*I,O) CALL PNUMBR(2.7,YORG-.3,.1,PIMAXO., & "'Reconstructed Incidence: Magnitude =",WE1.l*',0) CALL PLTOFS(TMINDXO.,DYXORGYORG) CALL PLINE(TPINNT,2,OO,1) C CALL PNUMBR(2.7,YORG+1.55,.1,EPO., & "'Total Transmitted Energy =",WEl.2*,0) CALL PNUMBR(2.7,YORG+1.7,.1,PPMAXO., & "'Transmitted Pulse: Maximum Magnitude =",WEl.2*,0) CALL PLTOFS(TMINDXO.,DY,XORGYORG+2.) CALL PLINE(TPPNT,2,OO,1) C CALL PNUMBR(2.7,YORG+3.55,.l,ENO., & "'Total Returned Energy =",WE1.2*1,0) CALL PNUMBR(2.7,YORG+3.7,.lPNMAXO., & "'Returned Pulse: Maximum Magnitude =",WE1.2*',O) CALL PLTOFS(TMINDXO.DY,XORGYORG+4.) CALL PLINE(TPN,NT,2,OO,1) C JPLOT=JPLOT+l

83 C YORG=2. 5+FLOAT(JPLOT) *7.8 YT=YORG-1.5 CALL PNUMBR(2.7,YT-.5,.l,NARO., & "'Number of Harmonics Calculated - ",JI3*',O) C C C C C C C C C C C C C C C CALL PNUMBR(2.7,YT-.35,.lFLOAT(N),O., & '"Carrier Frequency = ",JF3.1*',O) CALL PNUMBR(-O.,-O.,.lFLOAT(M),O., & '', Period = ",JF2.1*',O) CALL PNUMBR(2.7,YT-.2,.lAKO., & "'A(O) / Wavelength = ",JF2.3*",O) CALL PNUMBR(-O.,-O.,.lLK,O., & '", L / Wavelength = ",JF2.*",0) CALL PSYM(XT,YT,.15,TLO,10, O) CALL PNUMBR(2.7,YORG-.3,.1,PIMAXO., & "'Reconstructed Incidence: Magnitude CALL PLTOFS(TMIN,DXO.,DYXORGYORG) CALL PLINE(TPINNT,2, O,,1) CALL PNUMBR(2.7,YORG+1.7,.lPPMAXO., & "'By Modified PEM: Maximum Magnitude CALL PLTOFS(TMINDXO.DYXORGYORG+2.) CALL PLINE(T,PP,NT,2,OrO,l) CALL PNUMBR(2.7,YORG+3.7r.1,PSMAXO., & "'By Standard PEM: Maximum Magnitude CALL PLTOFS(TMINDXO.,DYXORGYORG+4.) CALL PLINE(TPSNT,2,rO,O,) =tv WE1.1*1,0) = I WEI..2 * J' ', 0 =11 WEI..2*'rO) RETURN END SUBROUTINE GAMA(GPGNUK,VK) Calculate the two-frequency mutual coherence functions GP and GN for the forward- and backscattered fields. UK and VK are inputs of this subroutine. VK must be non-zero. IMPLICIT REAL*8(A-HOZr$) COMPLEX*16 GP,GNLA1,LA2,E1,E2,E3,E4,CD1, & CD2,CD3,CD4,CD5,Cl,C2,S1,S2,UPUN,DOWN REAL*8 IT REAL*4 LDA COMMON /CORR/CDlCD2,CD3,CD4,CD5, & D1,D2,ANGLEV2,LDA U=DABS(UK) V=DABS(VK) USQ=U**2

84 D3=DSQRT( (USQ/V+3.*V) **2/64.+l) C ZE=V*DSQRT(-.5")+.5*D3) IT=.-V*DSQRT(.5+.5*D3) LA2=DCMPLX(ZErIT) ZE=ZE*LDA IT=IT*LDA 1 IF(IT.LT.-ANGLE) IT=IT+ANGLE IF(IT.LT.-ANGLE) GOTO 1 C IF(U.NE.O.DO) GOTO 2 FACT=DEXP (-ZE) E3=CDEXP(DCMPLX(O.DO, IT)) E4=CDEXP(DCMPLX(-2 *ZE, -IT)) C2=.5* (E3+E4) S2=.5*(E3-E4) UN=V2*(FACT-C2)+LA2/16.*S2 DOWN=V2 *FACT+CD2 *C2 - & LA2*DCMlPLX(.1875DOjKl./V)*S2 GOTO 4 C 2 CONTINUE AL=U*DJl BE=U*D2 LA1=DCMPLX(AL, BE) AL=AL *LDA BE=BE*LDA 3 IF(BE.LT.-ANGLE) BE=BE+ANGLE IF(BE.LT.-ANGLE) GOTO 3 C ARG=DMAX(AL, ZE) FACT=DEXP (-ARG) El=CDEXP(DCMPLX(AL-ARG, BE)) E2=CDEXP(DCMPLX(-AL-ARG, -BE)) E3=CDEXP(DCMPLX(ZE-ARG, IT)) E4=CDEXP(DCMPLX(-ZE-ARG, -IT)) C1=.*EE2 S1=.*EE2 C2=.5* (E3+E4) S2=.5* (E3-E4) UN=V2*(Cl-C2)+DCMPLX(O.DOUSQ*V/64.) & /LA1*Sl+(USQ*CD4-CD5)/LA2*S2 DOWN=UN+CD1*C2+CD3/LA2 *S2 C 4 UP=CD1*FACT GP=UP/DOWN GN=UN/DOWN IF(VK.LT.O.DO) GP=DCONJG(GP) IF(VK.LT.O.DO) GN=DCONJG(GN) C RETURN END

REFERENCES 1. P. R. Jedrzejewski, Backscatter from a Random Medium, Ph.D dissertation, The University of Michigan, 1982. 2. A. Ishimaru, "Pulse Propagation, Scattering and Diffusion in Scatterers and Turbulence," Radio Science, Volumn 14, No.2, 1979. 3. V. I. Tatarskii, The Effects of the Turbulent Atmosphere on Wave Propagation, National Technical Information Service, U.S. Department of Commerce, 1971. 4. Roman. Introduction to Quantum Field Theory, New York, Wiley, 1969. 5. A. Ishimaru, Wave Propagation and Scattering In Random Media, Volumn II, Academic Press, 1978. 6. A. Friedman, Partial Differential Equations of Parabolic Type. Prentice-Hall, Inc., 1964. 7. B. Noble and J. W. Daniel, Applied Linear Algebra, Prentice-Hall, Inc., 1977. 8. C. T. Tai, Dyadic Green's Functions in Electromagnetic Theory, Intext Educational Publishers, 1971. 9. J. D. Jackson, Classical Electrodynamics, 2nd ed., John Wiley & Sons, Inc., 1975. 10. A. Ishimaru, "Theory and Application of Wave Propagation and Scattering in Random Media," Proceedings of the IEEE, Vol. 65, No. 7, July 1977. 11. I. Sreenivasiah, A. Ishimaru, and S. T. Hong, "Two-frequency mutual coherence function and pulse propagation in a random medium: An analytic solution to the plane wave case," Radio Science, Vol. 11, No. 10, October 1976. 12. A. S. Gurvich and V. I. Tatarskii, "Coherence and intensity fluctuations of light in the turbulent atmosphere," Radio Science, Vol. 10, No. 1, January 1975. 13. A. Ishimaru, "Pulse propagation, scattering, and diffusion in scatterers and turbulence," Radio Science, Vol. 14, No. 2, March-April 1979. 85

86 14. S. T. Hong, I. Sreenivasiah, and A. Ishimaru, "Plane wave pulse propagation through random media," IEEE Trans. on Antennas and Propagation, Vol.AP-25, No.6, November 1977. 15. A. Heins, "Methods of Applied Mathematics," MATH. 556 & 557 classnotes, The University of Michigan, 1982. 16. R. H. Kraichnan, "The Closure Problem of Turbulence Theory," Proc. of Symposia in Applied Mathematics, Vol. XIII, American Mathematical Society, Providence, RI, 1962. 17. U. Frisch, "Wave Propagation in Random Media," Probability Methods in Applied Mathematics, Vol. I, Academic Press, 1968.