Report 390573-2-T J.L. Volakis and J.M. Jin On The Formulation and Solution of the E-Field Integral Equation The University of Michigan Department of Electrical Engineering and Computer Science 1301 Beal Avenue Ann Arbor, MI 48109-2122 June 1991 Northrop Corporation 8900 E. Washington Blvd. Pico Rivera, CA 90660-3737

J unto e: p^ i " I

1 Introduction In this report the discretization and solution of the E-field integral equation is considered for scattering by open and closed resistive boundaries. Two specific forms of this integral equation are studied and their solutions are compared. One form of the integral equation involves the usual unknown surface current density and the other is in terms of the charge distribution. Provided the charge and current expansions are chosen to satisfy the continuity equation, it is shown that the two integral equations result in identical systems. This is analytically verified and to do so it was necessary to employ linear weighting functions for both integral equations. The conclusion of this study is that by resorting to higher order basis for expansion and weighting, one can readily reduce the singularity of the kernel without a need to reformulate the integral equation. 2 Derivation of the Integral Equations Let us consider the two dimensional surface shown in figure 1. This surface may be open and it is therefore a curved strip whose unit tangent and normal shall be denoted by s' and n, respectively. The strip has a normalized resistivity R and is illuminated by an H-polarized (TE) plane wave Hi -= eJkO(x cOs o+ysin)o) E' _ jZ = =Z0 ( sin q, -'cos O,) ejko(xcosso+y sinko,) ko Oy ax (1) where ko denotes the free space wavenumber and Zo is the free space intrinsic impedance. To derive an integral equation for the current supported by the resistive strip, we recall the boundary condition E8= ZoRJS (2)

where E, = s.-E is the s component of the total electric field at the surface of the resistive strip, J, denotes the surface current density along the surface at (x, y) and R is the normalized resistivity of the strip also at (x, y). The total electric field is given by the sum of the incident and scattered field, i.e. E = E'+ E' (3) and since E' is caused by the strip current we have Es = - k0Zo s J, (p') Ho2) (kor) ds' -' 4kL SJ' (p') VVH(2) (kor) ds' (4) In this H(2)(.) is the zeroth order Hankel function of the second kind and r =Ip - = - (') + (- yl)2 is the distance between the observation and integration points. Also, by refering to figure 2 we observe that s^' is the unit tangent vector to the strip at the integration point (x', y'). From (4) we can now extract the s component of E' to be used in (2). We have EJ= SEs*ES= - 4 k4Z J T (p') S{ S'HO2) (kor)} ds' - 4ko J J(p {(S V)(S^ V)Ho(2) (kr)}ds' (5) and upon noting that s.' V = -'. V' = -V,, we can rewrite (5) as - k-Z4 {| J, (') S^ s'H(2) (kr)] ds' 1 r 2~~~~~~~~~~6

Substituting this into (3) and then into (2) we obtain the integral equation YoE'(s) = RJ(s) + 4 f J:,(s') {(s. S') Ho2) (kor)} ds' 4 ko J1( s' s (2 (7) where YO = 1/ZO and for convenience we have replaced J,,(p') by Ji,(s')in which s' denotes the distance along the strip. This is the standard E-field integral equation [1,2]. By applying integration by parts we can obtain other integral equations which are equivalent to (7) and more suitable for numerical implementation. Noting that J8(s) vanishes at the strip ends upon applying integration by parts to the second integral of (7) we obtain Ic I 2) YEs,(s) = RJ8(s) + - J8,(s'){( S') H2(kor)}ds' + L-j dJ,(s') HO2)(kor) ds' (8) To do the same for the first integral we first define the function Gi(kor)= j (s. ") 2)(kr') ds (9) where a is an arbitrary constant and r' /(x -x"2 + (y -y2 with the integration being with respect to the double primed variables. Then Gl(kor) satisfies the identity dG(kor) = (.') Ho2) (kor') (10) ds te permitting us to rewrite (8) as

YoE(s) = RJs(s)- - (s')Gj(kor)ds' + O a H'(2' ) 8 (ko') )ds' (11) in which O(S) da(s) (12) ds is a quantity proportional to the electric charge on the strip. It can also be shown (see Appendix) that J(s) = j (L )26 [(s + 6) - (s - )] d (13a) where L denotes the length of the curved strip and q(L - s) =k(s) (13b) should be used in evaluating the integral. Substituting this into (11) yields the integral equation 621 (r)Jcd ( L [(+- -2s')] YoE'(- - CR(s) ds' (S + S') - q(s - ko k- | k (s')Gj(kor)ds' + H(2)(kor)ds (14) + 4ko Ic ('o, (1 Clearly, (14) is solely in terms of the charge distribution on the strip whereas (8) is the tranditional E-field integral equation involves only the current density. Consequently we shall refer to (8) as the (E-field) current integral equation whereas (14) shall be referred to as the charge integral equation. A standard procedure for solving (8) or (14) is the method of weighted residuals. In applying this technique the integral equation is multiplied

by a testing functions wi(s) and then integrated over the domain of the weighting function. Doing so, (8) and (14) become, YoJ E(s)wj(s)ds J R(s)J.(s)wi(s)ds + 4.. + 4 Jc; wi(s) Jc Ji(S') {(S * sI)H(2)(k r)} ds'ds - E, 1, dw;(s)ds L3- H~dJ8~s')2)(kor)ds'ds (15) 4ko i J )ws)ds - 1 ds' Y E()w2(s)ds =- R(s)w(s) ( L - 2s') [4(s + s') - (s - s')] ds'ds - JC Wi(S) JC (s')G(kr)ds'ds 4ko ~ ds (s') 2)(kor)ds' (16) for i = 1,2,...,N such that Z=1 Ci = C. That is, Ci represents the i th segment of the discretized strip. Also, in obtaining the last integral in (15) and (16) we employed integration by parts. Thus, in their present form both integral equations are associated with kernels which have an integrable logarithmic singularity. In contrast the kernel for the original integral equation (6) has a non-integrable singularity which precludes its implementation in any rigorous mathematical fashion. Nevertheless, it has been implemented [1.] using pulse basis expansion function with reasonably good success. The implementation of (15) and (16) requires that we first expand the current/charge as N J((s)= JjLj(s) (17) j=1 N j()= ZlXjj(S) (18) j=l

where Li(s) and /kj(s) are subdomain basis functions to be chosen. That is they are non-zero over the segment Cj and vanish elsewhere on the surface of the strip. Substituting these expansions into (15) and (16) yields YO E,(s)ws (s)ds =, j R(s)Lj(s)wi(s)ds j=l 1 N dwi(s) dLj(s') H(2)((kor)ds'ds (19) 4k0., ic ds ds' YOf E.(s)wj(s)ds = -2 R(s)wj(s)j ( -.s) [O(s + s') - (s - s')] ds'ds 2 j=l i ( 5 -i,0 O d Wi(,) _ j(s')G)(kor)d('d1 - j=l i 1 () dwi(s) (21)(SI) HO(2) (korrdssds (20) o j=1 i tand it remains to specify the testing and expansion basis functions for a solution of the urrent/charge. Since (19(s) is equal to the derivativegral of the current J,(s), it is logical to choose'j(s) so that Oj (S) = Lj(s) (21) ds That is, if Lj(s) is chosen to be the triangle function Aj(s) as shown in figure 3, then in accordance with (21), Oj(s) becomes the doublet function Dj(s) also shown in figure 3. With this choice of pj(s), it is seen that Jj is equal to Oj and if wi(s) are chosen to be the same for both integral equations, it follows that the last integral of (19) is identical to the last integral of (20). It can also be shown that the second integral of (19) can be readily written (via integration by parts) in a form that matches the corresponding integral in (20). This simply proves that by using appropriate choices for

the expansion and weighting functions, the resulting systems from (19) and (20) are identical. Of course, in principle, one could choose the expansion and weighting functions for (19) and (20) to be different. However, in view of the charge conservation requirements, the above choice for,bj(s) is the most appropriate. Expansion basis other than the triangle function Aj(s) were also considered. However, we have found that if pulse basis are chosen for expanding the charge along with point-matching, the resulting system is unstable. After going through various implementations of (20) with different combinations of testing and expansion functions, it was found that the best lowest order choice for the weighting function is Aj(s) and the same holds for the expansion functions Lj(s). In accordance with (21) this also implies that an appropriate choice for bjj(s) is the doublet function Dj(s). Next we consider the evaluation of the matrix elements associated with current and charge integral equations. First we compute the matrix elements for a flat perfectly conducting strip. 3 Solution of the current integral equation for the flat strip Consider the perfectly conducting flat strip shown in figure 4. Choosing the weighting and expansion functions as z(il)x z (j - 1)Ax < x < jAx wj(x) = Lj(x) = Aj(x) = (j+l)Ax-x jAx < x < (j + 1)Ax from (19) we obtain the system [Aij] {Jj} = {bi} (22) In this, [Aij] is a square N x N matrix whereas {bi} and {Jj} denote column matrices for the excitation and unknown current density. Directly from (19), the elements of the excitation column are f(i+l)A bi = Yoi i Ai(x)Ei(x)ds (23)

YEx(iAx)Ax and those for the impedance matrix [Aij] can be written as A3= 4-(aii-2aij) (24) in which J (,_,)AX (j+1)Ax ~aij =()A itA(x) j)A (j(xz)Ho2) (ko Ix - x'l) dx'dx (25) a iji)Ax dAi(x) (j+)A d(z)H2) (ko I - x'l) dx'd (26) i (i-1)a dx (j-1)ax dx' To evaluate aij we may employ the 5-point integration formula. We have 2 2 aij- a cn E a,,minjm (27) n=-2 m=-2 where It, =i+2n-)Ax (j 2 1) Ho(ko Ix x') d'dx (28) and n -2, 2 an = n =-1,1 1 n-O The integral Ii'5m can be readily evaluated numerically using the mid-point integration formula provided i: j or n ~ m. When i = j and n = m the integrand has a logarithmic singularity and we must then evaluate i" analytically. In this case ko Ix - x'I < 2k0o- << 1 and we can thus employ the small argument approximation for the Hankel function, H12)(z) = 1 - j-ln (29) 2

Substituting this into (29) we obtain zInn(=) L 2(i 2ii k0(x x') 6 6( - L - In -) dx'dx +62 _ J 2 ~ ( in- ) ddxd = 62 j{j _ )iak(x -,)dx' + j k In dxd = 62- j { _ [(x)l ~( )] XI lfl kox_-_x')P2 r ~= 62 _ f(z + yk 2 () +2 e 2 - x) 62- j 2 - (x - )n yk, x + x) 7 716 2e +( - x)ln7 ( jdx = 62 j yln7 jdy (30) 4 (2e I~ (6) =52-_j42 2 jezlnzdz ~' yyo j yn tky y 90

- 2 - ( z -2r -Inr -lZ r - 7ko) 2 4]0 2 2 2e 2 y(7k6\2 yk26 ir J7jo 2e 2e = [1 j2In (7k )2 (31) We now turn our attention to the evaluation of aUj. Upon substituting the expressions for Aij(x) and differentiating we obtain t iAx 1 (i+l)Ax 1 i jar 1 (+1 1 aij LJ(i- )a Az Jix x J(jl)a Jj a = [i-lx,&x f~x Ax iJfj-1)Ax Ax 7Ax Ax] 2) (ko I - x'j) dx'dx} _ 1 f rfitx f(i+l)AX] [f3AX -(f+1)Ax (ax)2 [tJ(i-1)Ax Jiax J -)jaAx JAx * H12) (ko x - x'I) dx'dx} (/1 f)2 {o Ax [2Ho2) (ko t(i -j)/Ax + x - xJ) - Ho2) (ko I(i - j - 1)Ax + x - x'1) - H(2) (ko (i -j + 1)Ax + x - x'I)] dxdx' (33) which can be evaluated numerically provided the integrands are not singular there. This gives a' - 2Ho(2) (ko li -jl Ax) - Ho(2) (ko li -j -11 Ax) -Ho2) (ko i - j + 11 x) (34) 10

When the integrands are singular, which occurs for the first term when i - j, for the second term when i = j + 1, and for the third term when i - j - 1, we may use the result of (31). In particular, we have j j H12) (ko Ix - x') ddZ = [ jln 2e (X)2 (35) 4 Solution of Charge Integral Equation For the solution of the charge integral equation we again choose the same weighting functions defined in (23) but in view of the requirement for charge conservation, the expansion functions are chosen to be 1-j x (j - 1)A x < x < jx Oj (x) = Dj (x) = 1 jAx < x < (j + )x (36) Substituting these into (20) with R = 0 we obtain the system [Aij] {)} = {bi} (37) where bi is again given by (24). The elements of the matrix [Aij] are defined by Aij iij (38) 4 k2 with ~= - I Ai(x) ] Dj(x')H22) (ko Ix - x"1) dx"dx'dx J(i —1)Ax J (j-1)Ax = —i( h,\ Ai(x) - dhj(' )(ko Ix - xtl) dx'dx i'1)A A [(. j-il)Ar Aj(x)H o) (ko Ix -'1) dx'd = aj (39) 11

and,: (,+,A d,(x) /(j+.')Ax a = (i1)A dA (X Dj((x')H(2) (ko Ix - x'j) dx'dx (i-1)Ax dx -(j-1)Ax (i+l)Ax dAi(x) j(+1)Ax dAj(x)(2) (k x - x'I) dx'dx (i-1)Ax dx (j-1)Ax dx' - ai (40) That is, the elements of the matrix [Aij] are identical to those of the matrix [Aij] associated with the current integral equation. This was, of course, to be expected since as noted in the previous section Oj = Jj. 12

5 Appendix: The strip current in terms of charges Mitzner [3] derived that for a closed two-dimensional surface whose circumferential length is L (see figure Al), the surface current density satisfies the relation Ks(s) = Ks - ( L) [g(s + 6) - g(s - 6)] d6 (Al) where Ks =K I(s)ds = average current on C (A2) ds g(s): da, (A3) and g(s)ds = K(L)- (O) = 0 (A4) since K,(s) = K,(s + L), i.e. K, is a periodic function having a period L. We like to derive an expression similar to (Al) for the current on a curved strip. Equation (Al) is still applicable for the strip except that C is now made-up of the top and bottom surfaces of the strip (see figure Al). The net current on the strip is given by J.a=K - K K- = K.(s)- Ks(L- s) (A5) and from (Al) on letting s L - s we get,(L ~- s) = K8-, - L 26) [g(L - +6) 2 L -g(L - s - 6)] d (A6) Combining (Al), (A5) and (A6) yields J = K.(s)- K:(L- s) 2 L [g(s + ) - g(s - 6) - g(L - s + ) + g(L - s - )] d (A7) 13

However, since dJ, d [K,(s) - K,(L - s)] dK,(s) dK,(L - s) ds ds ds ds dK8(s) dKI 8(L - s) ds d(L - s) it follows that dJ8 s= g(s)+ g(L - s) = q(s) (A8) ds where q>(s) is proportional to the net charge on the surface of the strip. Using now (A8) into (A7) gives J=~j (L L) [q(s + )-q(s-S)] d (A9) and it should be noted that,(s) = g(s) + g(L - s) = g(L - s) + g(L - L + s) = q(L - s) (A10) Expression (A9) was implemented and found to hold when J, = 0 at the ends of the strip. An example calculation of the relation between the current and charge on a 3A flat perfectly conducting strip is given in figure A2. Also figure A3 shows the same quantities for a 1A flat strip. It is clear from both of these figures that the charge distribution has a much larger dynamic range than the current density. This probably implies that the computation of the charge density is more difficult. To test the validity of (A9) we can substitute the charge density shown in figures A2 or A3 and attempt to recover the original current density distribution. The results of our first attempt is shown in figure A4 and as seen the re-calculated current, albeit close, is not precisely in agreement with the original current density (see routine in figure A6 used to implement (A9)). As it turned out, this is because the original current density was not exactly zero at the ends of the strip (because the MoM code does not sample at the strip ends) which violated the conditions for which (A9) holds. Since the re-calculated current density does vanish at the ends of the strip, this can then be used to validate (A9). Differentiating the current given in figure A4 and substituting the result in (A9) it was found that the original "modified" current density is recovered exactly as demonstrated in figure A5.

References [1] E.F. Knott and T.B.A. Senior, "Non-Specular Radar Cross Section Study," University of Michigan Radiation Laboratory Report 117641-T (also U.S. Air Force Report AFAL-TR-73-422), Jan.1974. [2] M.A. Ricoy and J.L. Volakis, "Integral Equations with Reduced Unknowns for the Simulation of Two Dimensional Composite Structures," University of Michigan technical Report 389492-2-T, Nov. 1987. [3] K. Mitzner, personal communication.

Js = Ks+ - Ks Ks+ O~~~~~~~ Ks' Figure 1. Geometry of the 2D curved surface (strip). n _A~~~~ a~~~~~~ Ix, Y) S~ ~~~~~S / ~~r 0 Figure 2. Illustration of the observation and integration point parameters.

A j(s) J I i Sj- A Sj+A S (a) D j(s) J S j- A Sj.Sj+ A -1 (b) Figure 3. Expansion functions. (a) Linear expansion functions for the current (b) corresponding expansion functions for the charge.

SSi C (a) (b) Figure Al. Geometry of a closed surface. (a) Original geometry of circumference L (b) Currents on a closed surface collapsed to a strip.

CURRENT AND CHARGE FOR A 3 WAVE. STRIP ANGLE OF INCIDENCE = 45 DEG OFF GRAZING 20 RE (J) M(J)............ RE GE) )......... IM(CHARGE) 10 -2.,- 0 0o 1 2 3 Figure A2. Current and charge densities on a 3X flat strip. Routine to compute charge (CHARC) from current (PHI) DO 66 I=1,M IF (I. EQ. 1) THEN ELSE IF (I.EQ.M) THEN CHARG (I)= (3. *PHI (I) -4. *PHI (I-1) +PHI (I-2) ) / (2. *DSQQ) ELSE IF((I.NE.1).AND.(I.NE.M)) THEN CHARG (I) = (PHI (I+1)-PHI (I-1) ) / (2.*DSQQ) ENDIF 66 CONTINUE

Mag. (J) Phase (J) Mag. (Charge) Phase (Charge)..... 8 200 7 *"...... 150........................k................................... t.............................................. 1 0 0 6 l,,. -\'_}...................... CZg 6''1 0 6 _t o-.,. -' " u.'' I i' 0r 02 04 060 1.6............'.........................................I......................................................... I. - 50 1 00a 2' - 150 -1~~~~, ~-200 0 0.2 0.4 0.6 0.8 1 x Figure A3. Current and charge densities on a 1, flat strip.

Solid lines ~ Current Density (original from MoM code) Dashed: Recalculated using (A9) 3 - 1 1 I, 200.";'"-..s _. ~.~, - 150..I......... I.............. Z........................... 00 X A ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~I. I O =. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~.................'.'.................... I......... I... c, ~ ~~~:,_ J: 5 =~ 1.5................................ -.. -................................................'............... ~ ~ ~ ~ ~~~~~~~~~.: i...~......~.~.. I~..~...~. iJ'~t.................................................. 0~~~~~'2 " L\ 0 0.5 1 1.5 2 2.5 3: Figure A4. Recalculation of Js from using the formula ('). As seen, Js is not precisely zero at the ends.

ABS(J) - -PHAS(J) - - -ABS(J)-RECALC ---— PHAS(J)-RECALC 3 1 80 2.5 ------.................................................y 21~~~~~~~~~.5 i- A'X + 9 1.5, —................-.... —0.5 -1.80.............. 0 0.5 1 1.5 2 2.5 3 x Figure A5. Recalculation of Js from Q using equation (A9). In this case the original Js is given by the dashed curve in figure A4.

Routine which Implements Equation (A9) Js = - ds (L) [+) s-8)] U2 = M * DSQQ PHI = J s C REcover current density CHARG=e DO 68 J=1,M PHJ(J)=(O.,O.) DO 67 I-1,M RII=I-1. +. 5 RJJ=J-1. +. 5 RI 1=RI I+RJJ RI2=RJJ-RI I I1=I+J-1 I2=I-J IF (I1.LT. 0) I1=2*M+I1 i 1 i+2 IF (I2. LT. 0) I2=2*M+I2 IF (I1. GT.M) Ii=2*M-I1 IF (I2.GT.M)I2=2*M-I2 ( zi+1 ) IF(I1.EQ.0) THEN CHP=CHARG ( 1 ) -0.25* ( 4. *CHARG ( 2 ) -3. *CHARG ( 1 ) -CHARG ( 3 ) ) ELSE IF (I1.EQ.M) THEN CHP=CHARG (M) +0. 25* (3. *CHARG (M) -4. *CHARG (M-1) +CHARG (M-2) ELSE CHP= (CHARG (I1) +CHARG (I1+1)) /2. ENDIF IF(I2.EQ.0) THEN CHM=CHARG (1) -0. 25 (4. *CHARG (2) -3. *CHARG (1) -CHARG (3)) ELSE IF(I2.EQ.M) THEN CHM=CHARG (M) +0. 25* (3. *CHARG (M)-4. *CHARG (M-1) +CHARG (M-2) ) ELSE CHM= (CHARG (I2) +CHARG(I2+1) )/2. ENDIF PHJ(J) = ((2. *M-2 *RII) * (CHP-CHM) / (2. *M) ) +PHJ (J) 67 CONTINUE PHJ (J) =-.5*DSQQ*PHJ (J) PHIA=CABS (PHI (J)) PHASI=dig*atan3 (aimag (PHI (J) ),real (PHI (J) ) PHJA=CABS (charg (J)) PHASJ=dig*atan3 (aimag (charg (J)), real(charg (J))) WRITE (8, 950) PHJ (J) WRITE (9,952) (J-0.5) *DSQQ, TAB, PHIA, TAB, PHASI, TAB, PHJA, TAB &, PHASJ 68 CONTINUE Figure A6. Routine to recover the current from charge.

UNIVERSITY OF MICHIGAN 3i1 9111 011121111111 111111111 3 9015 02229 1580