RL-586 EXCITATION AND PROPAGATION OF WAVES BETWEEN TWO PLANAR INTERFACES by Yu-Ping Liu A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1973 Doctoral Committee: Professor Chen-To Tai, Co-Chairman Professor Chiao-Min Chu, Co-Chairman Professor Vi-Cheng Liu Professor Charles B. Sharpe Dr. Dipak L. Sengupta RL-586 = RL-586

ABSTRACT EXCITATION AND PROPAGATION OF WAVES BETWEEN TWO PLANAR INTERFACES by Yu-Ping Liu Co-Chairmen: Chen-To Tai, Chiao-Min Chu The theoretical study of wave propagation in the free space region between two parallel planar interfaces is presented. The upper part of the region is bounded by a perfectly conducting plate while the lower part is bounded by a lossy dielectric earth. The dyadic Green's function for this geometry is developed as an Ohm-Rayleigh type of expansion in terms of the Hansen vector wave functions. This dyadic Green's function is used to obtain the field expressions for Hertzian dipoles of electric and magnetic type with different orientations. Methods of evaluating the integral expression for the fields are studied in great detail. The pole contributions from the integral representation of the field yield the modal picture of the surface wave. The excitation factor and the mode function for the wave expressions are derived. The dependence of the real and imaginary parts of the poles on the property of the lossy dielectric medium is discussed. By means of the deformation of integration contour and the saddle point method, the far field expressions for the waves excited by different elementary sources are derived. The proper choice of branch cut in performing the contour integral is discussed. Some numerical results are presented in graphical form. The results obtained by using the Leontovich impedance boundary condition on the lossy interface are compared with the exact results obtained by using the continuity of tangential fields across the interface. The symmetry properties of the dyadic Green's function in a two layer problem are investigated.

A CKNOWLEDGMENTS I wish to express my appreciation to all the members of my doctoral committee for their comments. My deep indebtedness is given to Professor Chen-To Tai who suggested the area of research and provided continuous inspiration, guidance and encouragement. Also special thanks are due Professor Chiao-Min Chu for his advice and valuable criticism. Thanks are also due Dr. Harold E. Foster, Dr. Adel Mohsen and Mr. Gerard Desjardins for their helpful discussion. The careful typing of Mrs. Karen L. Christiansen is appreciated. Finally, I would like to take this chance to express my appreciation to my wife, Chi-ling, for her understanding and patient consideration during the course of study. ii

TABLE OF CONTENTS ACKNOWLEDGMENTS ii LIST OF TABLES v LIST OF ILLUSTRATIONS vi LIST OF APPENDICES vii I INTRODUCTION 1 II DYADIC GREEN'S FUNCTIONS 4 2.1 General Remarks 4 2.2 Free Space Green's Function 4 2. 3 Composition of Green's Function for Electric Dipoles 6 2.4 Composition of Green's Function for Magnetic Dipoles 9 2.5 Symmetry Property of Dyadic Green's Function 12 2.6 Impedance Boundary Condition 17 2. 7 Approximate Dyadic Green's Function Satisfying the 17 Impedance Boundary Condition 2.8 Relationship Between Exact and Approximate Green's 20 Function mi FIELDS DUE TO HERTZIAN DIPOLES 21 3.1 Introduction 21 3.2 z-directed Dipole 21 3.2.1 Integral Representation 21 3.2.2 Asymptotic Solution 23 3.2. 3 Comparison with Mode Theory Results 27 3.2.4 Leontovich Approximate Solution 29 3.3 x-directed Dipole 30 3.3.1 Exact Solution 30 3.3. 2 Leontovich Approximate Solution 32 3.4 Horizontal Current Loop 32 3.5 Vertical Current Loop 33 IV CALCULATION AND DISCUSSION OF RESULTS 34 4.1 General Remarks, - 34 4.2 Percentage Error by Using Impedance Boundary Condition 34 4.3 Far-.Zone Field 40 4.4 Distribution of Poles in Complex 0 Domain 46 4.5 Excitation Factor for Various Sources 49 iii

V CONCLUSIONS AND RECOMMENDATIONS 54 5.1 Conclusions 54 5.2 Areas for Future Study 54 REFERENCES 56 APPENDICES 58 iv

LIST OF TABLES Table Page 4-1 Refractive Index for Different Frequency 35 42 Percentage Error of Reflection Coefficients 36 4-3 Percentage Error for Far Field 46 4-4 Table of Excitation Factor for Wave 48 Components for Different Sources 4-5 Location of Poles for Good Earth 50 4-6 Location of Poles for Poor Earth 50 v

LIST OF APPENDICES APPENDIX A: General Discussion on Saddle Point Integration 58 APPENDIX B: Newton-Raphson's Method 64 APPENDIX C: Discussion of Pole Distribution 66 APPENDIX D: Computer Program 70 vii

LIST OF ILLUSTRATIONS Figure Page 1-1 Geometric Configuration 1 3-1 Path of Integration in Complex X Plane 25 3-2 Path of Integration in Complex f3 Plane 25 4-1 Magnitude of Reflection Coefficient at 6 MHz 37 4-2 Magnitude of Reflection Coefficient at 600 KHz 38 4-3 Magnitude of Reflection Coefficient at 60 KHz 39 4-4 Far-Zone Electric Field of Unit Dipole Moment 41 at 6 MHz 4-5 Far-Zone Electric Field of Unit Dipole Moment 42 at 3600 KHz 4-6 Far-Zone Electric Field of Unit Dipole Moment 43 at 60 KHz 4-7 Far-Zone Electric Field of Unit Dipole Moment 44 at 6( MHz (z 5m) o 4-8 Far-Zone Electric Field of Unit Dipole Moment 45 at (3 MHz (z = 25m) o 4-9 Pole Distribution for Good Earth 51 4-10 Pole Distribution for Poor Earth 52 2 2 A-1 Branch Cut Along Re k - X 0 58 A-2 Branch Cut Along Im k- X =0 59 A-3 Branch Cut Parallel to Imaginary Axis 60 A-4 Integration Path in Complex 13 Plane 61 C-1 Pole Locus at Cut-Off Frequency 67 C-2 Pole Locus Above Cut-Off Frequency 69 vi

ABSTRACT EXCITATION AND PROPAGATION OF WAVES BETWEEN TWO PLANAR INTERFACES by Yu-Ping Liu Co-Chairmen: (hen-To Tai, Chiao-Min Chu The theoretical study of wave propagation in the free space region between two parallel planar interfaces is presented. The upper part of the region. is bounded by a perfectly conducting plate while the lower part is bounded by a lossy dielectric earth. The dyadic Green's function for this geometry is developed as an Ohm-Rayleigh type of expansion in terms of the Hansen vector wave functions. This dyadic Green's function is used to obtain the field expressions for Hertzian dipoles of electric and magnetic type with different orientations. Methods of evaluating the integral expression for the fields are studied in great detail. The pole contributions from the integral representation of the field yield the modal picture of the surface wave. The excitation factor and the mode function for the wave expressions are derived. The dependence of the real and imaginary parts of the poles on the property of the lossy dielectric medium is discussed. By means of the deformation of integration contour and the saddle point method, the far field expressions for the waves excited by different elementary sources are derived. The proper choice of branch cut in performing the contour integral is discussed. Some numerical results are presented in graphical form. The results obtained by using the Leontovich impedance boundary condition on the lossy interface are compared with the exact results obtained by using the continuity of tangential fields across the interface. The symmetry properties of the dyadic Green's function in a two layer problem are investigated.

CHAPTER I INTRODUCTION The purpose of this thesis is to theoretically study electromagnetic wave propagation in the vicinity of two parallel plates of different material. In particular, one plate considered is the lossy ground plane; the other is the perfectly conducting screen above it. Two major steps are involved in the calculation of the far-zone field of a transmitting antenna. First, the integral expressions for the electromagnetic fields need to be derived by using the dyadic Greents function technique. Secondly, the integration involved is carried out by the saddle point integration method, __________ The geometry of the problem considered is illustrated in Fig. 1-1. The FIG. 1-1: GEOMETRIC CONFIGURATION dipole is located at a height z above the earth. The earth is considered to be a homogeneous medium with finite conductivity as in Sommerfeld's problem [16]. The distance between the perfectly conducting top plate and the earth is L which 1

2 is larger than z. The medium in between the plates is free space. 0 The problem of finding the electromagnetic wave propagation between the ionosphere and the earth has been discussed by many authors. Watson II] considered it as a propagating waveguide of electromagnetic wave. He employed a waveguide approach (between an idealized homogeneous earth and a concentric reflecting layer) in which he devised the celebrated transformation in rapid convergence. The same model with certain refinements has been discussed by many others; in particular, Bremmer [2], Budden [3,4], Wait [5] compiled many useful results in book form, utilizing either mode theory or Hertzian potential technique. In general, they considered the earth as highly conductive. The problem discussed here, unlike the ionosphere problem of low frequency and large distance, involves high frequencies and a very short separation distance. However, both have similar configurations. In this dissertation the method used to calculate the fields is that of the dyadic Green's function. The derivation follows the scattering superposition scheme suggested by Tai 76, 7]. The use of the dyadic Green's function technique to solve the wave propagation problem provides a general treatment of both excitation and wave mode at the same time. The exact integral represention of the field is obtained by using composite dyadic Green's functions. The integration is treated by the saddle point method. Thus, the solution at a large distance from the source is investigated. In contrast, the mode theory approach requires consideration of the excitation factor separately or a special source function must be chosen to obtain the right kind of potential. In Chapter II the dyadic Green's functions for mixed boundary conditions are found. Two kinds of Green's functions are needed for two different sources. When the material is highly conductive, an approximate impedance boundary condition is imposed on the lossy earth. The usage of the Leontovich boundary condition replacing the exact tangential continuous condition is termed as the approximate solution. A symmetry property is derived for the two layer prob

3 lem in our research. In Chapter III the fields due to Hertzian dipoles are derived. The dipole sources considered are of electric and magnetic type with both horizontal and vertical orientation. The far-field is evaluated by means of the saddle point method. The results obtained by the method of Green's functions are compared with the miode theory of the waveguide. Chapter IV contains the results of various computer aided numerical calculations. The ground parameters were chosen to be representative of those actually encountered in practice. Several frequencies were chosen to study the nature of the ground. Some overlooked comments on wave propagation in the literature are discussed. The error in the results when using the impedance boundary condition is found. The advantages of using the dyadic Green's function technique are shown by eliminating the determination work on the choice of Hertzian potentials and the excitation factors. A general discussion of mode distribution is included. In Chapter V the conclusions as well as suggested areas where future efforts may be productive are discussed.

CHAPTER II DYADIC GREEN'S FUNCTIONS 2.1 General Remarks The dyadic Green's function technique for treating electromagnetic boundary-value problems was first formulated by Schwinger [8] in 1943. Levine and Schwinger [9] used this method to solve the problem of diffraction by an aperture. In their book, Morse and Feshbach [10] introduced the freespace dyadic Green's functions by way of the vector Helmholtz equation. In 1954, Tai [LI] collected his work on these functions into a report which was developed further and published in book form in 1972 [6]. The general procedure for finding the Green's function related to a particular problem is discussed in great detail in Tai's book. In this thesis, we apply the method of Ohm-Rayleigh to a configuration as shown in Fig. 1-1. We consider only harmonically varying fields. The time factor e is applied in this work. 2.2 Free Space Green's Functions In dealing with the integral representation of fields in free space, the cylindrical vector wave function used has continuous eigenvalues in the X-domain as well as in the h-domain, where X is related to the radial eigenfunction and h is related to the longitudinal eigenfunction. The free space Green's function for the same type of problems is developed by Tai [6]. The method used is the Ohm-Rayleigh or eigenfunction expansion technique. In cylindrical coordinates, the free-space dyadic Green's function is given by [6], [7]: +i N N(h)X Mt (-h) + o 46r X 3! X n d/ oaX h M (-h) (h) + n=O on nX + N (h) Ne (-h) z > z' e _ 0 nX 0 nX 6( - Th') + N (-h)YN' (h) It k2 z <zl enX enX 0 0 (2.1) 4

5 7 AA! where I= zz 2i k:the free space wave number = +h X: eigenvalue pertaining to the radial function h: eigenvalue pertaining to the longitudinal function r1 nw O 6 = O Kronecker delta function 0 0 n 0O and the vector wave functions are defined as: (+ cos + ihZA M (+h)= Vx J (Xr) cos no e nX sn n sin A n cos + ihz nB f n e ---s-os r cos ar sin N (+h) = VxM (+ h) e k e onX onX _iE aJn(Xr) cos i (+h) nc ne (.J(r):n~ X J(Xn). n. e- (2.2) - ar sin r n cos n sin The subscript e corresponds to even angular functions, and (+h) corresponds o odd + ihz to e- and J (X(r) is the nth order Bessel function. n -_____ _ -- ---- A simplified notation for the' sum has been adopted: M (h) M' (-h)= M- (h) M'- h) + M-(h) M - (-h) e e en eno enn on on A, onA N (h) N (-h) ' N-(h) N — (-h) + (-h) N' h) eon onenX e The M' and N' functions are defined in the primed coordinate system, the source coordinates R'. For these vector wave functions, the sign of 2 2 h =+ 2- is chosen so that the imaginary part is always: positive throughout the calculations; this insures that the radiation condition is always satisfied. The free space Green's function is the basic building block in the method

6 of scattering superposition [6], [7]. 2. 3 Composition of Green's Function for Electric Dipoles The plane earth is assumed to be a homogeneous lossy dielectric with a permeability pu equal to that of free space, a relative permittivity eE, and a finite conductivity C E. Associated with the lossy earth is the following propagation constant. - 2 2 2 E E k = k -+ i i — where k 2= p e E o e W o o o The complex index of refraction is defined as 2 E I CE n- +i0 0 The new Green's function can be found in the following way. For convenience, the coordinates are picked so the perfect conductor is at z = 0 while the lossy earth is located at z = L. Applying the technique of scattering superposition, the dyadic Green's function of the first kind with perfect conductor at z = 0 is: G (R/R/1) + ) (=/R') z > 0 (2.3) where G is as previously defined and o G = a40 M e(h) t (h) +bN (h) (h) is 4 I h e h e e, n=0 on 0n nX ~nX J Here a and b are two unknown constants which are specified by the boundary conditions. The posterior elements are guided by the expression of G which has to be the same as the part of G when z < z' and the anterior ele-l 0 0 ments are dictated by the nature of the reflected outgoing wave. Notice that the singular term at source point has been dropped since the singular nature of G1 is taken care of by G. The boundary condition at z x 0 requires that:

7 z x G. (2.4) Substituting Eqs. (2. 1), (2. 2) and (2. 3) into (2.4) yields: a =- b= 1 Thus the complete expression for the dyadic Green's function of the first kind is:. D 2 -6 (h) M(-h)-M(h)M' (h)+N(h) N'(-h) + h 1 4 n= ( (- h) (h)M' (h) N(- (h) M(h) + + N(h) N' (h) z > z' + N(h)N' (h) O < z < z where the abbreviated form of M (h) and N (h) represented by M'(h) and e e onX n; ------ On n N(h) have been used to simplify the notation. Now the dyadic Green's function of the third kind needed for the problem can be constructed from that of the first kind, using the method of scattering superposition. In order to satisfy the radiation condition and the boundary conditions at the interfaces, the Green's function is assumed to have the following form - - r( ) - - - 2 -6 (h) L '(-h h) 2+N(M(h) N'-h(h +N(+N'(h)+ n3 O M(-h) - M(h) M'(h)+ N(-h)+N(h) N'(h) + +R eM(-h) - M(h) IM'(-h) - M'(h) +R -h)+ ( (-h) + N(h) L>z>z e _m J In - +Re [-h)M(h [M'(-h ) - M' i +Rm N(-h) + N(h) N'(-h)+ N '(h) J z'>z>O (2.5) co oo 2 - n0 - - + Tm N(hE) [N'(-h)+ N(h z > L In E [

8 where hE + k -X E- E The first superscript (2) implies that R is in medium 2 and the second superscript 1 means that R' is in medium 1. R R, T and T are four une In e m known coefficients to be determined. The boundary conditions at the dielectric surface which are used to determine the unknown coefficients are: A ((11) A (21) x G (R/R) x 3 (R/R) 3 3 zxL, xVx (R/R) = z xVx (R/R z = L We have assumed /0= 1,= 2. Substituting (2. 5) into the above equations we obtain: - i h E ihL+ -ihL ihL E e +R e e = T e e e k:-e + e T -kE ihihL he.+R -he he = T h e e (2. 6)e r -ihL ke +R e +ke T kEe These linear equations yield ihL (h-hE) e e 2 h cos(hL-) -i hE sin (hL) -ihEL T he (2. 7) e hcos(hL)-ihE sin (hL) Ei

9 (h kE2 hEk2) eihL (hk -h -ke( R E 2- - (2.8) 2 k h-cos(hL)- ikh sin(hk -ihEL hkkEe T k -. (2.9) k hE cos (hL) - ikEh sin (hL) Knowing the complete expression of G3, we can calculate the electric field due to a known current source by E(R)= iwm0 1G (R/R)- J-(R')dV' L> z> 0, E(R) = iW o f G( )(R/R) -J(R')dV' z >L 2.4 Composition of Green's Function for Magnetic Dipoles In principle, once the appropriate dyadic Green's function is known, the electromagnetic field due to any arbitrary current distribution can be found by evaluating the definite integrals involving the current distribution. In practice, for currents in the form of small loops it is more convenient to use the alternative choice of boundary condition. By introducing an equivalent magnetization vector m in place of the current J, the vector wave equation for H becomes VxVxH-k H= kI. (2.10) The magnetization vector m is defined as m= IA, where I denotes an equivalent loop current and A the vectorial area of the loop. To solve Eq. (2. 10) directly, it is convenient to construct a G function which satisfies the boundary condition n x Vx G2= 0. The procedure is exactly the same as before. The dyadic Green's function of the second kind

10 with a perfect conductor at z = 0 is 2(R/R) = G (R/R') + (R/R') 0 2S z >0 (2.11) where 2 i5 4 G2S- 4 2S' 47J dX x =0 n=0 2 -6 r 6 0 c M (h) M'(h)+dN(h)N'(h) h % — The boundary condition at z = 0 requires that z xVx G2= 0 (2.12) Substituting Eqs. (2. 1), (2.2) and (2. 11) into (2. 12) yields the constants c and d: c=1, d=-1. Thus, the G function is given by 2 4 00X 2-6 M(h) M'(-h) + M(h) M'(h) + N(h) N'(-h) - N(h) N'(h) hM(-h)M'(h)+M(h)M'(h) + (-h) N'(h) - N(h)N(h) n=0 z > Z' J z' >z>0. Now, with the dielectric earth on top at z = L, the Green's function of the fourth kind must have the form 4G( 11 )(R-^r ) i J L- N'-, n=O 0 M(-h) + M(h) M'(h) + [N(-h) -N(h) N'(h) +S (-h) +M(h '(-h)+ M'(h) + S L(-h) - N(h (-h) -N'(h) L>z>z' +S.M(-h)+M(h)[ + M'( (h) + S IN-h-)+'(h) -h)N(h) '(-h) -N'(h) z >z (2.13)

11 G 2-6 G (2/R,)= i d jm(hE) M'(-h)+M'(h +PN(h) N'(-h)-N(h) G4 ze E ( n=0 z > L where E= 2 E X The boundary conditions at z = L require: A =(11) zx G 4 = zx z=L 4 z=L 1 A (11) 1 A =(21) -z xVxG - z xVx G E - 4 4 z=L Substituting (2.13) into these equations we obtain: ihL e + S m -ihL e + e =P e -m heihL+S [ h ihL -- e +S En m 1 h -ihL h ihL - e — e = k k J hE ihEL P e ekE r h - i -ihL h ihL] e +-e =P * 1 - m hE ihEL E2 C2 Hence: k- ih -ihL k ih L kE ihEL -^ e +S -e - e - P e E1 e 1 E2 e 2 ihL (hE n h) e S- E 2 i n h sin (hL) - hE cos (hL 2 -ihEL n he h cos(hL)-in sin (hL) - E _o i hL (2.14) (2.15)

12 ihL (h-hE) eh e 2Lhcos(hL)-ihEsin(hL)4 (2.16) -ihEL nhe P - nh e (2.17) e [hcos(hL)-i hEsin(hL (2.17) 2 j where n = 2/1. Thus, knowing the expression for G, an application of Green's identity 4' gives an integral form of the magnetic field due a constant current loop distri- / bution (R)= k2j JG(l1).m(R')dV' L>z>0 2 2 (21)_ H(R)= n~k G(RI')dV' z>L -4 2,5 Symmetry Property of Dyadic Green's Function Analogous to the scalar Green's functions, the dyadic Green's functions have certain symmetry properties which are described in detail by Tai [6]. In his book, only problems involving two media have been discussed. It is not difficult to extend the multi-media problems. In this section, we extend one of Tai's relationships to a two layer problem. =(11) Let us consider the problem in our research. The functions G 3 G3, 2 and G satisfy the following equations: _ 4 4__ _ ___ _/11\ 2 - 1 1 -VxVx G )(R/R)-k2G( )(R/ )= 6 (R-R ) 3 a a a =(21) 2:5/21) xVxVG3 ()(R/R )-kE(3 (R/R )= o a )k 3 a =x x(11 Vx G4 )-k G4 R, (R6 Rh)

13 zx-(2 1) 2==(21) — VxVx 4 (R/) -kEG4 (R/) = 0 where R and Rb are assumed to be in region 1. At the interface z = L, a _... -------- these functions satisfy the boundary conditions that n x G = n x G 3 3 -1 A =(11) 1 A -(21) -nxVxG 3 -nxVx G3 (2.18) A =(11) A -(21) nx G =x G 4 4 1 A =(11) 1 A =(21) -- nxVx G nxVx G E1 4 E2 4 The vector Green's first identity has the form (fL VxP). (VxQ)-Q. VxVx PdV= xVxP).dS V S We form two vector functions defined by =(11) P1 G4b Q1 VxG (R/R ).a An application of the Green's first identity in region 1 yields: - -=(11) bV(VxG(G11 l)ja-a-Vx R / R ).b [ + [ x 4 d VxG R/R) R/) - xG3 /R) * rl )(R/b]} dV =

14 V xG( 11(R/R ). a x Vx G 1 b* dS 3 a 4 1 where - as before means the transpose of the dyad. The volume integral in the above equation can be changed into a surface integral by means of Gauss' theorem since it is equal to k GV- 4 )(R/Rb)-b x G3 (R/Ra)-T As a result of the radiation condition and the boundary condition at the perfectly conducting plane, the resultant surface integral exists only on the interface SL. Thus, we have =(11) - - =(11) b- VxG 4 (R /Rb) * a -a VxG3 (R/Ra) 3 a4 a L + Vx 3 (R/R ). x Vx ). (./) dS (2.19) Similarly, by letting =(21) P2G4 (R/ - =-(21)Q2 VxG3 (R)/R ) a - 3 a and applying Green's first identity to region 2 covering the entire domain z >L, we obtain

15 Sif.i;'j,. [k1)G/4.*^]21) Vx G(dS-) i SEG -(2) =(2b1) ](2 1) i3 * [k G()(R/Rb) * b]}d0V = VxG32 (R/Ra) x Vx G421)(R/%) dS/ 3 a-4 2 As a result of the application of Gauss' theorem to the volume integral and the radiation condition at infinity, the only surface integral left is on the interface S. Thus: 2 =(2 1)1 g( 2 ) - SL VxG (R/R ).xVxG Rb/.'b) dS. (2.20) 3 a 4 SL 2 22 In view of Eqs. (2. 19) and (2. 20) and noticing that k- n k we have b- Vx1 /)-'.Vx RR.b b. x(4 (R /R). a-a. x(3 aR/Ra b 2 k1 ) -k J (R/Rb) b]x G3 (R/R ).' - dS 4 3 a + kf G fGb L' I> j dS L + VXG aR/R1 a x G] )x /b.b] dS SL

16 2 [VXG21 (R/R )* xVxG )(R/R) dS SL When the boundary conditions (2.18) are applied, the surface integral over the interface SL vanishes. Thus 4 a3 a -. Vx =. Vx(11)(R/R) * Putting VxG3 (R/R )= AB, gives aV. x311(Rb/Ra) b =a (B A) *= b. B) (Ab) = (-fA.) (B -):= V. (AB) -)a =F.VxG(3 (Rb/R )a Since a and b are arbitrary constant vectors, =( 11) =(11) VxG4 (R/Rb) = VxG3 (Rb/R Replacing R by R' and R by R, we obtain 'xG4 (R/R)=Vx3 (R/R'). (2.21) 4 - Here the '' covers the whole dyadic function according to Tai's notation. In examining the multi-reflection coefficients R, R, S and S, it is m____ me e found that the following relationships exist: S =R S = R m m e e Thus, the G3 and G4 functions derived for our problem have the symmetry 3J 4X

17 property as is clearly seen in Eq. (2.21). This procedure could easily be extended to other multi-media problems. 2.6 Impedance Boundary Condition During the early 1940's, a large number of Russian papers were published dealing with various aspects of wave propagation over the earth. An attempt was made to take into account the actual properties of ground materials by specifying an impedance boundary condition at the surface. These conditions are attributed to Leontovich (for example, see Fock [12]) and were described to Leontovich [13]. Application of these boundary conditions to propagation problems has been discussed by Feynberg [14]. A straightforward and comprehensive study was done by Hiatt et al [15]; they obtained the following results when I|n >>1: aE z ik -- - ""- E n z z=0 aH =-iknH az z z=O z=0 where the interface is located at z 5 0 and the positive z direction is upward toward the free space. 2. 7 Approximate Dyadic Green's Function Satisfying the Impedance Boundary Condition / Analytically, when the impedance boundary condition is used, the formulation becomes different from the one presented in previous sections. The impedance boundary conditions at the interface z = L imply that there is no wave propagating beyond the interface. The method of obtaining the Green's function remains the same. Using the dyadic Green's function of the first kind as the basis and applying scattering superposition technique, we obtain Gf 2= i M(h) M(-h) -V'(h) + N(h) N'(-h)+ N'(h) + 31 43 r '' hh) + M(-h)-M(h) jM'(h)+ N(-h)+K(h) _'(h) + n= ' j~

18 +ReI M(-h) M(h) I"(-h) -M(h +R N(-h)+N(h (-h)+N (h) L>z> z (-)- ()RI +R eI f(-h) - VWh) (-h)a +R T I ( -h)+'N(h) PN'(-h)+1'N(h) z'>z>0 6J. - I~ll -~ %.- -. % (2 where R and R are unknown coefficients and the subscript I is used to eI mI denote the modified coefficients under the impedance boundary condition. The Leontovich boundary conditions mentioned before correspond, at z = L, to aE -z + i E =0 az n z.22) I --- aH z - +iknH =0 az z Under the assumption that Po= =l /u2 the impedance boundary conditions in dyadic form are: A aG31 + ik.(-a +- n )= z n 31 zL z = L A a ' (- VxG +iknVxG3I) z=L (2.23) Substituting Eq. (2. 2) and (2. 22) into (2. 23) yields: ihL -ihL ihL) ik r ihL+ -ihL ih Li -ihe +R (ihe -ihe )+ e R (e e =0 ml n ml J ihL -ihL ihL ihL -ihL ihLl -ihe +Re (ihe +ihe )+ikn e + (e -e = \Hence, - ihL ) R _ e (k-nh) Rml 2 nh sin (hL) - k cos (hL (2.24)

19 ihL e (nk-h) el 2 k n sin (hL)- h cos (hL)) (225) The same procedure carries over for G4i. It yields: i 2-6o M(h) M '(-h) +M'(h) -+N(h) N'(-h)-N' (h) + 4I 47r~ h M(-h)+ M(h M (h)+ (-h)-N(h) N'(h) + n=0 -+S M M(-h) + (h+(-h)+ M (h)+ SeI -h) - N(Nh (-h)- N'(-) L> z> z' +S I M(-h) + M(h) '(-h) + M(h)+S eI N(-h) - N -h) - h >z> (2.26) The impedance boundary conditions at z = L are: A 34I *(-g +iknG ) 0 (2.27) a4i 4 z=L l A + ik z. ( —VxG +-VxG ). (2.28) E az 41 n 41 Substituting Eq. (2. 2) and (2. 26) into (2.27) and (2. 28) we obtain ihe:L -e ihL ihL]+ ik r^ihL ihL ihL -ihe h eih + ihL + in ee +S (e -e 0 ihL, ihL ihL + ik eihL -ihL ihL] -ihe +Sm ihe -ihe e +S (e +e 0 MI n L mI The above two equations can be solved for the two unknowns: ihL (h-kn) e eI= 2h cos (hL)- i kn sin(hL) (2.29) ihL (nh-k) e (2.30) mI1 2Lk cos (hL) - i h n sin (hL)J

20 Thus we have obtained the expressions for two approximate dyadic Green's functions by applying the Leontovich impedance boundary condition to one problem. 2.8 Relationship Between Exact and Approximate Green's Function The difference between the exact Green's function and the approximate Green's function is that the former uses the continuous property of electric and magnetic fields at interfaces while the latter imposes the Leontovich boundary conditions at interfaces. Thus, the structure of Eqs. (2. 5) and (2. 22) is identical with the exception that the multiple reflectionl coefficients are different. To determine how close the Leontovich boundary condition approximates the exact boundary condition, it is necessary to study the physical meaning of the impedance condition. The fact that such an impedance exists at the interface between two medium can be clearly seen by considering a plane wave incident on the boundary from the direction of free space. Application of Snell's law shows that the transmitted wave is deflected toward the normal because of the refractive index n. For a fixed direction of incidence, the angle between the direction of the transmitted field and the normal of the interface is in the order of /j which implies that the radial eigenvalue X of the transmitted * j- - 1- - i -- - - -- =- - = _ - field is in the order of sin (ni). For a value of the refractive index I |n| large compared to unity, the eigenvalue X is negligible. Thus hE c k -nk (2.31) The same relationship holds in our problem. When the approximation h ~ n k is put into R and R, simple algebraic manipulation gives rise to E e m the expressions for ReI and R I. Similar arguments hold for the relationship ____ ejl __ ___ between S, S and S, S. Thus, it is clearly seen that the larger the value of n, the closer the two solutions. In other words, when the refractive index is in the order of unity it is not suitable to use impedance boundary condition.

CHAPTER m FIELDS DUE TO HERTZIAN DIPOLES 3.1 Introduction The dyadic Green's functions needed in our problem have been given explicitly in Chapter II. These Green's functions can be used to find an integral form for the fields in the region in between two interfaces for any kind of current source. In general, any current source can be decomposed into elementary Hertzian dipoles. So only the fields due to x-directed and z-directed Hertzian electric and magnetic dipoles are considered. In general, the exact form of the field is rather difficult to find. However, an asymptotic solution can be found using saddle point integration which has been discussed in great detail by many authors [13], [14]. 3. 2 z-directed Dipole 3. 2.1 Integral Representation For a z-directed infinitesimal electric dipole located at a height z above the interface z = 0, the current distribution can be expressed as J(R') = z 6 (x -0) 6 (y' -0 () 6(z' - z )I. (3. 1) The antenna current distribution is assumed to be uniform with a magnitude I. If the length of the antenna is I then in the limit as I goes to zero, the current I goes to infinity and the dipole moment IS. is assumed to have a finite value. Since the Green's functions derived before contain posterior terms of the Mt (+ h) and N' ( +h), the following two relationships, which arise from enX enX - o o the explicit expressions for the M and N functions given in Eq. (2.2), enable us to write: - IEMt (+h)- 7(R) dV'= 0 2.(3.2) rrr A 6X +ihz ' (+h) J(Rt')dV' = I-k -e JJJ en k 21

22 Equation (3. 2) implies that the final expression for E will contain only terms of NO type, because the M functions have no z component. eOX _ ---We recall that the electric field may be found from a knowledge of the current distribution and the Green's function. The specific relationship is: E= iwp O1 G (R/R) -J (RT) dV Substituting Eq. (3.2) into the above equation, we obtain [ -ihz ihz odihz 47Tkf e N h)+N o h);he0 + 4 e e ~ - -I, RI -ihz ihz l + R N (-h) +N (h) e h+e h L>z>z m eOX e >> o (3.3) r h [-ihz ihz] [ ( + R N e(-h)+ Ne (h) e ~+e ~ z >z>0 M eOX\ e OX o - The electric field representation is in the form of a semi-infinite integral. Under certain conditions to be stated later, we can find the asymptotic expression by the method of saddle point integration. To derive this expression we shall first transform the integral as given by (3. 3) from semi-infinite path to an infinite path in the \-plane. Applying Sommerfeld's half-circuit relation [16] between the two kinds of Hankel functions and the circulation relation of Hankel functions: (2) (1) H (Xr) -H ((-Xr) 0 0 to our integral, we obtain

23 00,u I! ^ 47rk -00 where — (1) h) 2 cos (hz Neu)(h)2 os(hz)+ [N.(-h)+ N (hi 2cos(hz ) e o mL eOX euOX N (-h)+ N (h)he +R [eo) (-h)+N(h)12cos(hz ) LeOX eOX m eOX e&\ J (3.4) -N(1 h)- (1) eihz] N (+h)=Vx H (Xr)e- z eOX'- L 0 This is the exact integral representation of the electric field due to a z-directed infinitesimal electric dipole. The final task is to evaluate this integral. 3.2.2 Asymptotic Solution When Xr is large compared to unity, corresponding to a far-zone field, n-(1) -(1) the Hankel function in M and N can be approximated by its asymptotic expansion, 1I (1) 2 - Hn (rr) = r (-i) n n rlr 2 irlr e Accordingly, the functions M) and N(1) become Accordingly, the functions M and N become: ~(1) n+ 2 M(1) (+h1) (-i) en0r o / 2 i(rr+hz)cos, r) -r i. n V trr sin -s n+1 N1) (+h) (-i) en - nr0 o 2 i(r6r+hz) cos rl \ esin V -rre sin n0 Then, the far-zone electric field is represented by: -Wpm II 0 47rk E c (- i X\ 2rX cos(hzo) k e (h)(1+R )+ -o0 +T 9i (At - hz) (A + hr) +R 2cos(hz ) e ) (Xh m o k L> z > z -

24 l fr ( i)/2 X { + 2R cos(hz x xei(r-hz)( +h)+i(r+hz) (X -h,. k i (r - he z ] >z>0 (3.5) k O ~ It is desirable to introduce a new complex variable j3 to execute the integration. This variable is much easier to interpret physically than the eigenvalue ) in the presentation of the wave propagation problem. The technique used here is the same as the one used by Budden [3], Felsen [17] and Nomura [18]. / To facilitate the function-theoretic manipulations involving integrals of the type occuring in Eq. (3. 5), the new variable 13 is introduced through the transformation t = ksing3 which makes \= +k a pair of regular points in the 13 plane rather than two branch points in the X plane, since the transcendental function sin/3 is single valued. The paths of integration in both complex X and 13 planes are shown in Figures 3-1 and 3-2. The reflection coefficient R in Eq. (3. 5) is a function 2 2 22 2 of h = /k -X and hE +n k2- X. Corresponding to the four combina2 2 22 2 tions of signs of + k - X and + n k - X, the integrand is four-valued and its Riemann surface has four sheets. These sheets are connected with 2 2 2 one another by branch cuts along the line Re/k - = 0 and Re k -X =0 from X = +k or X = +n k to infinity as shown in Fig. 3-1. The choice of branchcut is discussed in detail in Appendix A. The branch cut starting at the origin is due to the Hankel function of the first kind. Because of this particular choice / of branch cut, the diagram of Fig. 3-1 is the upper-most sheet among the four /... s s I 2 2 22 Riemann sheets. It is then defined as Re - > 0 and Re k > 0 \/i7-X2>Oand Rein k -; >0,

25 complex X plane 411-~ /, / _/pole brcbranch cut ReX -kIG. 3 PATH OF INTEGTIO N IN COMPLEX X PLANE. -nk ~ FIG. 3 -1; PATH OF INTEGRATION IN COMPLEX \ PLANE. \ k) \ FIG 32: PATH OF INTEGRATION IN COMPLEX PLANE. FIG. 3 —2: PATH OF INTEGRATION IN COMPLEX 9 PLANE.

26 where the shaded region represents that imaginary part of the exponential term which is greater than zero. The transformed diagram is shown in Fig. 3-2. Through the transformation x = ksin/3, the four Riemann sheets in the X domain become four adjacent sections in 13 domain with a period of 27r. The region of concern is that between - 7r/2 and r/2 as shown in Fig. 3-2 which corresponds to the uppermost sheet as shown in Fig. 3-1. The branch cuts starting from x = +k in X domain are shown in dashed lines as regular lines starting from 3 = +7r/2 in 13 domain. Since the integral path C no longer circles around these dashed lines, the introduction of the transformation is now easily realized. There are poles arising from the denominator of the reflection coefficient R. These poles are closely related to the surface wave. A detailed study of m the pole characteristic will be given later. In applying the transformation X = ksin13 to Eq. (3. 5) the counter transformation pair h = +k cos 3 exists. The choice of the sign again follows the rule that the integral must be convergent in order to satisfy the radiation condition. For the exponential term in the integrand having i(Rrr+hz) the positive sign has to be chosen, while for i (Xr -hz), the negative sign has to be used. Thus, the integral in the complex 13 domain becomes W I1 A CO —S o( - 4rk 3k 2 sin 2i cos(kz cosP)(sinBz-cos13r)(l+R)4wk 7Rsin0 o m \ 2 i-i ao A[ A1 ikRcos (0 —) 2 C - R 2 Cos (-kz cos 3) (sin z-cos 1r) e ( rn o i L>z>z () — i 0 (3.6) ____i - _ikz cos - s -cosin P +R c+ o(kz (si k cos 4 rk nRsin 0 e 0 2ikz cosp (kzos1 Fi(sink Z cOr)o +2R cos(+kz cos z >z > 0 -

27 Now the saddle point method is applicable in the complex /3 plane. The path of steepest descent is shown in Fig. 3-2 as CS. The deformed path of integration passed over the saddle point at cos (0-)0 -lr where 0 = tan - z which gives /3= 0 as a saddle point. If the path Cs does not pass near a singularity nor cross singularities when deformed, the asymptotic solution of the integration is W 0I I ikR c E- ^- ~ isin0 — cos(kz cos 0)(l+R - R )x (z), V R I om m - z kcos 0 k cos 0 x (sin0 z-cos 0r) L> z>O (3. 7),-2 h k2. ihL (hkE-h k )e h E E ) -, 2 k~h cos (hL)- i k~hsin (hL) h= +kcos 0 E E h +kcos 0 This is the far-zone electric field solution for the problem of concern. Under our assumption that no singularity is crossed, this far field contains the geometrical optics reflected and direct waves. However, when we deform the path from C to Cs, it is possible that it sweeps through either poles or branch cuts. These singularities contribute the so-called surface waves and lateral waves respectively. They are discussed in more detail in Appendix A. 3.2.3 Comparison with Mode Theory Results Previously mentioned references [2, 3, 5] approach this waveguide problem using mode theory. Section 3.2. 2 mentioned that the poles contribute to the surface wave. Thus, the relationship between the dyadic Green's function

28 technique and the mode theory approach should be pointed out. The characteristic equation for the existence of a self-consistent mode in the mode theory is - R(0) RE () exp(2 i k L cos 0) a 1 (3.8) where 0 is the angle between the plane wave normal and the interface, R(0) is the reflection coefficient at one plate z = 0, and RE(0) is the reflection coefficient cient at the other plate z = L. For the TM waves travelling in the waveguide, the Fresnel reflection coefficients at the perfect conductor plate and the earth plate are given by: R(0) = 1 2 2 2 1/2 n cose- (n -sin 0) Putting the transformation X = k sin 0, h= k cos 0 into the above equations and substituting the results into (3.8) yields n2n - (n2k2-t2)/2 -2ihL 2 22 21/2 n h-(nk -X ) -2ihL. (n2+ (nak2"O2)X2 which is equivalent to hE -- itan(hL) (3.9) n2h where hE (n2k2- X2) /2. The mode equation (3. 9) is identical to the equation of the poles in Eqs. (2.14)

29 and (2. 15). Similarly, for the TE waves, R(0) = -1 2 21/ cos0-(n - sin ) 62 2 2 /9 cos 0+ (n2- sin 0) 2 and the mode equation can be found as h -h itan(hL). (3.10) hE This mode condition is identical to the equation of the poles in Eqs. (2.6) and (2.7). Hence, the poles of the dyadic Green's functions are directly related to waveguide modes. Thus, the discussion of the mode theory only constitutes the pole contribution part of our approach. This important fact implies that when j the pole contribution, which is the surface wave, is not dominant among the waves, mode theory alone cannot describe the wave propagation problem completely. 3. 2.4 Leontovich Approximate Solution When the impedance boundary condition is applied at the lossy dielectric surface, the solution is referred to as the approximate solution. As derived previously, the approximate Green's function has the same form as the exact Green's function with the exception that the multi-reflection coefficient has a different value, so the above procedure can be used without change in this section. There is no need to repeat the entire discussion here. The asymptotic solution is as follows: o/ L It ikR E C --- ~-i sinO e cos(kz cos0)(l+R -R+ ) x kcos 0 -k cos 0 x (sin z -cos ) (3.11) -

30 ihL (k- nh) e where R (k-nh) e mI k cos e 2[i nh sin (hL) -k cos (hL, kcos- +kcos 0 Again this field contains only the direct and reflected waves defined in geometrical optics. 3. 3 x-directed Dipole 3. 3. 1 Exact Solution For an x-directed infinitesimal electric dipole with current moment It and placed at z = z, the current distribution may be written as (R') = x 6 (x - 0) 6 (y -0) 6 (z'- z ) I o Then ihz JJJ X- ni 2 (3.12) +ih ~ +ihz J 1N (+h) J(R) dV' I 1 2 e even parts J - ni k 2 1 nx 1 where 6n n nl= 0 n01 The terms involved in the calculation of the field are of the type M_ 0(+h) and Ne l^(h). It implies that for a horizontal dipole, it can no longer represent the electric field with one vector wave function alone. This result agrees with Sommerfeld's derivation [16]. Thus, the electric field due to a horizontal electric dipole is:

31 00' ( li s; (h)isin (hz +R M -h)+N2isin(hz 4 J ho1 -oJ eel M X J oJ 1{h) 2isin(hsOZ +R VNeXh)+N h F2i sin (hz k ei L e 1X 0] (x L > z > z | H0 J ~LM h1 l -h)-M. (h3 ee +Re [Mol: h) -Mol(h) 3 [ 1 n ihz il ( h -h) - + (h)-M ~2isin(hz 4wj hA 1o olX olx + Nel-h)+lk(h e +RF l -Nh)+NelJ! sin(hz k el e m e ea J[ e o} z >z>0 (3.13) where subscript (x) means that the field is due to an x-directed dipole. The basic structure of Eq. (3.13) is the same as Eq. (3. 6), hence the procedure used for the vertical dipole case is applicable to the horizontal dipole case. After the application of half-circuit relationship, the path of integration changes to an infinite path in the X-plane. Substituting the asymptotic expansion of Hankel function and changing the variable of integration into 3 we can again apply the saddle point method. Thus, the asymptotic solution for an x-directed dipole is found to be wM0 Ie ikR A Ex) - er R sin(kz cos 0)(1-RR + ) sin 00+ X 7r R o e e o kcos e -kcos0 + cos 0 cos sin (kz cos 0) (1( RR ) (sin 0 - cos >O 0 mk cos 0 - k cosO L>z> 0 (3.14) (h-hE) eihL where R = (h- hE^eihL where Re 2 h cos(hL) - i hE sin (hL) h=+kcosO h=+kcos

32 3, 3. 2 Leontovich Approximate Solution Again, following the previous derivation, it is apparent that with the substitution of ReI and R I into (3. 14), we obtain the approximate solution due to a horizontal dipole: wp I.! ikR E( in(kz cos )(1 -R l +R e x)I 7T R o e 0 k cos 0 k cos +i -l- p~ns sin )sin o ~ + + cos 0 cos 0 sin (kz cos 0) (1 +R -Rmi kcos (sinz - cosr) o mmI kcos 0 - k os0 0 (3.15) where R = he h= + k cos 0 ihL (h- nk) e 2[h cos (hL) -i n k sin (hL)] h=+kcos 0 3.4 Horizontal Current Loop As shown by Stratton [20], a convenient measure of the strength of a loop carrying a uniform current is in terms of the magnetic dipole moment. Hence, a small horizontal current loop carrying a current I with area A located at (0, 0, z ) can be represented by a magnetic dipole with moment m(R)= z6(x'-0) 6(y'-0)6(z'-z )IA 0 We recall that the magnetic field can be determined with the aid of the following equation: H'= k2 G(R/R) - m-(R) dV' I The analysis is exactly the same as that for a z-directed electric dipole. With the changing of Green's function and excitation factor in Eq. (3.11) the result is given.

33 ikR NH,- i- sinO R sin (kz cos 0)(l+ S | -S ) (h) w R o elk eelko | ekcos0 e_kcos0 (sin0 -cos r) (3.16) where S is represented by Eq. (2.16). The approximate solution is obtained e by replacing S with Se in Eq. (2. 29). e ein 3. 5 Vertical Current Loop For a small vertical current loop carrying current I with area A located at (0, 0, zo) 0 (R') = x 6 (x - ) 6 (y - 0) 6 (z'- z ) IA 0 Thus: ikR - A ^ - IAk e H - os (kz C 0) (1 - +S ) sin 0 0 + (V) 7T R 0 o m o vk cos R -k os0 +cos 0cos0cos(kz os 0) (1 +S -S ) (sin z- cos r) kcose |-k cos 0 L >z> (3.17) where S is represented by Eq. (2.14). m The replacing of S by S gives the approximate solution. Thus, the asymptotic expressions of the fields due to various simple dipoles are derived. The similarity of the fields of the electric and magnetic dipoles was expected because of the duality principle. The fields patterns and the detailed results will be presented in the next chapter.

CHAPTER IV CALCULATION AND DISCUSSION OF RESULTS 4. 1 General Remarks In this chapter we present a computational study of a particular problem involving the wave propagation above lossy earth. The height L of the particular structure is chosen to be 30 meters. The environment of the ground will tvary from good to poor earth. The refractive indices associated with these grounds at different frequencies will be cited later. The results based on the Leontovich impedance boundary condition at the lossy interface are studied. The errors associated with these approximate conditions are calculated. For a given elementary point source, the excitation factor is investigated. The advantages of using the dyadic Greents function approach are discussed. In the following sections, three different types of curves are presented. They are the reflection coefficient modulus, position of the poles in complex 0 domain and the far-field pattern for different earth material and frequencies. 4. 2 Percentage Error by Using Impedance Boundary Condition At the very beginning of Chapter m, the Leontovich impedance boundary condition is discussed. This condition is valid only under the assumption that n >> 1. However, no quantitative criterion was set for the values of n. The bound of n can not be ascertained in general, since the exact criterion depends on the geometry and the material constant of the earth. In this section we attempt to find the difference between the approximate solution and the exact solution for different values of n. Let us assume that aXE and EE are independent of frequency. The complex index of refraction n is represented as n xE l+i-.. E In general, crE and EE are functions of the water and the salt contents of the earth. Here we categorized earth into two kinds: good earth and poor earth. 34

35 The refractive indices as functions of frequency are given below. frequency used good earth poor earth for calculation 2 4 for calculation (= 10, rE 10-2 /m) (e 4, aE 104 /) r E r E 2 2 f n n 60KHz 10 + i 3000 4 + i30 600 KHz 10 + i300 4 + i3 6 MHz 10 + i30 4 + i0.3 Table 4-1: Refractive Index for Different Frequency In our problem the difference between the exact and approximate solutions is due to the different values of the multi-reflection coefficient. Since the calculations for different types of sources involve different expressions for the reflections it is therefore suitable to plot and discuss the coefficient modulus first. The coefficients R and Rm of Eq. (2.8) and (2. 24) which relate to both electric dipoles can be written in the form: 2 f77 ikLcos (n cos 0 - Vn-sin2 ) e R = m 222 2 [ -sin 0 cos (kL cos 0) - i n cos 0sin (kLcos 0 ikL cos 0. (n cos 0 - 1) eikLcos I2 cos(kLcos 0) - i n cos 0 sin(kLcos 0 The coefficients R and Re of Eqs. (2. 6) and (2.25), relating to the horie eI _' _ _ -...' zontal electric dipole can be written in the form I \2 " ikL cose _ (cos 0- V- sin e) e cos 2 [cos (kLcos 0) - i i -sn sin (kL cos 0 ikL cos 0 *R = -(cos 0 - n) e eI2 — os cos(kLcos )- sin(kLcos 2 Ios 0 cos (kL cos 0) - i n sin (kL cos

36 Before showing the curves for the reflection coefficients, it is worthwhile to mention that the cut-off frequency of the guide is 5 MHz, if both plates are perfectly conducting. Thus, 60 KHz is far below this cut-off value. The spatial region of concern is the far zone, thus the direction of observation 9 is greater than 82 degrees and less than 90 degrees, hence, the significant range of cos 0 is extremely small. Studying Figs. 4-1, 4-2, and 4-3 it is found that the reflection coefficient R line gets more perpendicular with both interfaces as the operating frem -- quency gets lower. At low frequencies the earth acts like a perfect conductor implying that the two interfaces look much alike. Thus the R curves are m more symmetrical with respect to both interfaces. The difference in R Is for __m good earth and poor earth becomes smaller as the frequency is lowered. / However, these properties do not apply to R. The reason is that there is a e pole at 0 = 7r/2 as shown in the figures. The difference in the R es gets much e larger as the frequency is lowered. Even though the earth surface has a veiy large refractive index as is the case for f= 60 KHz (Fig. 4-3), there is no symmetry property. Thus, the geometry of the structure as well as the nature of the material affects the results. The following table gives the percentage error in the results when using the impedance boundary condition. f n 2 cos 0 0o error of R 0 /0 error of R 6 MHz 4 + i 0.3 0.039 1.38 0.84 6 MHz 10+i30 0.039 0.35 0.07 600 KHz 4 + i3 0.004 0.092 6.925 600 KHz 10 + i300 0.004 0.0084 0.016 60 KHz 4 + i30 0.0004 0.002 0.358 60 KHz 10 + i3000 0.0004 0.0 0.0041 Table 4 —2: Percentage Error of Reflection Coefficients

z(m) earth good earth 10-1 - poor earth R m 1.6 "~ 4 - \~ \ 101,2 RI 2 \ e FIG. 4-1: MAGNITUDE OF REFLECTION good COEFFICIENT AT 6 MHz..0 poor 8 6 4 - mag. of R 2 and R 2- e 0 1 2 3 4 5 6 7 8 9

I 24 -I I 1 o-1 20 R good k.-poor 16 12 \ FIG. 4-2: MAGNITUDE OF REFLECTION 8 \ Re COEFFICIENT AT 600 KHz. good poor mag. of Rm -— ____ and R e CO3 00 0 I I I 0 1 2 3 4 5 6 7 8 9

CD FIG. 4-3: MAGNITUDE OF REFLECTION COEFFICIENT AT 60 KHz. R x104 I- e -- mag. of R —.and R onnd htonr, I I I e m 0 1 2 3 4 5 6 7 8 9

40 The percentage error for R I at n - 4 + i 3 is very large. The reason el for this large value is the fact that the reflection coefficient Re has a pole near this particular combination of n and cos 0. While for the same combination the approximate function ReI has the corresponding pole not as close as the exact one. Thus, for all points of observation near the poles of Re v the error of using the impedance condition is large. Other than those particular points in general for 1n 2 > 5, the error is less than 1 percent. For very large values of n the difference is negligible. The far-zone field is directly proportional to these coefficients so that the error in the far field in using RMI and R is of the same order of magnitude. They will be shown in the next section. 4. 3 Far-Zone Field The far-zone electric fields due to the vertical and horizontal electric dipoles are presented in this section. There are the Er and E components corresponding to the vertical dipole and only the E0 component corresponding to the horizontal dipole. They are shown in the same figure (with only 0 = 90 0 shown for the horizontal dipole case). Due to the use of the saddle point method, the fields obtained are in the far-zone. The definition of the far field in our research is r > 15X. Thus, the range of the far-zone varies with the frequency. Hence the magnitudes of the fields are not in the same order of magnitude for three frequencies. All curves of E, E and E tilted forward near the dielectric surface r z 0 compared to that at the conducting plate. Figures 4-4 and 4-7 show that E is stronger for the poor earth than for the good earth. As the frequency is lowered the difference in E0 with good earth and poor earth increases greatly. But the differences in E and E get smaller. This is directly related to the properties of R and Re as discussed before., Figures 4-5, 4-6 clearly show that the curves of Ez and E are almost perpendicular to both plates. In these two cases, the index of refraction of the earth is so large that it appears to be a perfect conductor.

24 20 16 12 4^ FAR-ZONE ELECTRIC FIELD OF UNIT DIPOLE MOMENT AT 6 MHz.

FAR-ZONE ELECTRIC IFIELD OF UNIT DIPOLE MOMENT AT 600 KHz.

good 24 20 16 12 8 4 10-8 E E r z 10-9 -good z =15m 0 FIG. 4-6: FAR-ZONE ELECTRIC FIELD OF UNIT DIPOLI MOMENT AT 60 KHz. mag( E mag. ( /m) 0 1 2 3 4 5 6 8 9

/ 'good z =5m 0 16 12 FIG. 4-7: FAR-ZONE ELECTRIC FIELD OF UNIT DIPOLE MOMENT AT 6 MHz. (mag. (Vm) 5 6 7 8 9

24 20 16 12 8 4 FcPh 0 V0 I 2 3 1234 5 6 7891 7 8 9 10

46 In Fig. 4-7, the dipole source is located near the perfect conductor at z = 5m, all three electric field curves tilted forward more than in the case of a symmetrically located source. Figure 4-8 shows the results when the source is located near the dielectric surface. A much stronger E is found. Since the separation distance L between the two plates is very small, there is no drastic change in the far field patterns as shown in the figures for different source locations. The percentage error in the far field when using the impedance boundary condition is given below. 2 ~/o error /O error /O error. C ofE ofE of Er z0 6 MHz 4+i 0.3 0.039 1.140 1.122 0.232 6MHz 10+i 30 0.039 0.240 0.221 0.065 600KHz 4+i 3 0.004 0.091 0.093 7.417 600KHz 10+ 4i300 0.004 0.0037 0.0027 0.0157 60 KHz 4+ i 30 0.0004 0.0019 0.0078 0.357 60 KHz 10+i3000 0.0004 0.0 0. 0 0.0043 Table 4-3: Percentage Error for Far Field 4.4 Excitation Factor for Various Sources Several difficulties have to be treated carefully when the problem involves lossy media. They are: how many components of the Hertz vectors are responsible for the field and what is the excitation factor for a given source. The advantages of using the dyadic Green's function technique is that it systematically solves electromagnetic problems in such a way that these questions are easily answered. The dyadic Green's functions are chosen so that when applying the vector Green's theorem to the region concerned, the surface integral equals zero as a result of the boundary conditions and/or the radiation condition. Thus, the two field equations are

47 Ex i Wu|||G (R/R) J T(t) dVt H k fG4(R/R) m() )dVt. The dot product in the integrand shows explicitly which Hertz vector is involved and what is the excitation factor. For example, in Eq. (3. 3), it is / clear that for a z-directed electric dipole, only the N function is needed to specify the field. This is in agreement with Sommerfeldts result; however, he used physical considerations to find this particular vector function. For an x-directed electric dipole, it is seen from Eq. (3.13) that the N.. and M 1X _______olX elX functions, both of higher order, are necessary. This again agrees with Sommerfeld's discussion on the choice of the Hertz vectors for a horizontal dipole. For a detailed discussion, see Chapter 8 of reference [6]. As discussed above, the excitation factor for a given source is built into the Greents function. Since the radiated wave is composed of several different component waves, the excitation factor has several components. These components are directly related to the results of the integration as derived before. The excitation factors of the electric field for different kinds of sources are shown in Table 4-4. The pole location, e +iOi, is discussed in the next section. For completeness, the propagation factor of the surface wave is derived in Eq. (A-1). Thus, both difficulties are solved at the same time. These features, which avoid all the determination works on how to choose the potential functions and the excitation factors, show that the dyadic Green's function technique is more systematic.

source type direct and reflected wave surface wave (mode factor) electric dipole i o0) -i If k sin 2(0R+i.) sin ckz os(O i0)j x xrsin0+0co0cos000 /o2i R LkZ ( R | l x/ — 8in0 0u+cos(0R+ii) cos00J xV 7Tr current loop 4 sinsin(zco) 0 iIAk 2 sin (0 +iO) cos cos(4 i x-directed W~/aoIl 3/ Ixk -0 I vertic dipole Ak 0) x i ksin2sinkz os(+ x current loop osk csin os sin o (R ii x sin+cosocos0 sinJ in0 + cos (0R+i.) cos 0 7 r R 1 horizontal IAk2 3 ' ' sin 0 sin (kk sinkz cose)x current loop i k 5 s Rios R 1 current loop 4 k1 R 1 x in 0 + cos e cos 0 A 2 i i~no + cos(eR+ ie.) Cos Table 4-4: TABLE OF EXCITATION FACTOR FOR WAVE COMPONENTS FOR DIFFERENT SOURCES.

49 4. 5 Distribution of Poles in the Complex 0 Domain When considering wave propagation, it is necessary to discuss the properties of mode poles. For guides with perfectly conducting walls, the poles are either on the real axis of the complex 0 domain, corresponding to propagating modes, or on the imaginary axis, corresponding to evanescent modes. These poles can be found from the fundamental equation of mode theory. For example, the equation satisfied by TM waves is kL cos 0 = + mr where m = integer. When one plate of the waveguide is not perfectly conducting, the mode equation is represented by Eq. (3.9) h = =itan(hL) nh By using Newton Raphson's method (Appendix B), the solutions for the propagation constants of quasi-evanescent modes can be obtained. The impedance boundary condition is used to get an approximate solution for the modes in Eq. (3. 10) which yields - itan(hL) nk The locations of these poles are also obtained here. They are shown in the following table for comparison. In general, for the same operating frequency, the higher the value of the refractive index n the closer the approximate solution is to the exact solution. This points out again the necessary assumption for the proper application of the impedance boundary condition. In the case of the poor earth at 60 KHz (n = 10+i 3000), due to the large value of n, it is logical to assume that the poles will be located very close to those of the perfectly conducting case. However, the imaginary part of the pole is very large compared to zero as shown in Table 4-4. Notice that particular frequency is very far from the cut-off frequency, 5 MHz in our study. Thus, the electric length between the guides plays an important role as well as the refractive index n.

50 frequency mode complex 0 (eR+ i.) number perfect conductor exact approximate 6 MHz 0 1.5708+i0 1. 447+i 0. 1930 1.448 +i 0. 1949 (n 10+i 30) 1 0.585+i 0 0.6457 + i 0.08229 0.6463 + i 0.0822 2 0+i 1.094 0. 0176+i 1.086 0.0178+i 1. 086 600 KHz 0 1.5708+i 0 1.425 + i 0. 3568 1.426 +i 0. 3569 n =10+i300) 1 0+i2.81 0.00144+i2. 808 0. 0016+i2. 808 60 KHz 0 1. 5708+i 0 1.346+i 0.619 1.346+i 0. 619 (n 10+li 3000) Table 4-4: Location of Poles for Good Earth (a = 102 /m) frequency mode complex 0 (0 t- ii) number perfect conductor exact approximate 6 MHz 0 1.5708+i0 1. 289+i0.1713 1.267+i 0.169 2 n =4+i0.3) 1 0.585+i0 0. 6021+i 0.2817 0. 5913+i0.2927 2 0 + i 1. 094 0.0745+i 1. 101 0.0611+i 1. 102 600 KHz 0 1.5708+i 0 1.106+i 0.8015 1.078+i 0.8752 - ' - -"" -~(n =4+i3) 1 0+ i2.81 0. 0575+i2. 775 0. 0166+i2.805 60 KHz 0 1.5708 + i0 1. 153+i 1.481 1. 180+i 1.514 n24+30) (n = 4 + i 30) A-4 Table 4-5: Location of Poles for Poor Earth (a= 10 P/m) Since the path of integration in the complex 0 plane is deformed when evaluating the integral, it is convenient to plot the pole distribution in the same plane for the discussion of surface waves. The pole locations for our problem are shown in Figs. 4-9 and 4-10. Note that the number 1 pole for the good earth case in Fig. 4-9 has a much smaller imaginary part than that of the number 0 pole. The opposite is true for the poor earth case as shown in Fig. 4-10. The steepest descent line I-!

ImO 1.3 1.2 1.11 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 6 MHz * 600 KHz * 60 KHz o perfect conductor at 6 MHz number 2 number 2 01 * number 0 FIG. 4-9: POLE DISTRIBUTION FOR GOOD EARTH. * number 0 \ number 0 \ s0 steepest descent line \ number 0 I I I I. I - /I I Ree number 1 I I I I I \number 1 J K I i I * * 1... us 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7

* number 0 ImO 1.3 1.2 1.1 1.01 0.9 0.8, 0.7 0.6 0.5 0.4 0.3 0.2 0.2 0.1 0o o 6 MHz * 600 KHz * 60KHz ~ perfect conductor at 6 MHz o number 2 number 2 * number 0 FIG. 4-10: POLE DISTRIBUTION FOR POOR EARTH. en tO number 1 number 1 I I - I l - I I \ I \ /steepest descent line number 0 o 0 1 - I 1 I I - I /2 m -. -.... ReO 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 00.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5

53 0 has an inclination of 45 with the real axis for any point of observation. Thus, when the observation point is beyond 0 41. 8, for the good earth case, the contribution of the number 1 pole has to be taken into account. However, since the number 0 pole is on the right hand side of the steepest descent line, this pole contributes nothing to the field for any observation angle. Thus, the strength of the surface wave depends only on the imaginary part of pole number one. The dependence of the surface wave upon the location of the pole is discussed in Appendix C. In the case of poor earth, both poles, number 0 and number 1, enter into the picture when the observation angle is larger than 49. 80 and less than 900. Because the imaginary part of the number 0 pole is much smaller than that of number 1, it is evident that the number 0 mode is the dominant one. However, if the observation angle is greater than 49. 80 and less than 83. 60, only mode number 1 contributes, so it is the dominant mode. Thus one has to be careful in discussing wave propagation when lossy media are involved. The zeroth order mode does not always supply the dominant contribution to the surface wave as suggested by some authors. It is possible that during the deformation of the integration path no pole is swept through until 0 has a certain value. Then it is clearly seen that no significant surface waves exist when the observation point is less than certain angle. For example, for 0 < 410 in our two cases, only the direct and the reflected waves exist even in the near zone.

CHAPTER V CONCLUSIONS AND RECOMMENDATIONS 5. 1 Conclusions In this thesis we have developed the technique of using the dyadic Green's function method to solve a particular problem in wave propagation. We have found the exact integral form of the fields due to various simple Hertzian dipoles in between two parallel plates with different material and the asymptotic solution associated with them. The use of the Leontovich impedance boundary condition at the lossy interface is discussed. The percentage error in using this condition compared to the use of the exact boundary condition is calculated. Except when the point of observation is near a mode pole, the higher the refractive index the lower the percentage of error. The pole contribution of the Green's function technique is shown to be the same as the mode contribution of the waveguide theory. A detailed study on the pole distribution in complex 0 domain leads to the conclusion that it is premature to predict which mode is dominant before actually examining the pole location. The symmetry property of the dyadic Green's function in a two layer problem is proved. The same procedure could easily be extended to other many-body systems. 5. 2 Areas for Future Study In the present investigation the exact integral formulation of a parallel plate waveguide problem is derived. The far-zone field is calculated by the saddle point method, which is only the asymptotic expression. When the time domain solution is needed, the high frequency, intermediate frequency and the low frequency solutions are required for the application of Fourier transform. Thus, the exact integration for all frequencies needs to be derived. The technique of performing this infinite integral is an interesting problem and worth further investigation. 54

55 Another area for future work would be the efficiency problem. According to the study in section 4. 3., the far-zone field is slightly affected by the source location z. However, the portion of energy entering into the lossy dielectric is directly related to the source location It is very useful for the simulation test designer to understand the distribution of power and the efficiency of a simulator L21].

REFERENCES [1] Watson, G. N. (1919), "The transmission of electric waves round the earth," Proc. Roy. Soc., 95, p. 546. [2] Bremmer, H. (1949), Terrestrial Radio Waves, Elsevier, Amsterdam, Netherlands. [3] Budden, K. G. (1961), The Wave-Guide Mode Theory of Wave Propagation. Prentice Hall, Englewood Cliffs, N. J. [4] Budden, K. G. (1962), Radio Waves in the Ionosphere, Cambridge University Press, London. [5] Wait, J. R. (1962), Electromagnetic Waves in Stratified Media, Permagon Press, New York. [6 Tai, Chen-To (1971), Dyadic Green's Functions in Electromagnetic Theory, International Textbook Company, Scranton, Penna. [7] Tai, Chen-To (1973), "On the eigenfunction expansion of dyadic Green's functions, " IEEE Proceedings, No. 4, p. 480. [8] Schwinger, J. (1943), "Fourier transform solution of integral equations," M. I.T. Radiation Laboratory Report 43, 44. [9] Levine, H. and J. Schwinger (1950), "On the theory of electromagnetic wave diffraction by an aperture in an infinite plane conducting screen," Comm. Pure and Appl. Math., m, No. 4, p. 355. [10] Morse P. M. and H. Feshbach (1953), Methods of Theoretical Physics II, McGraw-Hill, New York, chap. 13. [11i Tai, Chen —To (1954), "A glossary of dyadic Green's functions, " Stanford Research Institute Technical Report No. 46. [12] Fock, V. (1946), "The field of a plane wave near the surface of a conducting body, " Jour. of Phys., Vol. X, 5, p. 399. [13] Leontovich, M. A. (1948), Investigations of Propagation of Radiowaves. Moscow. [14] Feynberg, Ye. L. (1961), The Propagation of Radio Waves Along the Surface of the Earth, Clearinghouse Translation 1967. [15] Hiatt, R. E., T. B. A. Senior and V. H. Weston (1960), "Studies in radar cross sections XL-surface roughness and impedance boundary conditions, "The University of Michigan Radiation Laboratory Report 2500-2-T. [16] Sommerfeld, A. (1949), Partial Differential Equations, Academic Press, New York, p. 197. 56

57 [17] Felsen, L. B. and N. Marcuvitz (1973), Radiation and Scattering of Waves, Prentice Hall, Englewood Cliffs, N. J., p. 465. [18] Nomura, Y. (1953), "On the theory of propagation of electric waves over a plane surface of the homogeneous earth," Reports of the Research Institute of Electrical Communication, Tohoku University, Japan, Vol. 5, No. 3, 4, p. 203. [19] Van der Waerden, B. L. (1951), "On the method of saddle points, " Applied Sci. Res., B2, p. 33. [20] Stratton, J. A. (1941), Electromagnetic Theory, McGraw-Hill, New York, p. 235. [21] Hansen, P. M. (1970), "The Radiation Efficiency of a Dipole Antenna Above an Imperfectly Conducting Earth, " Ph. D. Thesis, University of Michigan. __ [22] Carnahan, B., H. A. Luther and J. 0. Wilkes (1969), Applied Numerical Methods, John Wiley and Sons, New York, p. 319.

APPENDIX A GENERAL DISCUSSION ON SADDLE POINT INTEGRATION In Chapter m, we discussed the situation where the path of steepest descent does not cross over a singularity. In general, when we deform the integration path to that of the steepest descent, it sweeps through some poles and branch cuts. In this appendix, we discuss the general characteristics of the field distribution due to these two kinds of singularities. It is found that the contribution to the field of these singularities under our assumptions is negligible provided that kr is greater than 207r. Since the saddle point integration always involves multi-valued functions it is worthwhile to consider several often neglected points about branch cut integrals and contour integration./ We follow closely the method discussed in reference 17 2 The reason we chose the branch cut along the lines of Re /k -X 0 22 2 and Re n k - X2 0 in section 3.2. 1 was that we could have unique analytic property in each quadrant of the complex X plane. The analytic property of 2k- X on the uppermost Riemann sheet in the complex X domain is.. a~^:^ i X v SSAd- -ai- f vI 1'L. 58

59 Recall that the integral converges only when the imaginary part is greater than zero. Thus, it is clear that the path of integration must pass through the second and fourth quadrants only. There is no way that we can close the contour along a semi-infinite circle; hence the Cauchy residue theorem can never be applied. In the discussion of multi-valued functions, in general, the roots of one equation have to be checked to assure that some roots are indeed the poles of / that equation. The advantage of this kind of choice of the branch cut is to insure that the roots of the reflections coefficient equation are indeed the poles on this particular Riemann sheet [18]. An alternative way of choosing the branch 22 22 cut is along the line of Im k - X = 0. The analytic properties of /- in the complex X domain are shown in Fig. A-2. Riemann sheet, the integration path could be everywhere in that sheet. Hence the contour can be closed and the Cauchy residue theorem so applied. In -;V - V -E 04i, >8 FTG. A-2- BRANCH CUT ALONG I k2-? 0. Again, each quadrant has the unique property in its whole quadrant range. Since the imaginary part of the function is greater than zero everywhere on this Riemann sheet, the integration path could be everywhere in that sheet. Hence the contour can be closed and the Cauchy residue theorem so applied. In

60 general, the roots of the equation will be the poles of that equation on this Riemann sheet 18. In some of the literature, the choice of the branch cut is parallel to the imaginary axis. This choice is neither along the real part being zero nor along the imaginary part being zero. The analytic properties in the complex X plane are shown in Fig. A-3..P. - 73^ ^ <:77-::-* l.':Z.^ 1 A3 A FIG. A-3: BRANCH CUT PARALLEL TO IMAGINARY AXIS. There are two regions corresponding to the first and the third quadrant which exhibit different analytic properties. Thus, again the closure of the contour j can not be used. Similarly, the transformation X = k sin3 maps this mixed property into the complex,3 plane. Hence, special attention has to be given to the region of convergence of the integral. Thus, there is no advantage at all for this particular choice. We recall in Section 3. 2. 1 that when the path of integration is deformed from C to C in the complex 3 domain, it is possible to cross some singularities. In this case, the path is plotted as follows:

61 FIG. A-4: INTEGRATION PATH IN COMPLEX /3 PLANE. Rewriting the integration in the complex /3 domain for the case of a vertical electric dipole: 7 E - dj3 k sin S 2cos(kz cos3) i (l+ S kc -S Ik 47rk J o Rsin 0 k cos 3 -kcosJ 7si ~-+ i cO 2 ikRcos (0-s) As - co x e (Z sin/3- r cos 3) 2 k2 ihL (h k -k h )e where S= E 2 k h cos (hL) - i kE h sin (hL)

62 2 r2 s I S 2 n cos cos(kLcos3)- iIVn -sin sin(kLcos ) cosn- ksin3 cos (kL cos ) - i n cos sin (kLcos3) The roots of the coefficient equation which are the poles of this equation on this Riemann sheet could be found by using the Newton-Raphson's method [22]. The pole location of the problem concerned was shown in Figs. 4-9 and 4-10. The numerical value of the pole is represented by 0 +ie. Thus, the contribution to the residue by the coefficient equation is: r2.~2 2 Res 1S -k os n cos cos(kL cos ) - ivnsin sin (kL cos 3) Re kcs/-osp-kcos/cs d r2 i ( n sin3) cos (kL cosl)-i n2cos sin (kL cosp P=, +ie9 j3 = eR+ ie. -------- ---—. JK The field due to the contribution of the pole is: E~ iwpoI k sin (0R+ i )cos kz cos (R+ i, P R i L c R isin x +Res lcost S-k cos z sin ( i i) rcos (0R+ i i) }kR 1 R i!. e ikRcos ( - R -iW.) R /oe From the exponential function: i Ros(- -ikR cos(- -i)= ikR cscos 0+sin sin) cosh 0. + lii L R + i (sin 0 cos R- cos 0 sinR) sin h] (A. 1) we find that the decaying factor is given by: Exp EkR sin h.i (sin 0 cos - cos 0 sin. (A. 2) From the discussion in section 4. 5, only one pole contributes significantly for the good earth case at 6 MHz. The other poles have a large imaginary

63 part corresponding to highly attenuated waves. To calculate the surface wave field due to the pole contribution, the values are read from Fig. 4-9 as 0 85 0, 0 +iO 0.64+i0.08. In the far-zone region, the minimum distance R i is 10X away, which implies that kR = 62. 8. Substituting these figures into (A. 1) gives the decaying factor of the order 1/55, which is less than two percent of the total field. For the case of good earth at 600 KHz, the number 0 pole is never crossed during the deformation of the path of integration. Only those poles with large imaginary part are crossed, which implies that these waves never travel too far. As shown in Fig. A-4, the deformed path of integration could sweep through a branch cut. It simply means that the path of integration has to travel along the cut to infinity and cross the cut onto the second Riemann sheet. Since this branch cut corresponds to the singularity at nk, which does not occur in the exponential term of the integrand. This function has identical decay characteristics on either Riemann sheet. The contribution due to the branch cut depends greatly on the value of n and kr. According to Felsen 17, the contribution behaves in the order of (kR) 3/2. For the problem with which are are concerned, it is of the order 1/490 in the far-zone. That means that the lateral wave field due to the branch cut is smaller than the surface wave field; hence, it too can be neglected. The previous discussion holds only for the far-zone field. If the nearzone field is considered, then some poles may contribute significantly to the field. The same applies to the contribution by the branch cut. When a pole singularity is near the saddle point of the integrand, the modified saddle point integration method has to be adopted. For detailed information, see reference [17] and [19].

APPENDIX B NEWTON-RAPHSON'S METHOD Newton-Raphson's method, an iterative method of solving nonlinear equations, is used here to find the zeroes of a complex valued function. To find the zeroes of the following equation F(z)=O where zxx+iy we separate the equation into real and imaginary parts. F (x, y) x fl(x, y) + i f2(x, y) X 0 where fl and f2 are real functions. Then we arbitrarily choose a pair of trial solutions x, y, and substitute them into fl and f2. After this iteration, a pair of increments Ax, Ay results, which is given by f -f f1 af2 Ax2 ---l — f ' f 2 y ay af2 af1 1 ax 2 ax A y D where af af2 afl af2 ax ay -y ax _ ~ __-.. --- —--. — --. - - Then the next approximations are: x X+Ax y1 y+Ay Substituting these values into the original equations f and f2, we can proceed in a similar fashion to determine the higher order approximation. ' 64

65 Under certain restrictions [22], both functions f and f will converge to a certain preassigned small value. Then, this is the pair of approximate solutions of the complex equation we wish to solve. The programming of this method to our problem in the FORTRAN language is given in Appendix D.

APPENDIX C DISCUSSION OF POLE DISTRIBUTION In mode theory, the value of 0 for a given mode can be represented as a point in the complex 0 plane. The real and imaginary components of 0 determine the phase velocity and the attenuation rate of the surface wave as shown in /Appendix A. More importantly, the location of the pole can be used to determine whether it; would be swept through during the deformation of the integration path. When a pole is swept through, there arises a Heaviside function which equals unity, implying that the pole contribution should be taken into account. Thus, the location of the pole can lead us to the delimited regions where the pole contribution could be observed. In this Appendix, we discuss the pole locus of the problem of wave propagation between two parallel plates, one of which is perfectly conducting. There are two loci of poles presented for the TM waves. For each case, the frequency/ and the separation distance are kept constant. _ Figure C-1 shows the particular case at; cut-off frequency corresponding to a guide formed by two perfectly conducting plates. The dotted line represents the case when one plate represents the surface of a pure dielectric media,'whose refractive index has zero imaginary component. The broken line represents the case where the real part of the refractive index is very close to unity which corresponds to free space. The direction of the arrow indicates that the imaginary part of the refractive index is decreasing. In the limit both curves end up at n = 1+ i 0, which is the free space case. This case implies the absence of the second plate, hence, an open region. Thus, there is no mode occuring when a wave propagates above an infinite perfectly conducting plate. It is clearly seen that the locus of the pole, when the other boundary represents an arbitrary kind of material, has to be inside the range defined by the above two limits. For number 0 mode, there is a turning point such 66

Ime 1.4 -1.3 / number 2 1.2 1.1 1.0 0.9 0.8 0.7 - number 1 0.6 -0.5 -0.4 poor earth number 0 0.3 / 0.2 /good earth 0.1 01 —/ReG 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 FIG. C-1: POLE LOCUS AT CUT-OFF FREQUENCY.

68 that the attenuation rate stops increasing as the refractive index becomes smaller. While the number 1 mode, on the contrary, rises rapidly as the refractive index decreases. Figure C-2 presents the case for a frequency larger than cut-off. Except for the displacement of the number 1 pole, the shape of the pole locus remains almost the same. Each locus behaves in a similar fashion as before. It is obvious that for a frequency much larger than cut-off, there are more than two pole loci located near the real axis. The number 2 or higher numbered pole will behave like number 1, which has no turning point. For frequency less than cut-off, there exists only the number 0 pole near the real axis around 9 = v/2. Special attention should be given to this particular pole because sometimes it will never be swept through, as the case of good earth at 6 MHz, hence there is no surface wave related to it. In the literature, there is the common concept that the lower order mode dominates, which must be stated with care. In general, for a fixed frequency, the number 1 mode will dominate for large refractive index material. However, when the refractive index is smaller than a cerfain value, the number 0 mode will dominate. It is necessary to study the locus of the poles to assure a proper understanding and correct application of the mode theory.

Ime 1.3 1.2 1.1 1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 2 o0 1 number 0 poor I ReO 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 FIG. C-2: POLE LOCUS ABOVE CUT-OFF FREQUENCY.

APPENDIX D COMPUTER PROGRAM C PROIGRAM FOR THF FAR-7ONE FIELD CALCULATION C \1 *-:c2=(AtiA)=REFRFACTIVF INDEX C 1_ =SFPtARATION DI STANCE!AFTWFEFN TWO INTERFACES CK= K=WAVF NIIMRER C, U=PERMEALFAILITY IN 4RFF SPACEF C 7fl=POSITION OF THF SOU(RCE C ' ~ R=fIRSERVAT IO(N POI NT I N R-DIlRECT IO(N C ~*~W=FREOUJENCY UJSED C ' FZ=FAR FIELD IN THE 7-DIRECTION C F~ R=FAR FIELD IN THE R-DIRECTInN C '' EMZ_=AM4PLITIIDEF OF EZ C EMR=AMiPLIT1J1DE OF ER C *EPHI=F:AR FIELD IN THE PHI-DIRECTION(ONLY FOR PHI= C 90 DEGREES) C F4 MPHI:=AMPLITIJDF OF EPHI C AE7=AP3PROXIMATE FAR FIELD IN THE 7-DIRECTION C * AER=AP3PROXIMATE FAR FIELD IN THE R-DIRECTION C -*AEPHI:=APPROXIMATE FAR FIELD IN THE PHI-DIRECTION C * RM-=REFLECTION COEFFICIENT C RF=REFLECTiON COEFFICIENT C *RMI AP3PROXIMATE REFLECTION C * REI=AiPPROxIMATE REFLECTION REAL K, LKLKZ0 CO0M PLE FX P9S M I JL, N E RE ZE CA SMWYE P HIAE RAE ZRMIJ LAR M *R 1,vR MI,9E X PRIR ER EIE X PCAE PHIPE C PI=3.1415926 A=4. i1=0.3 K=0. 125663 W = 6 * E6 `2 * 4P I R=750. L = 30. ti=. E —7/ (4.*P I K<L =K * L COS T=Wl"tJ/P I *1 0..*,4 R 2 =R **2 N=CSOR T (CMPLX (AR)) Donl 101 LL= 1,5,4 70=5. 'LLL K70= K *70 f)O 10 M=1,31 L L = m- -1 RR=SO.RT(R2+LL**2) RRKK=K*RR EXPR=C'MPLX(COS(RRKK),SINIRRKK)) C =.L L/ RR 70

71 S =R /RR C. KL K LC CC=C0SK(KL) SC=SI N(CKL) FXPC,=CMPLX (CC,SC) Cz=Cfls( K ZCV0,C PCS OR T(CMP LX (A -S c2d3) R - 2 PC CM Itsc) RMI =E-XPC,,: ( N)'cC-1 * 2 **(CC-CMPLX (0. 1 *:N:C~SC RMMI=CAHS(RM) R M 4 I =C A BS (R MT) WR I T ( 6,?20) 20 FO)RMAT('-') WR ITE (6?51 )R MRMI WRITE(6,102)RMM,qRMMI Jt'(L=(CMPLX(AR)*C"'CC-CMPLX(0.,l.P*P*SC)/(P)-CC.-CMP *LX (.-f9 A) 'CC'SC) SMMI4JL=C ABS ( SMJL) IR I TE (6 t107 ) SMIJLSMMtJL AMStM=CAIAS ( ASM) W-R I TF- (60,1 08 ) ASM,AMSM [F(I.EO.0)GI) T10 13 RFI=(C-N-)):,Expc"/(2.):((C*CC-GMpL-X(O,1.Tl,)**SC)) RMFI=CA~S (RET W-RI TE(6952)REREI W,-RITE (6j102 )RME,tRMEI RMIJL-=(CMPLX(0,,-1.)*C"'SC+P*CC)/(C*CC-CMPLX(0.,91.) RPV M IJL =C AFRS ( RMU JL) I4"RJ T ( 6,I 105 ) R M 1)L I R MM (JL 105 FO0R MA T ('SUtM O F R E'3F.1 2.,5) AMRM=CARS (ARM) WR I TE ( 6,106 ) ARMv AMRM PFC,=CU-S T*EXPR *S5 I N (KZO*C U RR EPHI =PEC-', ( i+RMIiL) EMPHI=CAFBS (EPHI AEPHI =PEC* ( 1+ARM) AFMPHI=CARS(AEPHI I WR I TE (6,101 ) EPHI,AEPHI WR I TE (6, 102 )EMPHI,AEMPHI I/

72 113 COn1T I NitIF 107 FORM AT( SHM (IF RM ',3EF12.5) 108 F)RM AT(IA SM I i3FE1?25) FC =COS T"'S"CZ "'CMP LX (0.,1.*FXPR/ RR E?.-EC I~ l+SMJL) FMZ=CABS( EZ) F MR C A S ( R) AFR FC. "C:c ( Il+A SM) AFM1Z=,Af1S ( AEZ) AEMR:CARS(4ER) 100 FORMAT (1.XqI2,2X,R8F12.5) 4R I TF(A6,102 )FMR IFMZ, qAEMR, AEMZ 102 FORMA T( 8X,E12.*5,w1OX? El12.5 v10XE12 5,912XtE 12 5) 51 'FV1R MAT( IMUL R EF c OF FF RM',4E12.5) 52 Fn RMAT( 'MUL REF COEFF RE',4E1295) 101 FnRMIAT('MEPHI 1,4EI2,5) 10 WRI TE(6,lIOO )LLEREZ,AERPAEZ STOP END

73 C, TH IS Pk,R GA M IsII u(G FS T F1)T I RF IE ~JS Fi (Th THE TFRMINA C L FHR THE SEARCHING; (IF PO LES 1IN THF COMfPLEX THFTA C ~PLANIF FOR 'TM' WA\VFS SINCEIPPRCHCEFS) CLOTIOINS MAY NOT CONfVFRGFEF MRPRC(IF C K=WAVP NljjMERR C. L' LSEPAR AT ION DI STANCE fAFIWEEFN TWO0 I NTERFACES C N "*N 2= ( APJ =RFEFRACTI VE I NDEX C I- LED=EI1GENVALIJE IN THE X-DIRECTION C **HL-=EIGENVALUJE IN I4-IE-%Z-DIRFCTION REAl.- KL COnMPLEX L-ED,,HL A=l.02 K=0. 125663 L-= 30. ARIAI TARI~ TRARILY CHOOSE A PAIR OF SOLUTIONS 20 k EADn( i lO0O)XY AL =K*:I. 1 SX= SIN(~lX) SHY =Si NH (Y) cx=COS C X) CHY=CO)SH (Y) C.S =CX'SHY S C =S X -- C, H Y C. C, C, CX C H Y S S =S X,SHY C= A-S, (*'-2 +C S 2 P=SnOR 1T(C-:-*2+)**2) PM=SOR.T(P-C) PP=SOR, T( P+C) PMO1=1./PM PPO =1. /pp PCR=-2.-VSX*CX* (SHY**2+CHYIN*2) PCI =2. *SHY —CHY*(CXX'1*2-SX**2) ACC= A L." CC, ASS= AL SS ACS=AL-."CS ASC=Al-:cSC F=l.4142i4*, (A""CC+PB*SS) TANH (ASS) (BCC-A""SS,)"TAN ~'(ACC) -PP+PM cTAN ( ACC),TjA'NH (ASS) *(ASS))-PM-PP*TAN(ACC))#;cTANHF(4SS) C ~.*PREASSIGN VALUES AS THE TOLENCE FOR THE CONVERGENCE IF(ABS.(F).LT* 0.001.AND. ABS(G) %LT. 0.001) Go *TO 11. FPR =1.*4142 144c (TANIH(ASS)*;'(-A*SC+B*CS ) +(A*CC '+MH*.SS)"( A CS)/((COSH(ASS A —

74 * ) ) 2) + ( AS C ):-c ( R.*CCC-A*SS ) (( (S ( ACC ) ) *`2 ) + #'TAN( ACC)* SCA*(, S )-O.5 c*-:PPO-*-( P,R+ O*(C* ---P(R-I)*-PC I ))+0.5,cPM 0* (-PCR+O *,-TAN\H( ASS )+PMi4*( T ANH (ASS ) ",(-ASC)/ (CUS (ACCH* *2+( ACS )*T AN( ACC)/ (JIS H (ASS)) )"*42 FP 1=1. 4142 14*I' (TAN( (ASS) *( A" CS +11*SC )+ (A *CC+B **2 ) C( CS) (*C C,- A,'S-S )/(GS (CO A C.C) )`2 -T AN (AC CC)4( - AC'.'SC) )-O.&5 *p pn*(pCI +O*, C P CI +D*,-,PC R)+ 0 5 *PMO 4( P CI +(*(:tC):-.P C1.+O P C R) TIA N(AC C*T AN H (A S S+ P M*(AC S ) "T ANH (A S S)/COS( *ACC) ) ",2+(ASC.) "-TAN( ACC/ ((COSH(ASS ) )~2) IPR = I. 4142 140c( (ASC)*(A*,,CC+R,:cSS)/((COS (ACC *))~~'NACC ) (-A fsc *+*C-S) +( ACS)* CCASS)/(COSH (ASS))-`2)+TA *NH ( AS S) -`,(-B4rSC-A"CS *)0. 5*,PMO0*(-PCR+O ( C*"PCR-D,',PC In -o * 5pPP0* *(PCR+O':I- C-*PCR-D-* PCI) *T AN(~ ACC C, *TA NH( A SS )P P( (A SC )T A NH ( A SS C/( (CS(A GPI=1.414214*((ACS)"'CA'*CC+B"'SS)/((COS(ACC) *)42) +TA N (AC C) (A C S *-+BF-:*SC) + (ASC) (p~,*CC-A,*SS) ((COSH (ASS.))4r2)+ *TANH ( ASS)*(`,*CS-A):~SC) ',+O. 5-*PfA0*(PC I -O*'(C"'PC I+D"PCR ))-O,5*PPO'(PC *I+0l (',C,:PC I +D0*PCR ) ),-*TAN *l(ACC)*:TANH(ASS)-PP*((ACS)*TANH(A4SS)/((COS%(ACC) )42) +(ASC)*TAN (ACC)/(COSH(ASS))**2) l= F PR *GPI-FPlI*GPR D)X= ( *FP I -F ':cGP I ) /DD F)y-( F*'CPR-G*,FPR )/DD WR I TE (6,r10) X,Y, F 9G X=X+DX Y=Y+D)Y GOn TO I 11 WRITE (6,10)XYFG LEO)=K`(SC+(O. l,1)*CS) HL=ACC+(0.,-. ) *ASS WlRI TF(6912 )LED WJRI TE(6112)HL GO To 20 10.F0RMAT(2XrE11.4,p2XgEl1.42XYE10.3,p2XjElO.3) 1000 FOR M A T (2 F1,5) 12 FORMAT(2X1,2E12.5) STOP END