U —M Technical Report: RL-864 Theore tical Development and Numerical Implementation of a Domain-Boundary Integral Equation for TwoDimensional Electromagnetic Scattering Jian-Ming Jin and Valdis V. Liepa Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109 April 6, 1989 RL-864 = RL-864

SUMMARY This report describes the theoretical development and numerical implementation of a domain-boundary integral equation (DBIE') for two-dimensional electromagnetic scattering. More specifically, it contains the following: * Derivation of the DBIE using two different approaches; * Interpretation of the DBIE from a physical point of view; * Moment method solution of the DBIE using pulse basis expansion and pointmatching technique (code included); * Moment method solution of the DBIE using isoparametric elements and point-matching technique (code included). Of particular importance are the two distinct properties of the DBIE: it is based on only one field component and does not involve derivatives of the field. The first property allows the use of a minimum number of discretized unknowns. The second allows an easy application of pulse basis expansion functions. The further development of isoparametric elements for the moment method results in a numerical solution of high accuracy. i

PREFACE This report consists of three papers and two numerical codes which represent our research work conducted during October 1987-October 1988 on the domainboundary integral equation. In order to maintain consistency with the frequentlyused phrases "Volume Integral Equation (VIE)" and "Surface Integral Equation (SIE)", we had originally named our new equation the "Volume-Surface Integral Equation (VSIE)". Though the phrases VIE and SIE are often used for the twodimensional case, strictly speaking they belong to the three-dimensional case. On the other hand, the phrase "Domain-Boundary Integral Equation (DBIE)" is suited for both two- and three-dimensional cases. Since this report primarily deals with two-dimensional scattering, we decided to use DBIE for its title. It is, therefore, understood that the DBIE used in the title and summary is equivalent to the VSIE in the text. The authors wish to acknowledge the discussions and suggestions of Professors C. T. Tai and J. L. Volakis in relation to this work. ii

TABLE OF CONTENTS Summary.............................................. i Preface...........................................................ii A Volume-Surface Integral Equation for Electromagnetic Scattering by Inhomogeneous Cylinders by J. M. Jin, V. V. Liepa, and C. T. Tai....... 1 A Simple Moment Method Program for Computing Scattering from Complex Cylindrical Obstacles by J. M. Jin and V. V. Liepa................. 25 A Moment Method Solution of a Volume-Surface Integral Equation Using Isoparametric Elements and Point-Matching by J. M. Jin, J. L. Volakis, and V. V. Liepa................................................... 58 Description and List of Code VSIEM by J. M. Jin.................... 78 Description and List of Code VSIE-ISO by J. M. Jin................ 110 111

A VOLUME-SURFACE INTEGRAL EQUATION FOR ELECTROMAGNETIC SCATTERING BY INHOMOGENEOUS CYLINDERS Jian-Ming Jin, Valdis V. Liepa, and Chen-To Tai Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, Michigan 48109 (published in The Journal of Electromagnetic Waves and Applications, vol. 2, no. 5/6) ABSTRACT This paper presents an integral equation formulation for electromagnetic scattering by inhomogeneous infinite cylinders having arbitrary scalar permittivities and permeabilities. The formulation involves both volume and surface integrals with only one unknown field component, and is applicable to both transverse electric and transverse magnetic cases. This volume-surface integral equation is wellsuited for numerical implementation. In this paper, the integral equation is first derived by integrating the wave equation with the aid of the free-space Green's function, and then analyzed from the physical point of view, resulting in a new interpretation for the scattering mechanism. Finally, the equation is programmed using the method of moments with pulse expansion functions and point-matching. 1

Numerical results are shown to demonstrate the validity of the formulation and the new interpretation of the scattering mechanism. 1. Introduction The problem of electromagnetic scattering by inhomogeneous cylinders can be formulated in terms of equivalent polarized currents (see, e.g., [1]-[4]). The formulation leads to a volume integral equation which involves three unknown current components for a scatterer having both permittivity e and permeability fu different from their free-space values. If the cross section of a cylinder is divided into N cells for numerical analysis, one has 3N unknowns; consequently a matrix equation of size 3N x 3N needs to be solved. However, the three unknown current components are not independent; they are related by Maxwell's equations. Using such relations, one should be able to reduce the number of unknowns. Recently, a compact formulation has been developed by modifying the equivalent currents, reducing the three current components to two and thus resulting in 2N unknowns for numerical analysis [5]. In this paper, the problem is formulated directly using the field concept, instead of the equivalent polarized currents, by integrating the wave equation with the aid of the free-space Green's function. The resultant equation contains both volume and surface integrals but involves only one unknown field component. This volume-surface integral equation provides an alternative physical interpretation for electromagnetic scattering, and is well-suited for numerical implementation. Using this equation, if the circumference of the cylinder is divided into Al segments, we 2

have only N+M unknowns and thus the resultant matrix size is (N+M) x (N+M). It is noted that a similar idea has been employed by Tai to formulate the scattering by a homogeneous electromagnetically permeable body and by an inhomogeneous dielectric body for the general three-dimensional case [6]. 2. Derivation of the Integral Equation Consider transverse electric (TE) and transverse magnetic (TM) wave scattering by an inhomogeneous cylinder having its infinite dimension along the zdirection (see Figure 1). The complex relative permittivity Cr(r) and permeability /t,(r) are continuous functions of position, where r is the position vector in the xy-plane. Let Q denote the region inside the cylinder, Qo the region outside the cylinder, and r the boundary of the cylinder separating Q and QO,. In the region Q, the z-directed electric or magnetic field, denoted by F', satisfies the wave equation (i o, P lo Source ~~ v E(r),| (r) Scatterer Figure 1: Geometry of wave scattering by an inhomogeneous cylinder. 3

V [u(r)VF(r)] + kov(r)F(r) = 0 (1) where for TM incidence 1 F(r) = E(r), u(r) (= -,(1) Tr(r) and for TE incidence 1 F(r) = H(r) u(r) =, v(r) -= (r) In the above, V is a two-dimensional operator defined in the rectangular coordinate system by V = x(a/0x) + y(a/ay) and ko is the free-space wavenumber. To find an integral equation for F, we integrate (1) over the region Q with the aid of the two-dimensional free-space Green's function Go(rli'), which satisfies the equation V2Go0(rr') + koGo(r') = -8( - r-) (2) This leads to the equation J J {Go(rr')V. [u(r)VF(r)] + -k2v(r)Go(rr')F(r)} ds = 0 (3) Then, using the identity V (GouVF) = GoV (uVF) + uVF * VGo and the two-dimensional divergence theorem, (3) can be written as J j [kl2V(r)Go( 1r')F(r) - u(T)VF(i) VGO(r-lr')] ds + Go( I')u(r) (9F(:)dl = 0 (4) J OdT 4

where ni is the outward unit vector normal to the cylinder's surface. Substituting the identity V (uFVGo) = uFV2Go + uVF VGo + Vu (FVGo) into (4) and applying again the two-dimensional divergence theorem, we obtain J {v(r)Go(rlrT')F(r) + (r)F(r)V2Go(rI r') -F Vu(r). [F()VGo(r I')]} ds + f G - ru( r)F^)aGll] di- 0 (5) Substitution of (2) into the above equation gives ko f j[v) -- u(r)]F(r)Go(r-r')ds + j Vu(r) [F(r )VGo(rr')]ds + 1 [()Go(r') (-()F(r)FG() r)a dl u(r')F(r') for r' e -= u(r')F(i') for r' on r (6) 10 for r' E QO Note that the line integral is understood to be a principal value integral when r' is on F, i.e., the singular point r = r' is excluded from the integration. Interchanging the unprimed and primed coordinates and using the symmetry property of the Green's function Go(rir'), the above equation can be written as k2 j[v(r') - (r')]F(r') Go(r I r')ds' + J Vu(r-'). [F(r')V'Go(rij?')]ds' + [u(')Go(rIr') a- ' -u(r-')F(r') a('")] dl' I u(i)F(i) for r 6E Q = u (u(r)F(r) for r on F (7) 0 f for r- E O 5

This is the equation resulting from the integration over the interior region Q. Now let us consider the exterior region QOO. In this infinite region the field is governed by the wave equation V2F(r) + k2F(r) = jwS(r) (8) and the radiation condition at infinity. The function S(r) is related to the zdirected electric or magnetic current by S(r) = 0oJz(r-) for E-polarization and S(r) = eoAfz(r) for H-polarization. Let us assume that the source distribution is confined within a region Qf. Integrating (8) over Qfc with the function Go(rlr'), we obtain F:JN(}) - fJ [G0o(r') F')( )-F(.') DGo(If ')] dl' 0 for r CQ = F(r) for r on r (9).F(r) for rE C Q0 where FINC is the incident field given by F'I C() = -jw J j S(r')Go(r i')ds' Equation (9) is the well-known surface integral equation. Since the boundary conditions require the field F and the quantity u(9F/&an) to be continuous across r, we can combine (7) and (9) to obtain an integral equation which does not involve the normal derivative of the field. The resultant equation is then FIC(rT) + ' / j/ [v(r -u(r/)]F(f)Go(/)ds 6

+ V() [F(')V'G (ri)]ds' + [1 - u(r')] F(r ') dl' u(r)F(r'i for r E Q - [1 + u(r)]F(r) for r on r (10) F(r) for r E Qoo and we call it the Vlolume-Surface Integral Equation (VSIE). It should be noted that the above methodl is not the only way to derive this equation; it can be also derived using the method of equivalent polarized currents. Such a derivation is given in Appendix I to show the difference between the two methods. n" Figure 2: A cylinder consisting of two inhomogeneous media. The above integral equation (10) is valid for problems with u(r) continuous in Q. For cylinders consisting of two or more continuous inhomogeneous media, but having step discontinuities in u(r) at their interfaces denoted by rd (see Figure 2), a line integral / [u ) - I')] F(r') ~ r dl (11) (11 7

where id points from the - side to the + side, must be included in the left-hand side of equation (10). n Conductor Figure 3: A conducting cylinder coated with inhomogeneous material. To deal with scatterers containing perfect conductors, such as coated cylinders (see Figure 3), a line integral eu(r')G(rr-) F(r') - SU(r)F(- G ') dl' (12) along the conducting surface Pa must be included, where nia points toward the inside of the conductor, and 6h and 6S are polarization factors defined as { Sh = 0 for E-polarization; 1 for H-polarization e = 1 for E-polarization; 0 for H-polarization 3. Physical Interpretation In this section, we use the concepts of electric currents and dipoles to gain insight into the individual terms of the volume-surface integral equation and interpret it from the physical point of view. For this, let us rewrite (10) for the TM 8

case as follows 0 j[Ir(i) -]Ez( )Go(rr')ds' + k ) 1 Erz(')Go(rlr)ds' + J vj ( 7] [Ez(r')V'Go(r r')]ds'+ j 1 EZ(r')V'Go(rl r ) ndlt [1//r()] Ez(r) for r C +ENTC(r) -= - [1 + 1/jt(r)]Ez(r) for r on F (13) E,(r) for rG C o Recalling that a z-directed two-dimensional electric current, denoted by Jz, produces a field Ez() = -jo / Jz (r')Go(r r')ds' (14) and a volume distribution of electric dipole moment, denoted by Pv, produces a field [7] Ez() =jW oJJ Pv(r') VGo(lr')ds' (15) and also that a surface distribution of electric dipole moment, denoted by ps, produces a field Ez(r) = jwo jps(r') VGo(rr')dl', (16) by comparing these with (13) we see that the scattered field results from the superposition of the fields produced by four sources. The first source term is an equivalent electric current J(1) jO(Cr - 1)Ez (17) 9

the second is another equivalent electric current jz(2) =j - 1 E, (18) /l r the third is an equivalent volume electric dipole PE=- VE (1) (19) and the fourth is an equivalent surface electric dipole p, = -Z-1 E (20) The second electric current J(2) differs from the first one J(1) in that it produces an additional term at the source; the Green's function, denoted by G(rlr'), applied to J(2) is G(r') = - + G- (r) (21) k(22 rather than Go(rlr') that is applied to J(1). From (17) - (20), we see clearly how each source is introduced. fie believe that the above interpretation offers a more physically realistic picture for wave scattering than the picture that follows from the conventional volume integral equations which interpret the scattered field as the superposition of the fields produced by equivalent electric and magnetic currents. To interpret equation (10) for the TE case, it is convenient to use magnetic current and magnetic dipole concepts. The interpretation is similar to that for the TMI case; the scattered field results from the superposition of the fields produced 10

by magnetic current sources and magnetic dipole moments of volume and surface distributions. Such an interpretation, however, falls back to the nonexistent magnetic sources again as do the conventional volume integral equations. 4. Numerical Implementation and Results This section deals with the programming of equation (10), which is solved by the method of moments in a straightforward manner. Since the equation does not involve any derivative operation on the unknown (as usual, for scatterers containing perfect conductors with TM incidence the quantity aF/an in (12) is treated as an unknown rather than the normal derivative of an unknown), the pulse expansion and point-matching technique [1]-[3] can be easily applied to give accurate results. This is in contrast to the conventional volume integral equations, which contain derivatives operating on unknowns for the general case, and thus pulse functions are inappropriate for use in expanding the unknowns [3]. Hence, the volume integral equation approach requires more effort in selecting expansion functions to obtain accurate results [8]. Since the method of moments is well documented, here we only give the result of discretization. We also only consider the problem illustrated in Figure 1; the problems of Figures 2) and 3 can be treated in a similar manner. Assume the region Q is divided into N cells and the boundary F is divided into M1 segments. If we use Xi to denote the discretized fields in Q and /b the fields on F, the application 11

of the pulse expansion and point-matching technique gives the matrix equation -N (22) q5INC where A, B, C, and D are submatrices with dimension of N x N, N x M, AM x N, and M x M, respectively. The expressions for their elements are given in Appendix II. On the right-hand side of (22) is the discretized incident field. The solution of the problem is obtained by solving (22) for X- and 5b. Once these discretized fields are found, the far field is calculated using (10) with the large-argument approximation for G0(rQr'). We next present two numerical examples. The first one is the bistatic scattering from a homogeneous circular cylinder with a radius of 0.32A (A is the free-space wavelength), a relative permittivity of Or = 2.5 -j0.5, and a relative permeability of btr = 1.5 - j0.5. The computed scattering cross sections, cr/A as defined in [9], are shown in Figure 4, along with the exact eigenfunction solution. In the figure, VSIEM stands for the formulation presented in this paper, i.e., the Volume-Surface Integral Equation Method. Note that the two results are almost identical. The second example is the backscattering from an inhomogeneous rectangular cylinder of size 0.4A x 0.8A. Its relative permittivity is a function of position given by 6r = [1 + cos(rx/0.4A)][1 + cos(7ry/0.8A)]; and its relative permeability is a constant, ur, = 1.5. The results are shown in Figure 5, and are compared with those obtained using the Hybrid Finite Element Method (HFEM) [10], [11]. Again excellent agreement is observed. 12

5. 0. E-,H -l \ -'5. -10. -10- -........ VSIEM -15. -20. 0. 20. 40. 60. 80. 100. 120. 140. 160. 180. ) (degrees) (a) --- EXACT -0. -------- VSJEM -10. -20. -30. -40.,,,,,. 0 0. 20. 40. 60. 80. 100. 120. 140. 160. 180. ) (degrees) (b) Figure 4: Bistatic scattering pattern of a homogeneous circular cylinder, a = 0.32A, 6, = 2.5 - j0.5, r = 1.5 - jO.5. (a) TE case; (b) TM case. 13

0. -2. - HFEM -------- - VSiEM -4. -6. T -8. - 7t X -10. - -._ -12. 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. ( (degrees) (a) 10. 0. /- -10. m- HFEM -----—. VSfIEM -20. -30. -40. 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. 4 (degrees) (b) Figure 5: Backscattering pattern of an inhomogeneous rectangular cylinder, a = 0.4A, b = 0.8A, e, = [1 + cos(7rz/a)][1 + cos(7ry/b)], tr = 1.5. (a) TE case; (b) TM case. 14

It should be noted that the above VSIEM results were obtained with a sampling criterion of 12 points per material wavelength. Although the program we developed can be applied to arbitrary two-dimensional geometries and configurations, the above two examples, we believe, are sufficient to demonstrate the validity of equation (10) and the corresponding computer program. To show the contribution of various terms in (10) to the scattered field, we also made a computation for scattering by an inhomogeneous circular cylinder, whose geometry is the same as that illustrated in Figure 4. The cylinder has a radius of a = 0.4A, a relative permittivity of r = 2 - (r,/a)2 which is that of a slice through the center of a Luneberg lens, and a relative permeability of ir = 1.2. For the TE case, the scattered field is due to an equivalent volume magnetic current (2nd term in (10)) aend an equivalent volume magnetic dipole moment (3rd term in (10)). There is no surface dipole moment, since Cr becomes one at the surface of the cylinder, making the 4th term vanish in (10). The results are shown in Figure 6 in terms of the amplitude and phase of the scattered far-field coefficient P(O) [9]. For the TM case, the scattered field is due to an equivalent volume electric current (2nd term in (10)) and an equivalent surface electric dipole moment (4th term in (10)). In this case there is no volume dipole moment, since /r is a constant. The results are shown in Figure 7. From both Figures 6 and 7, we see that in the forward scatter region the field due to the volume currents is dominant and in-phase with the field that is due to the volume and surface dipoles; while in 15

3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0 --- I * Total field Field due to volume currents ----- Field due to volume dipoles - - -- - -- - - D. L..". I.. I..-. I - -.. a. -. 20. 40. 60. 80. 100. 120. 140. 160. 180. ~ (degrees) (a) 360. 300. u 240. 180. Im 120. 60. 0. ~ 0. ) (degrees) (b) Figure 6: Bistatic scattered field of an inhomogeneous circular cylinder for the TE case, a = 0.4A, = 2 - (r/a)2, t, = 1.2. (a) Amplitude; (b) Phase. 16

3.5 - 3.0 -_ _ _ _ _ _ 5 --- Total field 2.5 ' 5 --- —--- Field due to volume currents 2.0 ------ Field due to surface dipoles _ 2.0 -7 1.5 - 1.0 0.5 - 0.0 ' 0. 20. 40. 60. 80. 100. 120. 140. 160. 180. ) (degrees) (a) 360. 300. ) 240. -- - -.- ~~~^ _.. -- - - - - - - ---- 180.,,", i '<, I120. 60. 0. 20. 40. 60. 80. 100. 120. 140. 160. 180. 0 (degrees) (b) Figure 7: Bistatic scattered field of an inhomogeneous circular cylinder for the TM case, a = 0.4A, r = 2 - (r/a)2, Plu = 1.2. (a) Amplitude; (b) Phase. 17

the backscatter region both fields due to currents and dipoles have comparable magnitudes and they are out of phase. 5. Conclusion In this paper, an integral equation is derived for the analysis of two-dimensional electromagnetic scattering by inhomogeneous cylinders. The equation contains both volume and surface integrals but involves only one field component. From its physical interpretation it is seen that the field scattered by a cylinder is produced by a set of equivalent electric (or magnetic) current and dipole sources for TM (or TE) wave scattering. In regards to numerical analysis, the equation has two advantages: first, the equation contains only one unknown component, and thus is numerically more efficient than the volume integral equation formulation; second, the equation does not contain any derivative operation on the unknown, and hence the simple pulse expansion and point-matching technique can be used without loss of accuracy. Appendix I Derivation of (10) using the Method of Equivalent Polarized Currents In the method of equivalent polarized currents, the scatterer is replaced by an equivalent polarized electric current = j6o(6- 1)E (23) and an equivalent polarized magnetic current M = jwuo((rt- 1)H (24) 18

where E and H are' the total electric and magnetic fields, which are related by Maxwell's equations (V x E = -jW/oturH (25) i V x H = jwcocrE Let Ee and He represent the fields produced by J, and Em and Hm the fields produced by 4M. They obey the free-space Maxwell's equations V x E = -jwi/oHe (26) V x He = jwcoEe + J and V x Es = -jw:oHm - AM (27) V x H = jiwoEm Let us consider rTM\A wave incidence first. For this case, the electric field has only a z-component. The field due to J is Ee(r) = — juWo f j J(r')Go(r r')ds' - k J I[r(r') 1- ]Ez(r')Go(r-r'ds and the field due to M is Em(r) = x M(r )Go(r-rds J [ t(/ ' G0(l')' xE(rds V XJj x1 {[1- ] Go(rlr))E(r)}d JJ F[ 1 -- V x Go(rr')V7' ( x E(')ds' Go~ir)V ( X (28) (29) 19

Noting that VGo(rlr') = -V'Go(ri'), we can show that (29) can be expressed as Em() JJz v { [1 Ez(r;')V'Go(-I') ds' -V/1 ---.-' —' Using equation (2) and the two-dimensional divergence theorem, we can also write [(, z [1)/ol] E() for G Qd + Z [-l (r)]() for ron r (31) t0 for r OO The total electric field is the superposition of the incident field, the field due to the electric current J, and the field due to the magnetic current al Ez{r) = >(r) + z. Ee(r) + z. Em(i) (32) Substituting (28) an (31) into the above equation, e thus obtain J+ k [E()V', Go(rV ')]ds ' Ez - Ez VC()+ r. E,(r)+. E (/) (32) Ez ~~I~(,~) + k ~(~lol ~(7I')-> (, Ez(~')ol')' 20

+ - Ez (r) &Go(lt') d1' 1- [, _J ^ — W [1l/ir(r)]Ez(r) for r C Q [1 + /(r)]() for r on F (33) E,(r) for r G Q This equation is the same as (10) for TM incidence. The derivation for TE incidence is similar. Note that the boundary conditions on F are automatically satisfied in this method, while in the derivation of Section 2 they are explicitly enforced. Appendix II Matrix Element Expressions for (22) The element expressions for submatrices A, B, C, and D in equation (22) are given below: Ai. = o[v((r u —t(r)]H(k01olj- -rf) 40 ' -(2)- - -i +Jk [u(TU 3 + U ri ll3( t 2( -i r2s 4 I r \. L1 ' H3 i= 1,2,3,.,N Bij= 4ko[l-u(r) - )] ( [nk~). i + Y n,(rb Z. ( - ri-r|)l 4 L3r 3 i=1,2,3,...,N; j=1,2,33,,M J_ \1 ~)_~(~ H~o(~ -~ _ 21

+ [ ( ifr + Uy 1b (2)(ko rbj ) ~k [ i - -J1 i - rj = 1,2,3,,M; j = 1,2,3,..,N -, = - u(i)] [n() |b i + n()- b b ( - ij = 1,2,3,,AM; i 2 j Di3 = [l+1 u(r)] i =1,2,3,-,M The superscripts i and b are used to distinguish those variables in Q and on r, respectively. In the above expressions, sj denotes the area of the jth cell, lj the length of the jth segment, a - =, ur and uy are defined by Vu = ux: + Uyy, and n? and ny by 7/ = njr + nyi. H(2) and H(2) are the Hankel functions of the second kind, of the zeroth and first order, respectively. Acknowledgments The authors wish to thank Mr. Leland Pierce for his many helpful suggestions, and the JEWA reviewers for their comments. Part of this work (J.M.J. and C.T.T.) was supported by the National Science Foundation under Grant ECS-8319595. References [1] Richmond, J. H., "Scattering by a dielectric cylinder of arbitrary cross-section shape," IEEE aTans. Antennas Propagat., vol. AP-13, pp. 334-341, 1965. 22

[2] Richmond, J. H., "TE-wave scattering by a dielectric cylinder of arbitrary cross-section shape," IEEE Trans. Antennas Propagat., vol. AP-14, pp. 460-464, 1966. [3] Harrington, R. F., Field Computation by AMoment Methods. New York: Macmillan, 1968. [4] Newman, E. H., "ITM and TE scattering by a dielectric/ferrite cylinder in the presence of a half-plane," IEEE Trans. Antennas Propagat., vol. AP-34, pp. 804-813, 1986. [5] Ricoy, M. A. and J. L. Volakis, "Integral equations with reduced unknowns for the simulation of two-dimensional composite structures," The University of Michigan Radiation Laboratory Report No. 389492-2-T, Nov. 1987. [6] Tai, C.-T., "A note on the integral equations for the scattering of a plane wave by an electromagnetically permeable body," Electromagnetics, vol. 5, pp. 79-88, 1985. [7] Harrington, R. F., Time-Harmonic Electromagnetic Fields. New York: McGrawHill, 1961. [8] Hill, S. C., C. H.. Durney, and D. A. Christensen, "Numerical calculations of low-frequency TE fields in arbitrarily shaped inhomogeneous lossy dielectric cylinders," Radio Science, vol. 18, pp. 328-336, 1983. [9] Bowman, J. J., T. B. A. Senior, and P. L. E. Uslenghi ed., Electromagnetic and 23

Acoustic Scattering by Simple Shapes. North-Holland Publishing CompanyAmsterdam, pp. 6-7, 1969. A revised printing is published by Hemisphere Publishing Corporation, New York, 1987. [10] Jin, J.-M. and V. V. Liepa, "Application of hybrid finite element method to electromagnetic scattering from coated cylinders," IEEE Trans. Antennas Propagat., vol. AP-36, pp. 50-54, 1988. [11] Jin, J.-M. and V. V. Liepa, "A note on hybrid finite element method for solving scattering problems," IEEE Trans. Antennas Propagat., vol. AP-36, pp. 1486-1490. 1988. 24

A SIMPLE MOMENT METHOD PROGRAM FOR COMPUTING SCATTERING FROM COMPLEX CYLINDRICAL OBSTACLES Jian-Ming Jin and Valdis V. Liepa Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, Michigan 48109 (accepted for publication in Proc. Inst. Elec. Eng., part H, 1989) ABSTRACT This paper describes a moment method program for computing electromagnetic scattering (for both E- and H-polarization) from multiple perfectly conducting and penetrable inhomogeneous cylinders (The latter can be lossy and have both permittivity and permeability different from their free-space values). The program is a numerical implementation of a recently derived volume-surface integral equation and uses the pulse expansion and point-matching technique. The main feature of this program is its simplicity in generating the model and transforming the integral equation into a linear system of equations. Using this program, we explore and 25

demonstrate the usefulness as well as the limitations of the volume-surface integral equation formulation. I. Introduction In the frequency domain the problem of open region electromagnetic scattering is usually treated using integral equation methods (see e.g., [1]-[4]), partial differential equation methods (see e.g., [5]), and hybrid techniques which combine the partial differential equation methods with a surface integral equation or an eigenfunction expansion (see e.g., [6]-[9]). The integral equation methods have the advantages of simple numerical implementation and minimum discretization region, since the formulations incorporate the radiation condition by means of the free-space Green's function. However, they have the disadvantage of rather difficult formulation for complex media and result in full matrices which are costly to solve when th.ey are very large. The partial differential equation methods have the advantages of simple formulation even for complex media and simple numerical implementation, and produce sparse matrices which are much easier to solve. Their major disadvantage, however, is the need of extending the discretization region to the far field region in order to enforce the radiation condition, which usually leads to an extremely large number of unknowns'. The hybrid techniques represent a synthesis of these two methods: their advantages remain while their disadvantages are eliminated. As a result, they are more powerful especially for treating large scatterers. However, the formulations and especially the numerical implementation 1Currently, much attention is aimed at employing more complicated radiation conditions to reduce the discretization region. 26

of the hybrid techniques usually require more effort. This paper deals with two-dimensional electromagnetic scattering via the integral equation approach. The purpose of this paper is to describe a simple moment method program which is based on the formulation of a recently derived volumesurface integral equation (VSIE) [10]. Using this program, we will explore and demonstrate the usefulness as well as the limitations of the VSIE, and answer some questions raised about this new formulation and its applications. Electromagnetic scattering from and interaction with complex cylindrical obstacles of arbitrary cross section has recently attracted much attention from researchers due to its applications such as the control of radar cross section and the design of microwave devices. A surface integral equation solution has been presented by Arvas et al. for multiple perfectly conducting and homogeneous, lossy dielectric cylinders for transverse magnetic (TM) wave scattering [11]. A similar solution has also been presented recently by Yuan et al. for a homogeneous dielectric cylinder partially covered by conductors for both TM and TE (transverse electric) wave scattering [12]. Earlier, Newman treated a special case involving TM and TE scattering from a dielectric cylinder in the presence of a perfectly conducting half-plane using a volume integral equation with the half-plane Green's function as its kernel [13]. We will show that the formulation and its numerical implementation presented in this paper can effectively handle the most general case - TM and TE scattering from multiple cylinders which include perfectly conducting cylinders and lossy inhomogeneous penetrable cylinders with both permittivity 27

and permeability different from their free-space values. In Section II, the VSIE is presented for both E-polarization (TM) and Hpolarization (TE) fields in such a way that the physical meaning of the terms in the equation is more obvious. The application of the method of moments with pulse basis functions and point-matching to the VSIE is described in detail in Section III. Section IV gives some numerical examples to demonstrate the validity and versatility of the computer program. Following this is Section V discussing the limitations of the VSIE and the program. It is the authors' intent that with the description given in this paper the interested reader can easily develop his/her own program and further explore the VSIE application from different perspectives. II. Volume-Surface Integral Equation The volurne-surface integral equation (VSIE) is an integral equation governing two-dimensional electromagnetic radiation and scattering phenomena. This equation has been recently derived by Jin, Liepa, and Tai [10]. It differs from the conventional integral equations in that: (i) the VSIE contains both line and area integrals while the conventional integral equations usually contain only area integrals for inhomogeneous obstacles [2], [3] or line integrals for homogeneous obstacles [4]; (ii) the NVSIE has only one scalar unknown variable regardless of the material properties of the obstacles and the polarization of the fields, while the conventional volume integral equations have one or two or three scalar unknown variables depending on the permittivity and permeability of the obstacles and the polarization of the fields (see Table 1); and more importantly (iii) the VSIE does 28

Table 1. Numbers of Scalar Variables in the Volume Integral Equations (VIE) and Volume-Surface Integral Equation (VSIE) Material Ce 7c o, [ = Po e-= o,/ /i Po e 6 co, /t 7[t/o Field E'-pol. H-pol. E-pol. H-pol. E-pol. H-pol. VIE 1 2 2 1 3 3 VSIE 1 1 1 1 1 1 not involve any derivative operation on the unknown while the conventional volume integral equations have derivatives operating on the unknowns for the general case. The second and third differences give the VSIE a notable advantage over the conventional volume integral equations. In this section, we will present, without derivation, the VSIE while the detailed formulation can be found in [10]. The problem under consideration is illustrated in Figure 1, where a source excites an electromagnetic field in the presence of some obstacles which consist of inhomogeneous dielectrics2 and perfect conductors. The source and the obstacles are of infinite extent in the z-direction, and have no variation of any kind with respect to the z-coordinate. Thus, the resultant field under such circumstances has no variation with respect to the z-coordinate. In the two-dimensional case which we consider here, an arbitrary electromagnetic field can be expressed as the sum of an E-polarized field which has only a z2By the word dielectric, in this paper we refer to the general penetrable material rather than a kind of particular material with r = 1 - 29

I —i S x 44 m -1 i- - - Source V/ /// Conductor 1:::j Dielectric Figure 1: An illustration of the problem under consideration. 30

component of electric field and an H-polarized field which has only a z-component of magnetic field. Therefore, we can solve the problem by treating the E-polarized and H-polarized fields separately. In the following, we state the VSIE for these two fields separately. A unified presentation of the VSIE for the two fields is possible if appropriate notations are used [10]. E-polarization (TM) The total z-directed electric field, Ez, is a superposition of the incident field, EZ, due to the excitation and the scattered field, E<, due to the presence of the obstacles: E(r) = E'(r) + E(r). (1) For the problem under consideration, the scattered field consists of five parts due to five different sources: 5 E(r)= E()(r). (2) i=l 1. The first source is a volume electric current due to the non-unity of the permittivity in the region, Q, occupied by the dielectric obstacles: Jz(') - jo[r(r')- -]Ez(r'). (3) The Green's function3 applied to this source is G(rr') = -jwaoGo(frt') (4) 3By the phrase Greern's function, here we mean the electric field produced by a unit source. 31

where Go(rfr') = -(j/4)H(2)(kor - r'l). Thus, the source produces a field ( ko [c (r') - 1]Ez(r')Go(I')ds'. (5) 2. The second source is also a volume electric current; but, it is due to the nonunity of the permeability: = o [ - (r E ' ). (6) The Green's function applied to this source is G(r) =- o + Go(rlr (7) and the field produced by the source is then Es(2)(r) 1() + k f - E(')Go(r)ds'. (8) V(1) Jr or) 3. The third source is a volume electric dipole moment due to the continuous variation of the permeability: W/~o. I(9) pV( ) = Ez(-')V' Kj F) (9) which is a transverse vector. Its Green's function is also a transverse vector: G(r r') = -jwoV'Go(rlr') (10) and gives E(3)() = Er' r(11) 32

4. The fourth source is a surface electric dipole moment induced on the interfaces denoted by rd where the permeability changes abruptly (including the air-dielectric interface): pS = / ) - () E,(), ) (12) where "n is a unit vector normal to Fd and pointing from the "-" side to the "+" side. The Green's function for this source is the same as (10). Thus, the source produces a field (r) '[ [ T() - (- ) Ez(r ')nd. V'Go(rr ')dl'. (13) 5. The last source is a surface electric current induced on the conductors' surfaces denoted by F,: z(') -. V'Ez(i ') (14) jWLO/r(r') where n' is the outward unit vector normal to r,. This current produces a field E()'(r') =G0( 1)1. V'z(')dl. (15) Note that in this case (E-polarization) the variation (either continuous or abrupt) of the permittivity does not induce any kind of dipole source. When all the above sources are combined, the VSIE for the E-polarization case becomes kJ [Cr - ] (') Go( )d'+ / J (i') V '' (V'G(l')d' 33

[1 11 aGnI'dl- [ 1 (I-)EZdl' I Ltr(+) -r(K ) ) nd C Pr() nr' c [1//i(r-)] Ez(f) for r in dielectrics +E(r) = 0 for r in conductors (16) 1 Ez(r) for r everywhere else. Il-polarization (TE) The formulation for the H-polarization case is dual to the above for the Epolarization case. For brevity, we only give the equations: H (r) = H(r) + Hz () (17) where 5 Hz(r)- Hs)(r) (18) i=1 with -(- (r) = 2 [ir') - l]Hz(r')Go(rjr')ds' (19) (volume magnetic current contribution) Hs(r)= [1- ( )]H(r) + k / [1 - +H(r')Go(rr')ds' (20) (volume magnetic current contribution) H )(r) =(r)V [ ] V'Go-(r ')ds' (21) (volume magnetic dipole contribution) 34

H r) () r (~ ) Hz( r' nd V'Go (r ldl' (22) (surface magnetic dipole contribution) H r= / ()Hz(r')n, VGo(rlr')dl' (23) (surface electric current contribution). In contrast to the E-polarization case, here the variation of the permeability does not introduce any kind of dipole source. The VSIE for the //-polarization case can then be written as [r(') () Hz(')Go(rI')d'+ J H(ri')V',) V'Go(r l')rds' / - 1 D GO(( ifI),, Jrd [Er(K+) - r(l')l H(I)Dnu Jrc dl'~rJ( H (ot f [1/6r(f)] Hz(r) for r in dielectrics +H(r) = < 0 for r in conductors (24) Hz() for r everywhere else. III. Application of the Method of Moments The VSIE of (16) for E-polarization and (24) for H-polarization can be numerically solved using the method of moments (MM) [14]. To ensure an easy generation of the model, we adopt the simple pulse expansion (basis) functions and delta testing (weighting) functions (point matching) for the MM solution. In this section, we will describe the procedure of applying this pulse expansion and point matching technique to the VSIE. 35

Discretization of Obstacles The first step of the application is to discretize the obstacles, i.e., to divide the obstacles into a set of cells and segments. Such is illustrated through an example shown in Figure 2. Figure 2(a) shows an obstacle consisting of a half-circular dielectric cylinder and a rectangular dielectric cylinder with a smaller rectangular conducting cylinder embedded inside. At the interface between the half-circular cylinder and the rectangular cylinder, there exists a step discontinuity in both permittivity and permeability. Figure 2(b) shows a discretization model of the obstacle. The region Q, occupied by the dielectrics, is divided into a set of triangular and quadrilateral cells, and the interface of discontinuity Fd and the conductor's surface Fc are divided into a set of line segments. The numb)er of discrete points determines the number of the discretized unknowns in the final system of equations. Assume the number of cells in Q is N, the number of segments on Fd is M, and the number of segments on Fc is L. The total number of discretized unknowns is given in Table 2 for different properties of the materials and different polarizations of the incident fields. The above discretization model is easy to generate. Since, as we will show later, the discretization of the VSIE only needs the area and the center-point position of each cell, one does not have to be concerned with the shape of the cell. This substantially simplifies the generation of models. 36

(a) 0 0 ~ *a; (b) Figure 2: (a) An original problem; (b) Its discretization model. 37

Table 2. Number of Discretized Unknowns in the System Matrix Equation Material e o, /t = o 6 = 60,: # to $e 7- eo, 7# [to0 Field E-pol. H-pol. E-pol. H-pol. E-pol. H-pol. Unknowns N+ L N+M+L N+M+L N++ L N+M+L N+M+L N - Number of cells in the regions (Q) occupied by dielectrics 1 - Number of segments on the interfaces (Fd) between dielectrics L -Number of segments on the surfaces (FU) of conductors To generate the points on the interfaces and the conductors' surfaces, we first break the curves of the interfaces and the surfaces into several arcs with straight lines as a special case. We then input the end-point coordinates of and the angle subtended by each arc, and the number of segments we wish to put on the arc. The program generates the length, the mid-point coordinates, and the normal sine and cosine of each segment. To generate the points inside the dielectrics, we first divide the cross-sections of the dielectrics into a number of curved or planar strips each of which has a constant thickness. Then we input the thickness of each strip, the end-point coordinates of and the angle subtended by the center-line of the strip, and the number of cells we wish to use. The program will generate the area and the center-point coordinates of each cell. Such a model generation can be accomplished by a rather small program. Wie found that the above input format is very convenient, and with care we can generate a quite accurate model for an arbitrary geometry. This is one of the reasons for our choice of pulse expansion and point matching. Other expansions 38

such as linear expansion requires generating more complicated models such as the triangular patch model. Such models are usually difficult to generate, since they need detailed information about the shape of each cell. Though there are some methods which are capable of doing such a job, they usually require quite large program packages. Discretization of the VSIE The second step of applying the MM is to discretize the VSIE of (16) for E-polarization and (24) for H-polarization. For brevity, we only consider the Epolarization case here. The discretization of (24) for the H-polarization case is similar. Let us consider the formulation for the discretization model shown in Figure 2(b). WVe first approximate the continuous unknown field in terms of a set of pulse functions. For the field in Q, we have N E,(r') = Pi (25) i=1 where Pi is a pulse function which equals 1 in cell si and equals 0 everywhere else, and a'i is an unknown coefficient. For the field on Fd, we have M Ez()= j (26) j=l where Pj is a pulse function which equals 1 on segment Ij and equals 0 everywhere else, and Oj is an unknown coefficient. Similarly, for the field on F~, we have 1 OE (W') L L= -7kPk (27) (rP ) On/ k=l 39

where Pk is a pulse function which equals 1 on segment Ck and equals 0 everywhere else, and 7k is an unknown coefficient. Substituting (25)-(27) into (16), we have ci oj { [r(r')- o G0(I') + V V'Go ') ds' + [7 [- - ( )dl- Ykj Goo(l')dlr [1//I(r)] EiN=i cPi for r in Q +EZ(r) -= (1/2) [l/I/r(r+) + 1//r(r-)] EjM= 3jPj for r on rd (28) 0 for r on Fc where f denotes the Cauchy principle value integral with singularities removed. By satisfying (28) at each center-point of N cells and mid-point of M and L segments, we obtain N + Ml + L linear algebraic equations whose solution gives values for the N + AM -+ L unknown cofficients. The scattered far field coefficient, P(4) defined in [15], is then obtained by P ( ). zi Ci (4) Z ii { [ r(Li) - tr(i) 4 [4OX r(r) = oi + g () sinj Je ^y }Or(r) scosb+yi sin(J) +- yX jLj1 [ (f - ~( ) (cos 6 cos + sin 0 dsin ) [[4r~rj) -- L.e Jko(x cos ++ yj sin )) + 3 )y Cejk o(xk COS +yksin fl) (2 ) k=l where ( is the observation angle, Si is the area of cell s, L- is the length of segment Ij, Ck is the length of the segment Ck, and 0d is defined by nd = x cos 0d + Y sin d. 40

For a two dimensional geometry the radar cross section (RCS) per wavelength is a(M lp( 12 (30) A 7r with A being the free-space wavelength. Before we close this section, we would like to address the evaluation of some of the integrals. To cast (28) into the final system matrix equation form, one needs to evaluate the following integrals 11= I ko [(r -' - ] G (,mr')ds' (31),= I[ -r()] 2O r K dl' (33) 14 = G0(rm r')dl' (34) Ok where rm = r (i = 1,2,..-,N), or rm = rj j = 1,2, *,M), or rm = k (k = 1,2,...,L). For evaluating the area integrals '1 and 12, if r-i $ ir,, we use the approximation I = [6r(ri) - tr() Go(r i) * Si (35) I2 = {V -(r] VG~o(rm Ir-) S-. (36) If Tr = im, we first approximate cell si as a circular cell of an area equal to that of the original cell and then evaluate the integrals analytically. The results are I - [6r,(i) - (/(i) + joai (koai) (37) '2 = 0 (38) 41

where as = (,S/7r)1/2. It may be worth pointing out that (35) and (36) are approxinmate from the perspective of the pulse function expansion. However, if we use the delta function expansion N Ez(r) = E Si6(r - r) (39) i=l for the non-self or off-diagonal elements, (35) and (36) become exact. Therefore, to be more precise, the expansion we use is a pulse function expansion for the diagonals (ri = rm) and a delta function expansion for the off-diagonals (rl: ri). The line integrals 13 and 14 can be approximately evaluated as = [1/r(;j+) -1/tr(j —)] aGo(rm\r)/O7zdlr=j Lj for rj m r ( I= - (40) 0t for rj = rm Go(r;mlrk) Ck for rk - rm 14 - (41) (-j/4) [1 - jlog (ykockl/4E)] Ck for rk = rm( where 7 = 1.781 and e = 2.718. It should be pointed out that the above integral evaluation is a crude approximation. However, it has the advantage of resulting in simple analytical expressions, and more importantly, the evaluation of A1 and 12 needs only the central point coordinates and the area of each cell, rather than detailed information about the shape of the cell. As pointed out earlier, such an approach substantially simplifies the work of the discretization of the obstacles. The discretization of (24) for H-polarization is similar to that described above, since (24) is dual to (16), except for the last integral which is along the conductors' 42

surfaces. This integral requires the evaluation of 15 1 X G0(imf') (42) ck aOnc and the result is 0 Go0(r-mr)/on\ c=k Ck for rk 7r rm I5 - (43) { 2 for rk = rm. IV. Numerical Results In this section, we present some numerical results computed using the moment method (MM) program based on the formulation and discussion given above. Seven different geometries are considered to demonstrate the validity and versatility of the program. The first geometry is a coated conducting circular cylinder, which can be solved analytically using an eigenfunction expansion technique. Such an exact solution can then be used to verify our VSIE MM results. In Figure 3 we present both eigenfunction and VSIE MM solutions for the backscattering radar cross section of the cylinder as a function of the coating thickness for three different coatings. Overall, the VSIE MM solutions agree quite well with the eigenfunction solutions. The second and third geometries are related to radar targets. The second geometry is an inhomogeneous dielectric circular cylinder backed by a conducting strip, as shown in Figure 4(a). The cylinder has the permittivity variation in radius the same as that of a spherical Luneberg lens. The third geometry, as shown in Figure 5(a), is an inhomogeneous dielectric half-circular cylinder backed 43

10... 5. 0. -5., = 2.4-jO.O, r = 1.2-jO.O -25. - --------- e,=2.4-jl.0,l = 1.2-jO.O ------ e,=2.4 -jl. =1. 2jl.O 0.00 0.05 0.10 0.15 0.20 0.25 0.30 tA (a) 10 T. -- 2.4-j0.,. = 1.2-jO.O --------- E, = 2.4-jl.0, r, = 1.2-jO.O ------ = 2.4-jl.0O,.=1.2.jl.O 0.00 0.05 0.10 0.15 0.20 (.25 0.30 0. -5. \ I/ -10. -15. ----- -20. -25. -30...... 0.00 0.05 0.10 0.15 0.20 0.25 0.30 (b) Figure 3: Backscattering cross section (o) of a coated conducting circular cylinder vs coating thickness (t). The radius of the conducting cylinder equals 0.2A. The lines represent the eigenfunction expansion solutions. The circles respresent the VSIE MM solutions. (a) E-polarization; (b) H-polarization. 44

10. 5. 0. ) -5..-10... -20.,.. I... I... 1.. I I...!... I... I... 0. 20. 40. 60. 80. 100. 120. 140. 160. 180. 0 (degrees) (a) 10. - Strip and cylinder ", \ \.......... Cylinder only O. ------ Strip only -20. - 0 0. 20. 40. 60. 80. 100. 120. 140. 160. 180. 0 (degrees) (b) Figure 4: Backscattering pattern of a cylindrical Luneberg lens backed by a conducting strip. a = 0.4A, oa = 90~, r - 2 - (r/a)2, ['r = 1. (a) E-polarization; (b) H-polarization. 45

10. 0 (degrees) (a) 10. ~., v * ';,I -20. - \ -30. 0. 20. 40. 60. 80. 100. 120. 140. 160. 180. 0 (degrees) (b) Figure 5: Backscattering pattern of a half cylindrical Maxwell fish eye backed by a conducting strip. a = 0.6A, a = 60~, Cr = 4/[1 + (r/a)2]2, r = 1. (a) E-polarization; (b) HI-polarization. 46

by a conducting strip. This cylinder has the permittivity variation in radius same as that of a Maxwell fish eye. In the figures we show the computed backscattering radar cross sections of the dielectric cylinder, conducting strip and the combined structure. The fourth geometry, shown in Figure 6(a), is a 3A wide conducting strip having its edges loaded with dielectric circular cylinders. The results presented are the backscattering radar cross section of the strip with and without loadings. 'We should note that for the E-polarization computations shown in Figures 4(a), 5(a) and 6(a) the conducting strips are infinitely thin, while for the Hpolarization computations of Figures 4(b), 5(b) and 6(b) the strips are assumed to have a finite thickness of 0.06A in order to allievate the difficulty encountered in the H-polarization formulation for thin strips of zero thickness. We should also note that the problem of electromagnetic scattering from microstrip structures, such as those treated in [11] and [12], is similar to that of Figures 4 and 5, and thus can be treated in the same manner. The last three geometries concern the application of the program to the calculation of the scattering from arbitrarily shaped conducting cylinders coated with dielectrics. The geometries are shown in Figures 7(a), 8(a) and 9(a), and they are coated square, ogival and wedge cylinders, respectively. The results shown are the backscattering radar cross sections of the cylinders for different coating thicknesses and for different permittivities and permeabilities of the coatings. Also given in Figures 7 and 8 are the results obtained using the hybrid finite element method 47

0 (degrees) (a) 20. 0 (degrees) (b) Figure 6: Backscattering pattern of an edge-loaded conducting strip. w = 3.0A, d 0.3A, r = 4.0 - j4.0, r = 2.0 - j2.0. (a) E-polarization; (b) H-polarization. 48

5. 0. -10... -15. - I I I I I 0. 5. 10. 15. 20. 25. 30. 35. 40. 45. * (degrees) (a) 5.. — 5. ---- t=0.03-. ------ t= 5.:..-..I................... I.. -5. 9 --10. = "......- t.00......... t=0.031. t=0.063. 0. 5. 10. 15. 20. 25. 30. 35. 40. 45. * (degrees) (b) Figure 7: Backscattering pattern of a coated square cylinder. a = 0.5A, r - 2 - j2,r = 2 - j2. The lines represent the VSIE MM solutions. The circles represent the HFEM solutions. (a) E-polarization; (b) H-polarization. 49

10. 5. 0. -5. -10. -15. 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. * (degrees) (a) 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. * (degrees) (b) Figure 8: Backscattering pattern of a coated ogival cylinder. w = 1.2A, a = 60~, t = 0.06A. The lines represent the VSIE MM solutions. The circles represent the UIFEM solutions. (a) E-polarization; (b) H-polarization. 50

10. -20. 0. 20. 40. 60. 80. 100. 120. 140. 160. 180. 0 (degrees) (a) 10. 0. -10. -:20. -:30. 0. 20. 40. 60. 80. 100. 120. 140. 160. 180. * (degrees) (b) Figure 9: Backscattering pattern of a coated wedge cylinder. a = 0.3A, a = 40~, t = O.G06A. (a) E-polarization; (b) H-polarization. 51

(HFEM) [9]. A generally good agreement between the VSIE MM and HFEM solutions is observed. Before closing this section, we would like to point out that though we have not encountered any major problem in obtaining the above results, we found that for the /I-polarization computations of thinly coated structures the number of discrete points or cells needed to produce a given accuracy is usually larger than that for the E-polarization computations for the same accuracy. This fact is due to the existence of surface waves for the H-polarization case which may produce a complicated field distribution along the circumference of the cylinder. As an example, for a coated circular cylinder with a coating thickness of 0.03A, a coating permittivity of cr = 2 and peremeability of j, = 2, computations show that for the E-polarization case a sampling rate of 12 points per free-space wavelength along the circumference is enough to obtain the backscattered far field with an accuracy of O.1dB in magnitude and 1 degree in phase, while for the H-polarization case one would increase the sampling rate to 22 points per free-space wavelength to achieve the same accuracy. V. Some Limitations Like any other integral equation and numerical method, the VSIE and the program described in this paper also have some limitations. Regarding the main limitations, we note the following four points. * For using the VSIE, one needs to specify the gradient of permeability for Epolarization and the gradient of permittivity for H-polarization. For some prob 52

lems where such a specification is difficult or impossible, the VSIE cannot be used. For such a case, one would turn to the volume integral equation method where only the values of permittivity and permeability are required. * Because of the presence of the line integral in the VSIE, the resulting linear system of equations does not inherit the convolutional form; hence, the conjugate gradient-fast Fourier transform method, which is an efficient means of solving the pulse expansion-moment method linear system, may not be applied easily. * The discretization model used here requires calculating the field values on the interfaces of all discontinuities as unknowns. For cylinders containing many internal step discontinuities, the number of unknowns will increase significantly. * The program is not suited for the calculation of internal fields of cylinders with large values of permittivity for the TE or I/-polarization case, though it is very efficient for the TM or E-polarization case. Of the above four points, the first and second and partly the fourth result from the integral equation itself. The third would not be a drawback if one uses patch models where unknowns are assigned on the patch boundaries. The fourth could be possibly removed if one adopts more accurate models and evaluates the integrals more accurately. More work needs to be done for such cases. VI. Conclusion In this paper, a moment method program is described for computing electromagnetic scattering from multiple perfectly conducting and penetrable inhomogeneous cylinders for both TM and TE wave incidence. The program is based on the 53

formulation of the volume-surface integral equation, and has the advantages of easy generation of the model and simple discretization of the integral equation. These two advantages result from the employment of the pulse and delta expansion and point-matching technique. The pulse and delta functions are in the domain of the VSIE operator, while they are not in the domain of the volume integral equation operator for the general case. The program is useful for solving the complicated two-dimensional scattering problems involving perfect conductors and lossy inhomogeneous materials, and may find applications in the areas of scattering control, radar identification and remote sensing. Concerning the three-dimensional case, we note that a similar equation has been formulated by Tai for scattering by homogeneous, electrically and magnetically permeable bodies and inhomogeneous bodies with bLr = 1 [16]. It would be beneficial if one would derive a similar VSIE for the general three-dimensional case. Acknowledgments The authors wish to thank Drs. D. G. Dudley, C. H. Durney and E. H. Newman for their valuable comments on this work, and Mr. Leland Pierce for many helpful discussions. References [1] Mei, K. K. and J. Van Bladel, "Scattering by perfectly conducting rectangular cylinder," IEEE Trans. Antennas Propagat., vol. AP-11, pp. 185-192, 1963. 54

[2] Richmond, J. H., "Scattering by a dielectric cylinder of arbitrary cross-section shape," IEEE Trans. Antennas Propagat., vol. AP-13, pp. 334-341, 1965. [3] Richmond, J. H., "TE-wave scattering by a dielectric cylinder of arbitrary cross-section shape," IEEE Trans. Antennas Propagat., vol. AP-14, pp. 460-464, 1966. [4] Wu, T. IK. and L. L. Tsai, "Scattering by arbitrarily cross sectioned layered lossy dielectric cylinders," IEEE Trans. Antennas Propagat., vol. AP-25, pp. 518-524, 1977. [5] Mason, J. L., "Finite element solution for electromagnetic scattering from twodimensional bodies," Ph.D. dissertation, The University of Michigan, Ann Arbor, 1982. [6] Chang, S. K. and K. K. Mei, "Application of the unimoment method to electromagnetic scattering of dielectric cylinders," IEEE Trans. Antennas Propagat., vol. AP-24, pp. 35-42, 1976. [7] Marin, S. P., "Computing scattering amplitudes for arbitrary cylinders under incident plane waves," IEEE Trans. Antennas Propagat., vol. AP-30, pp. 1045-1049, 1982. [8] Jeng, S. K. and C. H. Chen, "On variational electromagnetics: theory and application," IEEE Trans. Antennas Propagat., vol. AP-32, pp. 902-907, 1984. 55

[9] Jin, J. M. and V. V. Liepa, "Application of hybrid finite element method to electromagnetic scattering from coated cylinders," IEEE Trans. Antennas Propagat., vol. AP-36, pp. 50-54, 1988. [10] Jin, J. M., V. V. Liepa, and C. T. Tai, "A volume-surface integral equation for electromagnetic scattering by inhomogeneous cylinders," Journal of Electromagnetic Waves and Applications, vol. 2, no. 5/6, pp. 573-588, 1988. [11] Arvas, E., S. M. Rao, and T. K. Sarkar, "E-field solution of TM-scattering from multiple perfectly conducting and lossy dielectric cylinders of arbitrary cross-section," IEE proc. H, Microwaves, Antennas and Propagation, vol. 133, pp. 115-121, 1986. [12] Yuan, X., R. F. Harrington, and S. S. Lee, "Electromagnetic scattering by a dielectric cylinder partially covered by conductors," Journal of Electromagnetic Waves and Applications, vol. 2, no. 1, pp. 21-44, 1988. [13] Newman, E. H., "TM and TE scattering by a dielectric/ferrite cylinder in the presence of a half-plane," IEEE Trans. Antennas Propagat., vol. AP-34, pp. 804-813, 1986). [14] Harrington, R. F., Field Computation by Moment Methods. New York: Macmillan, 1968. [15] Bowman, J. J.,. B. A. Senior, and P. L. E. Uslenghi ed., Electromagnetic and Acoustic Scattering by Simple Shapes. North-Holland Publishing Company 56

Amsterdam, 1969. A revised printing is published by Hemisphere Publishing Corporation, New York, 1987. [16] Tai, C. T., "A note on the integral equations for the scattering of a plane wave by an electromagnetically permeable body," Electromagnetics, vol. 5, pp. 79-88, 1985. 57

A MOMENT METHOD SOLUTION OF A VOLUME-SURFACE INTEGRAL EQUATION USING ISOPARAMETRIC ELEMENTS AND POINT-MATCHING Jian-Ming Jin, John L. Volakis, and Valdis V. Liepa Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, AMichigan 48109 (submitted to IEEE Trans. Microwave Theory Tech. for publication) ABSTRACT It is shown that traditional subdomain elements such as rectangles and triangles with pulse expansion basis could lead to considerable inaccuracies when simulating biological scatterers having high permittivities. In this paper, isoparametric elements are used in amomentarin method implementation to remove modelling inaccuracies of fields and boundaries associated with traditional elements. Numerical results are also given that show the improvement achieved in the scattering solution for high contrast circular cylinders. 58

1. Introduction A Volume-Surface Integral Equation (VSIE) was recently presented [1], [2] for electromagnetic scattering by inhomogeneous cylinders. A moment method implementation of the VSIE was also considered using rectangular and/or triangular subdomains (elements) with pulse expansion basis functions and point-matching. Recently, however, we observed that such a moment method solution becomes inaccurate in the case of scatterers having high permittivity with transverse electric (TE) incidence or high permeability with transverse magnetic (TM) incidence. This is primarily true for cylindrical scatterers associated with curved boundaries such as those encountered in biological models. Under these circumstances, the simple rectangular and/or triangular elements with pulse expansion basis do not render an accurate modelling of the internal fields and the scatterers' boundaries. This model inaccuracy is generally forgiving in the case of scatterers with nominal values of refractive indices but it will be shown to be too compromising when the refractive index becomes large as is usually the case with biological scatterers. Alternative discretization elements are therefore required to improve the scatterer's modelling accuracy. The isoparametric elements [3] are found to be capable of serving this purpose and in this paper we present a numerical solution of the VSIE by employing such elements. The isoparametric elements were first introduced in finite element analysis [3] and refer to discretization elements of the geometry having the same order (shape) as the field representation within them. That is, the field expansion/basis over a 59

quadratic line element, for example, is also quadratic. The main advantage of using isoparametric elements is to allow an accurate modelling in the case of curved arbitrarily shaped geometries. However, it appears that only recently [4] they have been employed in applications related to solutions of integral equations for electromagnetics. Below we first discuss the inaccuracy associated with traditional solutions of integral equations for penetrable scatterers having high refractive indices. This is followed by the introduction of the isoparametric elements and the presentation of a moment method solution of the VSIE using such elements. Results are subsequently presented which reveal the stability of the solution in the case of cylindrical geometries associated with high refractive indices. 2. Discussion of Integral Equations for TE Scattering In this section we examine the three integral equations for the computation of the internal fields in a dielectric cylinder for the TE incidence. Assume that the cylinder has its infinite dimension along the z-axis. The VSIE for this case is [1] H Cr ~[ k() i H ') Go( ) ')d' + d V' [Hz (r')V' ( lr')]ds' + - H(') Go( ) dl' { [1/er(r)] H(r) for r not on r (1/2) [1/6r(r) - 1/r(r')] Hz(r) for r on (1) where HINC denotes the z-component of the incident magnetic field, Hz represents the z-component of the total magnetic field, Er denotes the relative permittivity of the scatterer and Go(rlr') is the two-dimensional free-space Green's function. In addition, Q denotes the region occupied by the cylinder, F denotes the interface 60

where c,(r) has a step discontinuity, n is the unit vector normal to F pointing from the "-" side to the "+" side and f- denotes the Cauchy principle value integral. Alternative integral equations for this problem were also given by Borup et al. [5] as Et(r) = Et (i) + ko J J [r(r') - ] E(r')Go(rlr')ds' - j v [(e )- I)Et(r')] V'Go(rrI)ds (2) and recently by Peterson and Klock [6] as () - HINC(r) + JJ[z V' x 7(r')]GO(r')ds' (3) with (r) = [1 - V x [zHz(-)] and Et denoting the electric field transverse to the z-axis. As usual, Et represents the incident electric field transverse to the z-axis. Let us now examine the numerical implications which may arise when either of the three integral equations are implemented in the case of large I |r. Referring first to (1), we may assume without loss of generality that the incident wave, HINC, has unity amplitude. As a result the total internal field, Hz, will also have an amplitude around unity and thus the value of the right hand side of (1) will be about 1/|erl for a field point inside the cylinder. For large eCr\ this implies that the integrals over Q and F must together give a value nearly equal to and negative of HINC. That is, if Ce Ir = 100 and we demand a 1% solution accuracy for the internal 61

fields, the integrals must be computed with a corresponding accuracy of 0.01%. In general, the presence of \Er\ acts as an error amplifier and we may conclude that for a solution accuracy of 6% the integral must be evaluated with a corresponding accuracy of 6/cr\l percent. This statement has been experimentally verified and, clearly, for large cEr a crude discretization of the scatterer is unlikely to produce the accuracy demanded here. A similar examination of (2) also reveals that for a desired solution accuracy of 8%, the integral must be evaluated with a corresponding accuracy of 6/m percent, where 7n = \EIN /I \E\ - I/v/4f. Thus, even though (2) is slightly less demanding in accuracy, it is still unstable for scatterers with large refractive indices. This exposition satisfactorily explains the difficulties encountered in [5] as well as the success of their modified approach. Also, the results given by Boyes and Kennedy [7] can be explained in this manner. At first glance, (3) appears to be the most attractive of all three integral equations. It is certainly more stable than (1) and (2); however, it involves a second derivative of the unknown field quantity and as a result demands higher order basis functions than (1) and (2). In particular, (1) can be implemented with pulse basis and (2) with linear basis but (3) requires quadratic basis1. Since the introduction of quadratic basis in (1) and (2) will improve their accuracy, there is no clear choice which of the three is more attractive for numerical implementation. Each is certainly associated with its own advantages and disadvantages for large refractive If ~r is uniform within each element, the corresponding basis functions can be one order lower. 62

indices but for nominal values of ~r the VSIE given in (1) is generally superior over the others. In the next section we consider an implementation of the VSIE using isoparametric elements. As stated in the introduction these provide a more accurate geometrical simulation of the scatterer and representation of its interior field distribution. They are, therefore, expected to provide the required accuracy for the simulation of high contrast dielectrics. 3. Formulation of Isoparametric Elements and Point-Matching Technique In this section we describe a numerical solution of the VSIE using isoparametric elements with point-matching. Let us consider both TE and TM wave scattering by an inhomogeneous cylinder having non-unity permittivity and permeability. From [1], the field satisfies the general VSIE FprN () + k J // [v( () - u(r')]F(r')Go(r r')d i' + j Vtu(r') I [F(r')V'Go(rr')]ds' + - [u(r+ - u()] F(r) Goi ' dl.n' u(r)F(r) for r not on F (4) t [u(r+) + u(r_)]F(r) for r on r where F(i) = Ez,(), u(r)- /= (), v(r) = (r), for TM incidence and F(T) = Z(r)), u(r) = ( ) v() t r(r). 63

for TE incidence, in which Ez represents the z-component of the electric field and ir denotes the relative permeability of the scatterer as usual. For convenience, we can write (4) as FINC(i) + J j A(rir')F(r')ds' + -IB(Frr')F(')dl' u(r')F(r) for r not on 1F (5) = [u(r+) + u(r_)]F(r) for r on r where the expressions for A(rlr') and B(rIr') are determined by direct comparison with (4). In a numerical solution of (4) or (5) using isoparametric elements, the twodimensional region Q is broken up into a number of quadrilaterals whose sides can be curved, as shown, for example, in Figure 1. A usual rule for this subdivision is that the interface F coincides the with boundries of the quadrilaterals. Assuming now that the discretization results in M2 quadrilaterals and that MI curved segments make up the contour(s) F, equation (5) becomes M.r2 n2.M1 nl F r) + E E J A(rr )V r )ds+ E> a y/1 B(rr')U(r')dl' e=:1 i J=i s=1 J= ' u(O)F(r) for F not on F (6) t [u(r+) + u(r_)]F(r) for r on F In (6), the field in the eth quadrilateral has been expanded as n2 Fe(i) = ZV e((F)r (7) i=1 where Vie(r) (i = 1, 2,...,n2) are a set of known n expanson basis functions and (i = 1,2,..., n2) are the corresponding unknown coefficients. Similarly, the field in 64

(a) (b) Figure 1: Two examples for breaking a circular region. (a) A 5 quadrilateral element model with 20 nodes. (b) A 12 quadrilateral element model with 45 nodes. 65

the sth segment was expanded as nl F(r) = U()s (8) i=i where Us(r) (i = 1,2,...,nl) represent the known expansion basis functions and q5 (i = 1, 2,..., nl) are the unknown constants. Before proceeding with a solution of (6) it is first necessary to choose the basis functions V1e(r) and U's(r) in accordance with the definition of the isoparametric elements. A possible method of constructing Us and Ve is to choose them so that qid (i = 1, 2,..., nl) represent the field values at nl nodes on the sth segment and likewise X (i = 1,2,..., n2) represent the field values at n2 nodes of the eth quadrilateral. A point-matching procedure applied to (6) will then lead to a matrix equation for the solution of the nodal field values. To find Us as described above it is convenient to introduce the transformation 3 3 x = ZL()x, y = L8()yi (9) i=l i=1 allowing a linearization of a quadratically curved segment as shown in Figure 2. In (9), Ls (i = 1,2, 3) are the shape functions which take the well-known form [3] 1. 1 L =-s (1-) ( +, i - -. 2 In addition, (9) implies the relation dl- ~ + < - IJJl`j (I0) and the unit vector n is now given by ( - - I I ( 11) 66

1 3 2 -1 o +1 x 0 (a) (b) Figure 2: A quadratically curved segment in the xy-p)lane (a) can be transformed into a straight segment lying on the (-axis (b). 67

It is observed that L'(0) is unity at the ith node and thus by choosing U^ = L 6s7 will coincide with the nodal values of the field FS(r). This actually is the definition of the isoparametric elements since Ls describe the geometrical shape of the discrete elements. When this expansion is substituted in (6), it leads to integrals of the form s = B( ')L ')dl' = J B(')d' (12) and these can be evaluated via a four point Gaussian integration formula giving 4 I = E WjB(r r-')Ls( J)Ijfs((j) (13) ^-E = ^ >3 ^ ^ U (13) j=1 where ~1 = -0.8611363116, 2 = -0.3399810436, =: 0.3399810436, 4 = 0.8611363116, TW = W4 = 0.3478548451, W2 = W3 = 0.6521451549, and r = x(gj)X + y(gj)y in which x( j) and y(gi) are given by (9). The treatment of the area integral over Q follows a similar procedure. We can again introduce the transformation 8 8 x - N~(, 71)xi, y = EN( (., 7)yi (14) i=l i=1 allowing the representation of an arbitrarily shaped quadrilateral with quadratically curved sides in the xy-plane to a square in the 67-plane as shown in Figure 3. The shape functions are now Nje (i = 1, 2,..., 8) and they take the known form [3] 1 1 Ne =- -(1 -- )(1- /)( + r + 1), N: = (1 + )(1-r -)(-ri/-1), 4 4 1 1 A"3- (+)( )( r-), N4 ++ 4 4 68

y TI 4 (0,1) 7 (-1,1) 4 (1,1 ) 3 2 8 (-1,0o >-6- --, (1,0) w H 0 1 5 2 (-1,-1 ) (0,-1 ) (1,-1 ) (b) (a) Figure 3: A quadrilateral element with quadratically curved sides in the xy-plane (a) can be transformed into a square in the local ~rl-plane (b). 69

1 1 N5 = (l1-2 (1-r ), = 2(1+ )(1- 2), -( 7e = ( - 2 (1 + r), N8e = (1- (1 - ). 2 2 With this transformation the area element ds can be expressed as ds = IJelddI (15) in which DOx Dy ox ay is the determinant of the Jacobian transformation matrix. Similarly to the one dimensional case, the shape function Ne(,r) is unity when ( = Hj and 71 = 7i. Thus, in accordance with the definition of the isoparametric elements we choose Ve = Nie and as such the coefficients <i will coincide with the nodal values of the field F6e(). When this expansion is substituted in (6) we obtain integrals of the form Ie = J A(r;')NT(r')ds' = J A(r l')N^(i', 1') (' ') (16) and by using a nine point Gaussian integration formula, PI can be written as 3 3 Ie = S W T4/jWkA(r-r'k)N;(Qj, rk3)lJ(j, rk)L, (17) j=l k=l where 1 = r1i = -0.7745966692, 22 = r2 = 0.0, (3 = r/3 = 0.7745966692, TlW = 14/3 = 0.5555555556, /T2 = 0.8888888889, and rk = xgj, 1k)X + y(, )y with x(j, 17k) and y(j, r1k) as given in (14). 70

We remark that the above formulation overcomes difficulties associated with singularities in the integrand. The numerical implementation is also rather straightforward and substantially simpler than those employed in traditional formulations. Most importantly the accuracy of the solution will be seen to be remarkable and to render accurate results for the internal fields for a scatterer having a large value of u. 4. Numerical Results To validate the above moment method formulation and show its usefulness, we performed a few numerical experiments. For the results shown below, the incident field is assumed to be a plane wave propagating in the x-direction and the axis of the cylinder is coincident with the z-axis. We first consider the problem of plane wave scattering by a dielectric circular cylinder with a radius 0.05A (free-space wavelength) and relative permittivity as high as or = 72.0-jl 62.0 corresponding to a muscle cylinder at 100 MHz. Figures 4 and 5 show, respectively, the electric and magnetic fields inside the cylinder for the TI\ and TE cases. Also, Table 1 gives the bistatic radar cross section, all compared with the exact eigenfunction solutions. As seen, there is an excellent agreement when the element size is chosen sufficiently small. We note that our previous moment method codes [1],[2] employing rectangular and triangular elements with pulse basis functions are unable to predict the correct TE result shown in Figure 5. A similar conclusion was also reached by Peterson and Klock [6] when they examined Richmond's 18] formulation where pulse basis functions were employed 71

Table 1. Bistatic Radar Cross Section (a/A) for a Homogeneous Dielectric Circular Cylinder with Radius Equal to 0.05A and Relative Permittivity 72.0 -- j162.0 TM case TE case Angle (deg.) 20 unk 45 unk Exact 45 unk 76 unk Exact 0 -4.03 -3.96 -3.95 -20.44 -20.46 -20.54 30 -4.19 -4.12 -4.11 -22.23 -22.49 -22.26 60 -4.63 -4.56 -4.55 -30.28 -28.93 -28.99 90 -5.25 -5.19 -5.18 -28.35 -26.92 -26.73 120 -5.89 -5.83 -5.82 -19.99 -19.72 -19.65 150 -6.38 -6.31 -6.31 -16.77 -16.64 -16.60 180 -6.56 -6.49 -6.49 -15.85 -15.75 -15.71 for the solution of the electric field integral equation. In [6], Peterson and Klock presented an improved magnetic field integral equation formulation by employing triangular elements with linear basis functions in the TE case. Here we compare our results with those in [6] for a dielectric cylinder with a radius 0.05A having a relative permittivity c. = 4.0 - jl100.0. The results are shown in Figure 6 and as seen the isoparametric elements provide a higher accuracy. In addition to the extreme situation described above, we also compared results obtained using isoparametric elements with those attributed to pulse expansion 72

0.25 0.20 exact \ 20 unknowns 0. 15 - 45 unknowns t X 0,10 - 0,.05 - 0.00, -0.050 -0.025 0.000 0.025 0.050 x/X (a) 0.25 0.20 0.15 - 0.10 / 0.05 0.00 - 0.000 0.025 0.050 y/X (b) Figure 4: The electric field along (a) the x-axis and (b) y-axis inside a homogeneous dielectric circular cylinder of radius a - 0.05A and cr -= 72.0 -j 162.0 for TM wave incidence. 73

1.5 o 76 unknowns 1. 0.5 -0.0 -0.050 -0.025 0.000 0.025 0.050 x/X (a) 1.5 1.0 zN 0.5 0.0 0.000 0.025 0.050 y/X (b) Figure 5: The magnetic field along (a) the x-axis and (b) /-axis inside a homogeneous dielectric circular cylinder of radius a = 0.05A and or = 72.0 - j162.0 for TE wave incidence. 74

1.40 exact (eigenfunction) 1.20 - E 45 unknowns (this method) o 76 unknowns (this method) 1.00 A 85 unknowns (Peterson & Klock)4 0.80 \ l/ 0.60 - o 0.40 -0.050 -0.025 0.000 0.025 0.050 x/X Figure 6: The magnetic field along the x-axis inside a homogeneous dielectric circular cylinder of radius a = 0.05A and cr = 4.0 - j100. for TE wave incidence. A comparison with exact data and those given in [6]. 75

functions for the case of nominal refractive indices. Overall, it was observed that the solution with isoparametric elements is more accurate and efficient than those employing pulse expansion basis. The improvement in these cases, though, is not as pronounced. 5. Conclusion In this paper, we examined three integral equations for TE scattering by dielectric cylinders having large values of permittivity. A moment method solution of the volume-surface integral equation was then developed by employing isoparametric elements and point-matching. Differing from the traditional solutions using pulse basis, the one presented here was shown to be more accurate and stable, particularly in the case of scatterers having large refractive indices. References [1] Jin, J. M., V. V. Liepa, and C. T. Tai, "A volume-surface integral equation for electromagnetic scattering by inhomogeneous cylinders," The Journal of Electromagnetic Waves and Applications, vol. 2, pp. 573-588, 1988. [2] Ricoy, M. A., S. Kilberg, and J. L. Volakis, "A simple set of integral equations for two-dimensional scattering with further reduction on unknowns," submitted to Proc. Inst. Elec. Eng., part H. [3] Zienkiewicz, 0. C., The Finite Element Method, 3rd ed. New York: McGrawHill, 1977. 76

[4] Graglia, R. D., "The use of parametric elements in the moment method solution of static and dynamic volume integral equations," IEEE Trans. Antennas Propagat., vol. AP-36, pp. 636-646, 1988. [5] Borup, D. T., D. NM. Sullivan, and 0. P. Gandhi, "Comparison of the FFT conjugate gradient method and the finite-difference time-domain method for the 2-D absorption problem," IEEE Trans. Microwave Theory Tech., vol. MTT-35, pp. 383-395, 1987. [6] Peterson, A. F. and P. W. Ilock, "An improved MFIE formulation for TEscattering from lossy, inhomogeneous dielectric cylinders," IEEE Trans. Antennas Propagat., vol. AP-36, pp. 45-49, 1988. [7] Boyse, W. E. and W. D. Kennedy, "A study of two-dimensional cylinders illuminated by TE polarized microwaves using the k-space algorithm," 1988 IEEE AP-S International Symposium Digest, vol. 1, pp. 84-87, Syracuse, NY, June 1988. [8] Richmond, J. H., "TE-wave scattering by a dielectric cylinder of arbitrary cross-section shape," IEEE Trans. Antennas Propagat., vol. AP-14, pp. 460-464, 1966. 77

DESCRIPTION AND LIST OF CODE VSIEM Jian-Ming Jin Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, Michigan 48109 ABSTRACT A FORTRAN program named as VSIEM and developed by the author is described briefly and listed thereafter. 1. Objective of VSIEM VSIEM is a computer code written in FORTRAN language and developed for computing electromagnetic scattering from perfectly conducting and inhomogeneous penetrable cylinders for both TE and TM incidence. VSIEM uses the pulse expansion functions and point-matching technique to solve the volume-surface integral equation (VSIE). For the formulation of VSIE, please read the first paper in this report. For the description of the numerical implementation, please read the second paper. VSIEM was used to generate the data presented in the two papers. 78

2. Input Data All input data are provided by file INPUT except for relative permittivity and permeability and their derivatives which are specified by subroutines EMFUN and DRFUN. A detailed description of the input data is given in the main program. A sample file of INPUT is listed below. COATED CIRCULAR CYLINDER 0 2. 1 0 0. 180. 30. 0. 10 -0.4 0.0 0.4 0.0 180. 10 0.4 0.0 -0.4 0.0 180. 0 10 -0.44 0.0 0.44 0.0 180. 10 0.44 0.0 -0.44 0.0 180. 0 10 -0.42 0.0 0.42 0.0 180. 0.04 10 0.42 0.0 -0.42 0.0 180. 0.04 0 3. Output Data The output data are written in file OUTPUT. The data contain the geometry and discretization information and the computed result. Following is a sample output. COATED CIRCULAR CYLINDER SEG NUM ENDPOINTS OF THE SEGMENT SEGMENT PARAMETERS NUM CELLS XA YA XB YB ANGLE RADIUS LENGTH 1 10 -0.40000 0.00000 0.40000 0.00000 180.00 0.400 1.2566 2 10 0.40000 0.00000 -0.40000 0.00000 180.00 0.400 1.2566 SEG NUM ENDPOINTS OF THE SEGMENT SEGMENT PARAMETERS NUM CELLS XA YA XB YB ANGLE RADIUS LENGTH 1 10 -0.44000 0.00000 0.44000 0.00000 180.00 0.440 1.3823 2 10 0.44000 0.00000 -0.44000 0.00000 180.00 0.440 1.3823 LAY NUM ENDPOINTS OF THE LAYER LAYER PARAMETERS NUM CELLS XA YA XB YB ANGLE RADIUS LENGTH 1 10 -0.42000 0.00000 0.42000 0.00000 180.00 0.420 1.3195 2 10 0.42000 0.00000 -0.42000 0.00000 180.00 0.420 1.3195 79

KEY PARAMETERS INCIDENT POLARIZATION NUMBER OF POINTS ON SURFACE(C) NUMBER OF POINTS ON SURFACE(D) NUMBER OF POINTS INSIDE DIELE NUMBER OF INCIDENT FIELD DIRECTIONS NUMBER OF BISTATIC DIRECTIONS WAVELENGTH E 20 20 20 1 7 2.00000 BISTATIC SCATTERING CROSS SECTION FOR INCIDENT FIELD DIRECTION= 0.00 THETA 10*LOG(SIGMA/LAMBDA) PHASE 0.00 -1.33 -93.0 30.00 -1.39 -98.1 60.00 -1.63 -114.1 90.00 -1.70 -143.6 120.00 0.15 -179.9 150.00 2.91 158.0 180.00 4.05 151.7 4. Program Description Following is a brief description of the main program and subroutines developed for VSIEM. All others not described here are standard subroutines. VSIEM - EMFUM - DRFUN - GEOMS - GEOMI - MATX1 - MATX2 - FARF1 - FARF2 - Main program acting as a driver for subroutines. Specifies cr and /r at a given point (x, y). Specifies Vcr and VUy at a given point (x, y). Discretizes the surfaces of conducting and dielectric cylinders. Discretizes the cross section of the dielectric cylinder. Generates the system matrix for E-polarization. Generates the system matrix for H-polarization. Computes the far field and RCS for E-polarization. Computes the far field and RCS for H-polarization. 80

C*********************** ***********************************************C PROGRAM VSIEM C*********************************************************************** C INPUT FORMAT FOR PROGRAM VSIEM ----VERSION OF JAN 21, 1988 C C MODIFIED JAN 21, 1988 C C THIS IS A GENERALIZED VERSION OF VSIEM EXTENDED TO TREAT CONDUCTING C C SURFACE. --- APRIL 23, 1988 C C*********************************** C CARD 1 FORMAT (18A4) TITLE CARD; USE UP TO 72 COLUMNS C C**********************************************************************C C CARD 2 FORMAT (I2,F10.5) KODE,WAVE C C KODE=O COMPUTES BISTATIC SCATTERING PATTERN C C KODE=1 COMPUTES BACKSCATTERING PATTERN C C WAVE WAVELENGTH C C CARD 3 FORMAT (I2,I3,4F10.5) IPP,IOPT,FIRST,LAST,INK,CANG C C IPP=I E-POLARIZATION C C IPP=2 H-POLARIZATION C C IOPT=O PARAMETERS NOT PRINTED C C IOPT=1 PARAMETERS PRINTED C C FIRST INITIAL SCATTERING AND INCIDENCE ANGLE C C LAST FINAL ANGLE C C INK ANGULAR INCREMENT C C CANG ANGLE FOR BISTATIC COMPUTATION C C*********************************************** ***********************C C FORMAT FOR INPUT CONDUCTING SURFACE C C**********************************************************************C C CARD 4 FORMAT (I5,5F10.5) N,XA,YA,XB,YB,ANGLE C C N NUMBER OF SAMPLING POINTS ON THIS SEGMENT C C XA,YA,XB,YB SEGMENT ENDPOINTS C C ANGLE ANGLE SUBTENDED BY THE SEGMENT C C***************************************************************-*******C C REPEAT CARD 4 FOR EACH SEGMENT C C CARD 5 FORMAT (12) INTEGER ZERO IN COLUMN 5; SHUTS C C OFF READING OF SEGMENT PARAMETERS C C FORMAT FOR INPUT DIELECTRIC SURFACE C C**********************************************************************C C CARD 6 FORMAT (I5,5F10.5) N,XA,YA,XB,YB,ANGLE C C N NUMBER OF SAMPLING POINTS ON THIS SEGMENT C C XA,YA,XB,YB SEGMENT ENDPOINTS C C ANGLE ANGLE SUBTENDED BY THE SEGMENT C C*********************** ***********************************************C C REPEAT CARD 6 FOR EACH SEGMENT C C*********************** *********************************************C C CARD 7 FORMAT (12) INTEGER ZERO IN COLUMN 5; SHUTS C C OFF READING OF SEGMENT PARAMETERS C C***********************************************************************C C FORMAT FOR INPUT DIELECTRIC LAYERS C C************************************************************************C C CARD 8 FORMAT (I2,6F:L0.5) N,XA,YA,XB,YB,ANGLE,DL C C N NUMBER OF SAMPLING POINTS ON THIS LAYER C C XA,YA,XB,YB LAYER ENDPOINTS C C ANGLE ANGLE SUBTENDED BY THE LAYER C C DL THICKNESS OF THE LAYER C C**********************************************************************C 81

C REPEAT CARD 8 FOR EACH LAYER C C CARD 9 FORMAT (12) INTEGER ZERO IN COLUMN 5; SHUTS C C OFF READING OF LAYER PARAMETERS C COMPLEX A(500,500),FINC(500) COMPLEX EPS(100),14U(100),EPSI (300),MUI(300) COMPLEX DXEPS(300),DYEPS(300),DXMU(300),DYMU(300) REAL LASTINK DIMENSION XC(100),YC(100),XNC(100),YNC(100),SC(100),DSQC(100), a ANGC(10oc) DIMENSION XS(100),YS(100),XN(100),YN(100),S(100),DSQ(100),ANG(100) DIMENSION XI(300),YI(300),SI(300),DSQI(300),ANGI(300) DIMENSION ID(18),LUMPC(100,2),LUMPS(100,2),LUMPI(300,2),IPOL(2) DIMENSION THEA(361),SCATA(361),FASEA(361),KPVT(500) DATA RED,DIG,PIF,IPOL/0.01745329,57.29578,0.07957747151, &4HEEEE,4HHHHH/ OPEN(UNIT=5,FILE 'INPUT') OPEN(UNIT=6,FILE='OUTPUT') C.....READ INPUT DATA AND GENERATE BODY PROFILE S READ (5,100) ID READ (5,125) KODE,WAVE READ (5,200) IPP,IOPT,FIRST,LAST,INK,CANG WRITE (6,150) ID WRITE (6,300) CALL GEOMS(LUMPC,XC,YC,XNC,YNC,SC,DSQC,ANGC,EPSMUMCLLC) IF(IOPT.EQ.0) GOTO 8 WRITE (6,500) DO 6 I=1,MC 6 WRITE(6,250)(LUMPC(I,J),J=1,2),XC(I),YC(I),XNC(I),YNC(I), &SC(I),DSQC(I) 8 WRITE (6,300) CALL GEOMS(LUMPS,XS,YS,XN,YN,S,DSQ,ANG,EPS,MU,MSLLS) IF(IOPT.EQ.0) GOTO 20 WRITE (6,500) DO 10 I=1,MS 10 WRITE(6,250)(LUMPS(I,J),J=1,2),XS(I),YS(I),XN(I),YN(I), &S(I),DSQ(I),EPS(I),MU(I) 20 WRITE (6,310) CALL GEOMI(LUMPI,X'I1,YI,SI,DSQI,ANGI,EPSI,MUI,DXEPS,DYEPS, &DXMU,DYMU,MI,LLI) IF(IOPT.EQ.0) GOTO 40 WRITE (6,510) DO 30 I=1,MI 30 WRITE(6,260)(LUMPI(I, J),J=1,2),XI(I),YI(I),S(I),DSQI(I), &EPSI(I),MUI(I) 40 CONTINUE C.....GENERATE SYSTEM MATRIX XK=6.283185/WAVE IF(IPP.EQ.2) GOTO 50 CALL MATX1(A,MC,MS.,MI,XC,YC,XS,YS,XNYN,XI,YIDSQCDSQDSQI, &EPS,MU,EPSI,MUI,DXMIU,DYMU,XK) GOTO 60 50 CALL MATX2(A,MC,MS,MI,XC,YC,XS,YS,XNC,YNC,XNYNXIYIDSQC, &DSQ,DSQI,MI,EPS,MUJE,EPSI,DXEPS,DYEPS,XK) 60 CONTINUE N=MC+MS+MI CALL CGECO(A,S00,NKPVT,RC,FINC) IF (KODE.NE.0) GO TO 70 82

NINC=l NBIT=1+IFIX((LAST-FIRST)/INK) GO TO 80 70 NBIT=O NINC=1+IFIX ((LAST-FIRST)/INK) 80 WRITE (6,400) IPOL(IPP),MC,MS,MI,NINC,NBIT,WAVE IF(KODE.EQ.1) WRITE.(6,600) IF(KODE.EQ.1) CANG:=FIRST C.....COMPUTE INCIDENT FIELD 82 TETA=RED*CANG CT=COS (TETA) ST=SIN (TETA) DO 84 I=1,MI HOLD=XK* (CT*XI (I)+ST*YI (I)) 84 FINC(I)=CMPLX(COS(HOLD),SIN(HOLD)) DO 85 I=1,MS HOLD=XK* (CT*XS (I)+ST*YS (I)) 85 FINC(I+MI)=CMPLX(COS(HOLD),SIN(HOLD)) DO 851 I=1,MC HOLD=XK* (CT*XC(I)+ST*YC(I)) 851 FINC(I+MI+MS)=CMPLX (COS(HOLD),SIN(HOLD)) C.....COMPUTE NEAR FIELD (NODAL FIELD VALUES) CALL CGESL(A,500,N,.KPVT,FINC,O) C DO 98 I=1,N C 98 WRITE(6,99) I,FINC(I),CABS(FINC(I)) C 99 FORMAT(5X,I5,1OX,;!E16.6,10X,F16.8) C.....CALL SUBROUTINE TO COMPUTE FAR FIELD IF(KODE.EQ.0) GOTO 88 IF(IPP.EQ.2) GOTO 86 CALL FARF1(FINC,MC,MS,MI,XC,YC,XS,YS,XN,YN,XI,YI,DSQC,DSQ,DSQI, IEPS,MU,EPSI,MUI,DXMU,DYMU,XK,CANG,CANG,1) GOTO 87 86 CALL FARF2(FINC,MC,MS,MI,XC,YC,XS,YS,XNC,YNC,XN,YN,XI,YI,DSQC, &DSQ,DSQI,MU,EPS,MUI,EPSI,DXEPS,DYEPS,XK,CANG,CANG,1) 87 CANG=CANG+INK IF(CANG.LE.LAST) GOTO 82 GOTO 1000 88 WRITE(6,610) CANG IF(IPP.EQ.2) GOTO 90 CALL FARF1(FINC,MC,MS,MI,XC,YC,XS,YS,XN,YN,XI,YI,DSQC,DSQ,DSQI, &EPS,MU,EPSI,MUI,DXMU,DYMU,XK,FIRST,INK,NBIT) GOTO 1000 90 CALL FARF2(FINC,MC,MS,MI,XC,YC,XS,YS,XNC,YNC,XN,YN,XI,YI,DSQC, &DSQ,DSQI,MU,EPS,MUIl,EPSI,DXEPS,DYEPS,XK,FIRST,INK,NBIT) 100 FORMAT (18A4) 125 FORMAT(I2,FlO.5) 150 FORMAT (1X,18A4) 200 FORMAT (I2,I3,4F10.5) 250 FORMAT (215,10F10.5) 260 FORMAT (215, 8F10.5) 300 FORMAT (/,/lOl SEG NUM,11X,24HENDPOINTS OF THE SEGMENT,11X, &' SEGMENT PARAMETERS '/11H NUN CELLS,6X,2HXA,8X,2HYA,8X,2HXB,8X, &2HYB,6X,21HANGLE RADIUS LENGTH) 310 FORMAT (/,/10H LAY NUM,11X,24HENDPOINTS OF THE LAYER,11X, &' LAYER PARAMETERS '/11H NUM CELLS,6X,2HXA,8X,2HYA,8X,2HXB,8X, &2HYB,6X,21HANGLE RADIUS LENGTH) 400 FORMAT (//25X,14HKEY PARAMETERS// &1OX,21HINCIDENT POLARIZATION,22X,1A1/ &10X,34HNUMBER OF POINTS ON SURFACE(C),I10/ 83

&10X,34HNUMBER OF POINTS ON SURFACE(D),I10/ &10X,34HNUMBER OF POINTS INSIDE DIELE,I10/ &10X,35HNUMBER OF INCIDENT FIELD DIRECTIONS,I9/ &10X,29HNUMBER OF BISTATIC DIRECTIONS,I15/ &1OX,10HWAVELENGTH,F34.5) 500 FORMAT (/,/11H I SEG,4X,4HX(I),6X,4HY(I),5X,5HXN(I), &5X,5HYN(I),6X,4HS(I),5X,6HDSQ(I),3X,7HEPSR(I),3X,7HEPSI(I), &4X,6HMUR(I),4X,6HMUI(I)) 510 FORMAT (/,/11H I SEG,4X,4HX(I),6X,4HY(I),6X,4HS(I),4X, &7HAREA(I),3X,7HEPSR(I),3X,7HEPSI(I),4X,6HMUR(I),4X,6HMI(I)) 600 FORMAT (///,20X,28HBACKSCATTERING CROSS SECTION/17X, &36HTHETA 10*LOG(SIGMA/LAMBDA) PHASE/) 610 FORMAT (///,19X,33HBISTATIC SCATTERING CROSS SECTION/18X, &29HFOR INCIDENT FIELD DIRECTION=,F7.2/17X, &36HTHETA 10*LOG(SIGMA/LAMBDA) PHASE/) 700 FORMAT(///,18X,35HFIELD DISTRIBUTION INSIDE SCATTERER/15X, &53HPONIT NUM COMPLEX FIELD MAG(FIELD)/) 710 FORMAT(15X,I6,5X,E11.4,4X,E11.4,SX,F9.4) 1000 STOP END 84

SUBROUTINE EMFUN CX,Y, EPS,MU) COMPLEX EPS,MU C.....SPECIFY PERMITTIVITY AND PERMEABILITY EPS=(2.4,0.0) MU=C1. 2, 0. 0) cc EPS=2.-(X*X+Y*Y)/(.4*.4) cc EPS=4./(l+(X*X+YiY)/(O.36))**2 RETURN END SUBROUTINE DRFUN(X,Y,DXEPS,DYEPS,DXMU,DYMU) COMPLEX DXEPS,DYEPS, DXMU, DYMU C.....SPECIFY THE DERIVATIVES OF PERMITTIVITY AND PERMEABILITY DXEPS=(O.,O.) DYEPS=(O.,O.) cc DXEPS=-2.*X/(0.4,0.4) cc DYEPS=-2.*Y/(0.4,0.4) cc DXEPS=-16./((l+(Xi(*X+Y*Y)/(0.36))**3)*X/0.36 cc DYEPS=-16./((l+(X]*X+Y*Y)/(0.36))**3)*Y/0.36 DXMU= (0. 0) DYMU'(.,o.) RETURN END 85

c** * *********************************************************************c C TTHIS SUBROUTINE READS DATA AND GENERATE POINTS ON THE SURFACE C C OF THE CYLINDER C SUBROUTINE GEOMS(LUMP,X,Y,XN,YN,S,DSQ,ANG,EPS,MU,M,LL) COMPLEX EPS(100),MU(100) DIMENSION X(100),Y(100),XN(100),YN(100),DSQ(100),S(100),ANG(100) DIMENSION LUMP(100, 2) DATA RED/0.01745329/ I=O L=O C.....READ INPUT PARAMETERS AND PREPARE TO GENERATE SAMPLING POINTS 10 READ (5,200) N,XA,YA,XB,YB,ANGLE IF (N.EQ.O) GO TO 120 LIM=2*N-1 TX=XB-XA TY=YB-YA D=SQRT(TX*TX+TY*TY) IF (ANGLE.EQ.0.0) GO TO 20 T=0.5*RED*ANGLE TRX=TX+TY/TAN(T) TRY=TY-TX/TAN(T) RAD=0.5*D/SIN(T) ARC=2.0*RAD*T ALF=T/N DID=2.0*RAD*ALF GO TO 30 20 RAD=999.999 ARC=D DID=D/N C.....START GENERATING 30 CONTINUE L=L+i DO 100 J=1,LIM,2 1=1+1 LUMP(I,1)=I LUMP(I,2)=L IF (ANGLE.EQ.0.0) GO TO 40 SINQ=SIN(J*ALF) COSQ=COS (J*ALF) X(I)=XA+0.5*(TRX*(1.0-COSQ)-TRY*SINQ) Y(I)=YA+0.5*(TRX*SINQ+TRY*(1.0-COSQ)) XN(I)=-0.5*(TRX*CO:SQ+TRY*SINQ)/RAD YN(I)= 0.5*(TRX*SINQ-TRY*COSQ)/RAD ANG (I) =(ANGLE/N) *RED GO TO 50 40 X(I)=XA+0.5*J*TX/N Y(I)=YA+0.5*J*TY/N XN(I)=-TY/D YN(I)= TX/D ANG(I)=(ANGLE/N) *RED 50 S(I)=0.5*J*DID C.....COMPUTE THE PERMITTIVITY AND PERMEABILITY CALL EMFUN(X(I),Y(]:),EPS(I),MU(I)) 100 DSQ(I)=DID WRITE (6,300) L,N,XA,YA,XB,YB,ANGLE,RAD,ARC GO TO 10 120 M=I LL=L 86

200 FORMAT (IS,5F10.5) 250 FORMAT (5X,5F10.5) 300 FORMAT (I3,I6,3X,4F10.5,F8.2,F8.3,F8.4) RETURN END 87

C TTHIS SUBROUTINE READS DATA AND GENERATE POINTS INSIDE THE C C CYLINDER C C********************** ***********************************************C SUBROUTINE GEOMI(LUMP,X,Y,S,DSQ,ANG,EPS,MU,DXEPS,DYEPS, &DXMU,DYMU,M,LL) COMPLEX EPS(300),MU(300),DXEPS(300),DYEPS(300),DXMU(300),I)YMU(300) DIMENSION X(300),Y(300),DSQ(300),S(100),ANG(100) DIMENSION LUMP(300,2) DATA RED/0.01745329/ I=O L=O C.....READ INPUT PARAMETERS AND PREPARE TO GENERATE SAMPLING POINTS 10 READ (5,200) N,XA,YA,XB,YB,ANGLE,THICK IF (N.EQ.O) GO TO 120 LIM=2*N-1 TX=XB-XA TY=YB-YA D=SQRT(TX*TX+TY*TY) IF (ANGLE.EQ.O.O) GO TO 20 T=0.5*RED*ANGLE TRX=TX+TY/TAN(T) TRY=TY-TX/TAN(T) RAD=0.5*D/SIN(T) ARC=2.0*RAD*T ALF=T/N DID=2.0*RAD*ALF GO TO 30 20 RAD=999.999 ARC=D DID=D/N C.....START GENERATING 30 CONTINUE L=L+I DO 100 J=1,LIM,2 I=I+1 LUMP(I,1)=I LUMP(I,2)=L IF (ANGLE.EQ.O.O) GO TO 40 SINQ=SIN(J*ALF) COSQ=COS(J*ALF) X(I)=XA+. 5*(TRX*(1.0-COSQ)-TRY*SINQ) Y(I)=YA+0.5*(TRX*SINQ+TRY*(1.O-COSQ)) ANG(I)=(ANGLE/N)*RED GO TO 50 40 X(I)=XA+0.5*J*TX/N Y(I)=YA+0.5*J*TY/N ANG (I)=(ANGLE/N)*RED 50 S(I)=0.5*J*DID C.....COMPUTE THE PERMITTIVITY AND PERMEABILITY CALL EMFUN(X(I),Y(I),EPS(I),MU(I)) CALL DRFUN(X(I),Y(I),DXEPS(I),DYEPS(I),DXMU(I),DYMU(I)) 100 DSQ(I)=DID*THICK WRITE (6,300) L,N,XA,YA,XB,YB,ANGLE,RAD,ARC GO TO 10 120 M=I LL=L 200 FORMAT (I5,6F10.5) 250 FORMAT (5X,5F10.5) 88

300 FORMAT (I3,I6,3X,4F10.5,F8.2,F8.3,F8.4) RETURN END 89

C SUBROUTINE FORMS SYSTEM MATRIX [A] FOR E-POLARIZATION; C C FOR H-POLARIZATION,, USE MATX2. C SUBROUTINE MATX1 (AMC,MS,MI,XC,YC,XS,YS,XN,YN,XI,YI,DSQC,DSQ,DSQI, &EPS,MU,EPSI,MUI,DXMU,DYMU,XK) COMPLEX A(500,500),,EPS(100),MU(100),EPSI(300),MLJI(300),DXMIJ(300), &DYMU(300),HZERO0,Ho0rE DIMENSION XC(100),Y'C(100),XS(100),YS(100),XN(100),YN(100),XEI(300), &YI(300),DSQC(100),D)SQ(100),DSQI(300) PI=3. 141592654 DO 20 I=1,MI DO 20 J=1,MI IF(I.EQ.J) GOTO 10 XJIXI (J)-XI(I) YJI=YI CJ)-YI(I) RJI=SQRT (XJI*XJI+YJ'I*Y.JI) RKO=XK*RJI CALL HANKZ2(RKO,2,H'ZERO,HONE) A(I,J)=(XK*XK*(EPSI'(J)-il./MUI(J))*HZERO+(DXMU(J)*XJI/RJI+ &DYMU(J)*YJI/RJI)/(!'~UI(J)*MUI(J))*XK*HONE)*CMPLX(0.,0.25)*DSQI(J) GOTO 20 10 RKA=SQRT(DSQI(J)/PI)*XK CALL HANKZ2 (RKA, 1, H ZERO,HONE) A(I,I)=l./MUI(I)+0.5*(EPSI(I)-1l./MUICI))*(2.+CMPLX(0. 11.) &*PI*RKA*HONE) 20 CONTINUE DO 30 I=1,MI DO 30 J=1,MS XJI=XS (J)-XI(I) YJI=YS (J)-YI(I) RJI=SQRT (XJI*XJI+YJI*YJI) RKO=XK*RJI CALL HANKZ2 (RKO, 1,H'ZERO DRONE) 30 A(I,J+MI)=(1l./MU(J)-il.)*(XJI/RJI*XN(J)+YJI/RJI*YN(J)) &*CMPLX(0.,0.25)*XK*HONE*DSQ(J) DO 40 I=1,MS DO 40 J=1,MI XJI=XI CJ)-XS(I) YJIYI (J)-YS(I) RJI=SQRT (XJI*XJI+YJI*YJI) RKO=XK*RJI CALL HANKZ2(RKO,2,HZERO,HONE) 40 A(I+MI,J)=(XK*XK*(EPSI(J)-l./MUI(J))*HZERO+(DXMU(J)*XJI/RJI-+ &DYMU(J)*YJI/RJI)/(MUI(J)*MUI(J))*XK*HONE)*CMPLX(0.,0.25)*DSQI(J) DO 60 I=1,MS DO 60 J=1,MS IF(I.EQ.J) GOTO 50 XJI=XS (J)-XS(I) YJI=Y5(J)-YS(I) RJI=SQRT (XJI*XJI+YJI*YJI) RKO=XK*RJI CALL HANKZ2(, 1, HZERO,HONE) A(I+MI,J+MI)=(l./MU(J)-1.)*(XJI/RJI*XN(J)+YJI/RJI*YN(J)) &*CMPLX(0.,0.25)*XK*HONE*DSQ(J) GOTO 60 50 A(I+MI,J+MI)=0.5*(l.+l./MU(I)) 60 CONTINUE DO 70 I=1,MI 90

DO 70 J=1,MC XcJI=XC(J)-XI(I) YJI=YC(J)-YI(I) RJI=SQRT CXJI*XJI+YJr1*YJI) RKO=XK*RJI CALL HANKZ2 CRKO, 0, HZERO,HONE) 70 A(I,J+MI+MS)=-CMPL:(0.,0.25)*HZERO*DSQC(J) DO 80 I=1,MS DO 80 J=1,MC XJI=XC(J)-XSCI) YJI=YC(J)-YS(I) RJI=SQRT CXJI*XJI+YJrI*YJI) RKO=XK*RJI CALL HANKZ2 CRKO, 0, HZERO,HONE) 80 A(I+MI,J+MI+MS)=-CM[PLX(0.,0.25)*HZERO*DSQC(J) DO 90 I=1,MC DO 90 J=1,MI XJI=XI (J)-XC(I) YJI=YI (J)-YC(I) RJI=SQRT (XJI*XJI+YJ'I*YJI) RKO=XK*RJI CALL HANKZ2 CRKO, 2,HZERO,HONE) 90 A(I+MI+MS,J)=(XK*XK.*(EPSI(J)-l./MUJI(J))*HZERO+CDXNUCJ)*XJI/RJI+ &DYMU(J)*YJI/RJI)/(MfUI(J)*MUI(J))*XK*HONE)*CMPLX(0.,0.25)*DSQI(J) DO 100 I=1,NC DO 100 J=1,MS XJI=XS (J)-XC(I) YJI=YS (J)-YC(I) RJI=SQRT (XJI*XJI+YJ'I*YJI) RKO=XK*RJI CALL HANKZ2(RKO,1,EHZERO,HONE) 100 A(I+MI+MS,J+MI)=(l./MU(J)-l.)*(XJI/RJI*XN(J)+YJI/RJI*YN(J):) &*CMPLX(0.,0.25)*XK*:HONE*DSQ(J) DO 120 I=1,MC DO 120 J=1,MC IF(I.EQ.J) GOTO 110 XJI=XC(J)-XC(I) YJI=YC(J)-YCCI) RJI=SQRT (XJI*XJI+Y.J'I*YJI) RKO=XK*RJI CALL HANKZ2(RKO,0,B'ZERO,HONE) A(I+MI+MS,J+MI+MS)=-CMPLX(0.,0.25)*HZERO*DSQCCJ) GOTO 120 110 A(I+MI+MS,J+MI+MS)=-CMPLX(0.,0.25)*(1.0-CMPLX(0.,2.)/PI* &ALOG(0.1638*XK*DSQC (3))) *DSQC(J) 120 CONTINUE RETURN END 91

C SUBROUTINE FORMS SYSTEM MATRIX [A] FOR H-POLARIZATION; C SUBROUTINE MATX2(A.MC,MS,MI,XC,YC,XS,YS,XNC,YNCIN,YN,XI,YI,DSQC, &DSQ,DSQI,EPS,MU,EPS;I,MUI,DXMUJ,DYMU,XK) COMPLEX A(500,500),EPS(100),MU(100),EPSI(300),MUI(300),DXMU(300), &DYMUC300),HZERO,HOIIE DIMENSION XC(100),Y(C(100),XS(100),YS(100),XN(100),YN(100), &XNC(1oo),YNC(100),XI(300),YI(300),DSQC(100),DSQ(100),DSQI(.300) PI=3. 141592654 DO 20 I=1,MI DO 20 J=1,MI IF(I.EQ.J) GOTO 10 XJI=XI ()-XI(I) YJIYI (J)-YICI) RJI=SQRT (XJI*XJI+YJI*YJI) RKO=XK*RJI CALL HANKZ2 CRKO, 2, I[ZERO,HONE) A(I,J)=(XK*XK*CEPS]:-(J)-l./MUI(J))*HZERO+(DXMU(J)*XJI/RJI+ &DYMU(J)*YJI/RJI)/(MIUI(J)*MUI(J))*XK*HONE)*CMPLX(0.,0.25)*DSQI(J) GOTO 20 10 RKA=SQRTCDSQI(J)/PI')*XK CALL HANKZ2 CRKA, 1,1HZERO,HONE) A(I,I)=l./MUI(I)+0.5*(EPSI(I)-l./MIJI(I))*(2.+CMPLX(0.,l.) &*PI *RKA*HONE) 20 CONTINUE DO 30 I=1,MI DO 30 J=1,MS XJI=XS (J)-XI(I) YJI=YS (J)-YI(I) RJI=SQRT CXJI*X.JI+YJ'I*Y.JI) RKO=XK*RJI CALL HANKZ2(RKO,1,H'ZERO,HONE) 30 A(I,J+MI)=(1./MU(J)-1l.)*(XJI/RJI*XN(J)+YJI/RJI*YN(J)) &*CMPLX(0.,0.25)*XK*HONE*DSQCJ) DO 40 I=1,MS DO 40 J=1,MI XJI=XI CJ)-XS(I) YJI=YI (J)-YS(I) RJI=SQRT (XJI*XJI+YJI*YJI) RKO=XK*RJI CALL HANKZ2 CRKO,2, HZERO,HONE) 40 A(I+MI,J)=(XK*XK*(EPSI(J)-l./MUI(J))*HZERO+(DXMU(J)*XJI/RJI-+ &DYMU(J)*YJI/RJI)/CMUI(J)*MUI(J))*XK*HONE)*CMPLX(0.,0.25)*DSQI(J) DO 60 I=1,MS DO 60 J=1,MS IF(I.EQ.J) GOTO 5O XJI=XS (J)-XS(I) YJI=YS (J)-YS(I) RJI=SQRT (XJI*XJI+YJI*YJI) RKO=XK*RJI CALL HANKZ2(RKO,1,HZERO,HONE) A(I+MI,J+MI)=(l./MU(J)-l.)*(XJI/RJI*XN(J)+YJI/RJI*YN(J)) &*CMPLX(0.,0.25)*XK*HON~E*DSQ(J) GOTO 60 50 A(I+MI,J+MI)=0.5*(l.+l./MUJ(I)) 60 CONTINUE DO 70 I=1,MI DO 70 J=1,MC 92

XJI=XC(J)-XI(I) YJI=YC(J)-YI(I) RJI=SQRT (XJI*XJI+Y.JI*YJI) RKO=XK*RJI CALL HANKZ2 (RKO,1,]HZERO,HONE) 70 A(I,J+MI+NS)=-(XJI/tRJI*XNC(J)+YJI/RJT*YNC(J)) &*CMPLXCO.,0.25)*XK*~HONE*DSQCCJ) DO 80 I=1,MS DO 80 J=1,MC XJI=XCC.J)-XS(I) YJI=YC(J)-YS(i) RJI=SQRT (XJI*XJI+YJ'I*YJI) RKO=XK*RJI CALL HANKZ2 (RKO, 1, TIZERO,HONE) 80 A(I+MI,J+MI+MS)=-(XJI/RJI*XNC(J)+YJI/RJI*YNC(J)) &*CMPLXCO.,0.25)*XK*-HONE*DSQC(J) DO 90 I=1,MC DO 90 J=1,MI 13=X1 (3)-lCdI) YJi=Yi (3)-YC(I) RJI=SQRT (XJI*XJI.Y-Y'I*Y.JI) RKO=XK*RJI CALL HANKZ2 (RKO,2,HZERO,HOE 90 A(I+MI+MS,J)=(XK*XK:*(EPSI(J)-1./MUI(J))*HZERO+(DXMU(J)*XJI/IRJI+ &DYM.U(J)*YJI/RJI)/(P[UI(J)*NUI(J))*XK*HONE)*CIIPLX(0.,0.25)*DSQI(J) DO 100 I=1,NC DO 100 J=1,MS XJI=XS (J)-XCCI) YJI=YS (J)-YC(I) RJI=SQRT (XJI*XJI+YJ'I*YJI) RKO=XK*RJI CALL HANKZ2 (RKO, 1,EHZERO, HONE) 100 A(I+MI+MS,J+MI)=(l./M4U(J)-l.)*(XJI/RJI*XN(J)+YJI/RJI*YN(J):) &*CMPLX(0.,0.25)*XK*:HOEE*DSQ(J) DO 120 I=1,MC DO 120 J=1,MC IF(I.EQ.J) GOTO tic' XJI=XC(J)-XC(I) YJI=YC&J)-YC(I) RJI=SQRT (XJI*XJI+YJ I*YJI) RKO=XK*RJI CALL HANKZ2 (RKO, 1, H ZERO,HONE) A(I+MI+MS,J+MI+MS)=-(XJI/RJI*XNC(J)+YJI/RJI*YNC(J)) &*CMPLX(0.,0.25)*XK*HONE*DSQC(J) GOTO 120 110 ACI+MI+MS,J+MI+MS)=:0.5 120 CONTINUE RETURN END 93

C THIS SUBROUTINE CALCULATE FAR FIELD AND RADAR CROSS SECTION C SUBROUTINE FARF1(FINC,MC,MS,MI,XC,YC,XS,YS,XN,YN,XI,YI,DSQC,DSQ, &DSQI,EPS,MU,EPSI,MII,DXMU,DYMU,XK,FIRST,INK,NBIT) COMPLEX FINC(500),EPS(MS),MU(MS),EPSI(MI),MUI(MI) COMPLEX PSCAT,HZERO,HONE,DXMU(MI),DYMU(MI) REAL INK DIMENSION XC(MC),YC(MC),XS(MS),YS(MS),XN(MS),YN(MS),XI(MI),YI(MI), &DSQC(MC),DSQ(MS),DSQI (MI) PI=3.141592654 RED=O. 01745329 CANG=FIRST DO 100 I=1,NBIT PSCAT=CMPLX (0.,O.) TETA=RED*CANG COT=COS (TETA) SIT=SIN(TETA) DO 20 J=1,MI RP=SQRT(XI (3) *XI(J)+YI (J)*YI(J)) COP=XI (J)/RP SIP=YI (J)/RP HOLD=XK*RP* (COT*COP+SIT*SIP) HZERO=CMPLX(COS(HOl.D),SIN (HOLD)) HONE=CMPLX(0.,1. )*HZERO PSCAT=PSCAT-(XK*XK* (EPSI(J)-1./MUI(J))*HZERO-(DXMU(3)*COT+ &DYMU(J)*SIT)/(MUI( 3)*MUI(J))*XK*HONE)*CMPLX(0.,0.25)*DSQI((J) &*FINC(3) 20 CONTINUE DO 30 J=1,MS RP=SQRT(XS(3)*XS(J)+YS(3)*YS(J)) COP-XS (J)/RP SIP=YS (J)/RP HOLD=XK*RP* (COT*COP+SIT*SIP) HZERO=CMPLX(COS(HOLD), SIN (HOLD)) HONE=CMPLX (0.,1. )*EZERO 30 PSCAT=PSCAT+(1./MU(J)-1.)*(COT*XN(J)+SIT*YN(J)) &*CMPLX(0.,0.25)*XKHONE*DSQ(J)*FINC(MI+3) DO 40 3=1,MC RP=SQRT(XC(J)*XC(J)+YC(J)*YC(J)) COP=XC (3) /RP SIP=YC(3)/RP HOLD=XK*RP* (COT*COP+SIT*SIP) HZERO=CMPLX(COS(HOLD),SIN(HOLD)) 40 PSCAT=PSCAT+CMPLX(0.,0.25)*HZERO*DSQC(J)*FINC(MI+MS+3) PMAG=CABS (PSCAT) PPHS=ATAN2(AIMAG(PSCAT),REAL(PSCAT))*180./PI SIGMA=10. *LOG1O(2. *PMAG*PMAG/PI) WRITE(6,50) CANG,SIGMA,PPHS 50 FORMAT (9X,F13.2,F15.2,F16.1) CANG=CANG+INK 100 CONTINUE RETURN END 94

c*** **** ******* ************************** ** ********************** *******c C THIS SUBROUTINE CALCULATE FAR FIELD AND RADAR CROSS SECTION C SUBROUTINE FARF2(FINC,MC,MS,MI,XC,YC,XS,YS,XNC,YNC,XN,YN,XI,YI, &DSQC,DSQ, DSQI, EPS,MU,EPSI,MUI,DXMU,DYMU,XK,FIRST,INK,NBIT) COMPLEX FINC(500),EPS(MS),MU(MS),EPSI(MI),MUI(MI) COMPLEX PSCAT,HZERO,HONE,DXMU(MI),DYMU(MI) REAL INK DIMENSION XC(MC),YC(MC),XS(MS),YS(MS),XN(MS),YN(MS),XI(MI),YI(MI), &DSQC(Q(C),DSQ(MS),DS;QI(MI),XNC(M)),YNC(MC) PI=3.141592654 RED=0.01745329 CANG=FIRST DO 100 I=1,NBIT PSCAT=CMPLX(O.,0.) TETA=RED*CANG COT=COS (TETA) SIT=SIN(TETA) DO 20 J=1,MI RP=SQRT(XI(J)*XI (J) +YI(J)*YI(J)) COP=XI(J)/RP SIP=YI (J)/RP HOLD=XK*RP*(COT*COP+SIT*SIP) HZERO=CMPLX (COS(HOLD),SIN(HOLD)) HONE=CMPLX(0.,1. )*HZERO PSCAT=PSCAT-(XK*XK*(EPSI(J)-1./MUI(J))*HZERO-(DXMU(J)*COT+ &DYMU(J)*SIT)/(MUI (J)*MUI(J))*XK*HONE)*CMPLX(0.,0.25)*DSQI(J) &*FINC(J) 20 CONTINUE DO 30 J=1,MS RPSQRT(XS (J)*XS(J)+YS(J)*YS(J)) COP=XS (J)/RP SIP=YS(J)/RP HOLD=XK*RP*(COT*COP+SIT*SIP) HZERO=CMPLX(COS(HOLD), SIN(HOLD)) HONE=CMPLX(0.,1. )*HZERO 30 PSCAT=PSCAT+(1./MU(J)-1.)*(COT*XN(J)+SIT*YN(J)) &*CMPLX(0.,0.25)*XK*HONE*DSQ(J)*FINC(MI+J) DO 40 J=1,MC RP=SQRT(XC(J)*XC(J)+YC(J)*YC(J)) COP=XC(J)/RP $IP=YC(J)/RP HOLD=XK*RP* (COT*COP+SIT*SIP) HZERO=CMPLX(COS(HOLD), SIN(HOLD)) HONE=CMPLX(0.,1. )*HZERO 40 PSCAT=PSCAT-(COT*XNC(J)+SIT*YNC(J)) &*CMPLX(0.,0.25)*XK*HONE*DSQC(J)*FINC(MI+MS+J) PMAG=CABS(PSCAT) PPHS=ATAN2(AIMAG(PSCAT),REAL(PSCAT))*180./PI SIGMA=10. *LOG10(2. *PMAG*PMAG/PI) WRITE(6,50) CANG,SIGMA,PPHS 50 FORMAT (9X,F13.2,F15.2,F16.1) CANG=CANG+INK 100 CONTINUE RETURN END 95

C C SUBROUTINE HANKZ2CR,N,HZERO,HONE) C C C C C Called by subroutines IATXE and MATXH to compute Hankle C C functions of the second kind for orders one and zero. The C C argument is variable R and must be positive. C C C C.....HANKEL FUNCTIONS ARE OF FIRST KIND —J-iY C......1=O RETURNS HZERO (H-zero) C......1=1 RETURNS HONE (H-one) C......1=2 RETURNS HZERO AND HONE C.....SUBROUTINE REQUIRES R>0 C.....SUBROUTINE ADAM MUST BE SUPPLIED BY USER DIMENSION AT(7),BC7),,CC7),DC7),E(7),FC7),G(7),H(7) COMPLEX HZERO,HONE DATA A,B,C,D,E,F,G,H/1.0,-2.2499997,1.2656208,-0.3163866, &0.0444479,-0.0039444,0.00021,0.36746691,0.605S9366,-0.74350384, &0.25300117,-0.04261214,0.00427916,-0.00024846,0.5,-0.56249985, &0.21093573,-0.03954289,0.00443319,-0.00031761,0.00001109, &-0.6366198,0.2212091,2.1682709,-1.3164827,0.3123951,-0.0400976, &0.0027873,0.79788456,-0.00000077,-0.0055274,-0.00009512, &0.00137237,-0.00072805,0.00014476,-0.78539816,-0.04166397, &-0.00003954,0.00262573,-0.00054125,-0.00029333,0.00013558, &0.79788456,0. 00000156,0.01659661,0.00017105,-0.00249511, &0.00113653,-0.00020033,-2.35619449,0.12499612,0.0000565, &-0.00637879,0.00074348,0.00079824,-0.00029166/ IF CR.LE.O.0) GO TO 50 IF (N.LT.O.OR.N.GT.2) GO TO 50 IF CR.GT.3.0) GO TO 20 X=R*R/9.0 IF CN.EQ.1) GO TO 10 CALL ADAM(A,X,BJ) CALL ADAM(B,I,Y) BYO0.6366198*ALOG(O.5*R) *BJ+Y HZERO=CMPLX(BJ,-BY) IF (N.EQ.o) RETURN 10 CALL ADAMCC,X,Y) BJ=R*Y CALL ADAMCD,X,Y) BYO0.6366198*ALOGO.*R) *BJ+Y/R HONE=CMPLX(BJ,-BY) RETURN 20 X=3.0/R IF (N.EQ.1) GO TO 30 CALL ADAMCE,X,Y) FOOL=Y/SQRTCR) CALL ADAMCF,X,Y) T=R+Y BJ=FOOL*COSCT) BY=FOOL*SINCT) HZERO=CMPLXCBJ,-BY) IF CN.EQ.o) RETURN 30 CALL ADAMCG,X,Y) FOOL=Y/SQRTCR) CALL ADAMCH,X,Y) T=R+Y 96

BJ=FOOL*COS(T) BY=FOOL*SIN(T) HONE=CMPLX(BJ,-BY) RETURN 50 WRITE(6,90) N,R 90 FORMAT(32HOSICK DATA IN HANKZ2 *QUIT* N=,I2,2X,2HR=,E11.3) CALL SYSTEM END C$******* *** *******!****** ** ***** ** **** ************* ** ***** ******C C C SUBROUTINE ADAM(C,X,Y) C C C** *************** ************ * *************************************C C C C Called by subroutine HANKZ2 to compute the value of a 7th C C order polynomial whose argument is X and coefficients are C C contained in vector C. C DIMENSION C(7) Y=X*C(7) DO 10 I=1,5 Y=X*(C(7-I)+Y) 10 CONTINUE Y=Y+C(1) RETURN END 97

C** * ** ******** ********** * ********* *** ****** ******** ***** ***** ***** *C C C SUBROUTINE CGEFA(A,,LDA,N,IPVT,INFO) C C c* * ** ******** *** *******4 **** * ************************* ****** ***** *C C C C NAASA 2.1.043 CGEFA FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C INTEGER LDA,N,IPVT(C1),INFO COMPLEX A(LDA,1) C C CGEFA FACTORS A COMPLEX MATRIX BY GAUSSIAN ELIMINATION. C C CGEFA IS USUALLY CALLED BY CGECO, BUT IT CAN BE CALLED C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. C (TIME FOR CGECO) = (1 + 9/N)*(TIME FOR CGEFA) C C ON ENTRY C C A COMPLEX(LDA, N) C THE MATRIX TO BE FACTORED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A C C N INTEGER C THE ORDER OF THE MATRIX A C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS C WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C INFO INTEGER C = 0 NORMAL VALUE. C = K IF U(K,K).EQ. 0.0. THIS IS NOT AN ERROR C CONDITION FOR THIS SUBROUTINE, BUT IT DOES C INDICATE THAT CGESL OR CGEDI WILL DIVIDE BY ZERO C IF CALLED. USE RCOND IN CGECO FOR A RELIABLE C INDICATION OF SINGULARITY. C C LINPACK. THIS VERSION DATED 07/14/77 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. C C SUBROUTINES AND FUNCTIONS C C BLAS CAXPY,CSCAL,ICAMAX C FORTRAN ABS,AIMAG,CMPLX,REAL C C INTERNAL VARIABLES C COMPLEX T INTEGER ICAMAX,J,K,:KP1,L,NM1 98

c COMPLEX ZDUM REAL CABS1 CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM)) C CCC Gaussian elimination with partial pivoting C INFO = 0 NM1 = N - 1 IF (NM1.LT. 1) GO TO 70 DO 60 K = 1, NM1 KP1 = K + 1 C C FIND L = PIVOT INDEX C L = ICAMAX(N-K+1,A(K,K),l) + K - 1 IPVT(K) = L C CCC Zero pivot implies this column already triangularized C IF (CABS1(A(L,K)).EQ. O.OEO) GO TO 40 C CCC Interchange if necessary C IF (L.EQ. K) GO TO 10 T = A(L,K) A(L,K) = A(K,K) A(K,K) = T 10 CONTINUE C CCC Compute multipliers C T = -CMPLX(1.OEO,O.OEO)/A(K,K) CALL CSCAL(N-K,T,A(K+1,K),1) C CCC Row elimination with column indexing C DO 30 J = KP1, N T = A(L,J) IF (L.EQ. K) GO TO 20 A(L,J) = A(K,J) A(K,J) = T 20 CONTINUE CALL CAXPY(N-K,T,A(K+1,K),l,A(K+1,J),1) 30 CONTINUE GO TO 50 40 CONTINUE INFO = K 50 CONTINUE 60 CONTINUE 70 CONTINUE IPVT(N) = N IF (CABS1(A(N,N)).EQ. O.OEO) INFO = N RETURN END C C C*********************** ****************************e*****e****C C C SUBROUTINE CGESL(A,]LDA,N,IPVT,B,JOB) C C 99

c** ***** * ***************4'************* ******* ********** * ***** ***** *c c C C NAASA 2.1.044 CGESL FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C INTEGER LDA,N,IPVT(1),JOB COMPLEX A(LDA,1),B(:) C C CGESL SOLVES THE COMPLEX SYSTEM C A * X = B OR CTRANS(A) * X = B C USING THE FACTORS COMPUTED BY CGECO OR CGEFA. C C ON ENTRY C C A COMPLEX(LDA, N) C THE OUTPUT FROM CGECO OR CGEFA. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A C C N INTEGER C THE ORDER OF THE MATRIX A C C IPVT INTEGER(N) C THE PIVOT VECTOR FROM CGECO OR CGEFA. C C B COMPLEX(N) C THE RIGHT HAND SIDE VECTOR. C C JOB INTEGER C = 0 TO SOLVE A*X = B, C = NONZERO TO SOLVE CTRANS(A)*X = B WHERE C CTRANS(A) IS THE CONJUGATE TRANSPOSE. C C ON RETURN C C B THE SOLUTION VECTOR X C C ERROR CONDITION C C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGUI.ARITY C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER C SETTING OF LDA. IT WILL NOT OCCUR IF THE SUBROUTINES ARE C CALLED CORRECTLY AND IF CGECO HAS SET RCOND.GT. 0.0 C OR CGEFA HAS SET INFO.EQ. 0. C C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX C WITH P COLUMNS C CALL CGECO(A,LDA,N,IPVT,RCOND,Z) C IF (RCOND IS TOO SMALL) GO TO... C DO 10 J = 1, P C CALL CGESL(A,LDA,N,IPVT,C(1,J),O) C 10 CONTINUE C C LINPACK. THIS VERSION DATED 07/14/77 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. C C SUBROUTINES AND FUNCTIONS C C BLAS CAXPY,CDOTC 100

C FORTRAN CONJG C C INTERNAL VARIABLES C COMPLEX CDOTC,T INTEGER K,KB,L,NM1 C NM1 = N - 1 IF (JOB.NE. O) GO TO 50 C C JOB 0, SOLVE A * X = B C FIRST SOLVE L*Y = B C IF (NM1.LT. 1) GO TO 30 DO 20 K = 1, NM1 L = IPVT(K) T = B(L) IF (L.EQ. K) GO TO 10 B(L) = B(K) B(K) = T 10 CONTINUE CALL CAXPY(N-K,T,A(K+1,K),1,B(K+1),1) 20 CONTINUE 30 CONTINUE C C NOW SOLVE U*X = Y C DO 40 KB = i, N K = N + 1 - KB B(K) = B(K)/A(K,K) T = -B(K) CALL CAXPY(K-1,T,A(1,K),1,B(1),1) 40 CONTINUE GO TO 100 50 CONTINUE C C JOB = NONZERO, SOLVE CTRANS(A) * X = B C FIRST SOLVE CTRANS(U)*Y = B C DO 60 K = 1, N T = CDOTC(K-1,A(1,K),i,B(1),i) B(K) = (B(K) - T)/CONJG(A(K,K)) 60 CONTINUE C C NOW SOLVE CTRANS(L)*X = Y C IF (NM1.LT. 1) GO TO 90 DO 80 KB = 1, NM1 K = N - KB B(K) = B(K) + CDOTC(N-K,A(K+1,K),1,B(K+1),1) L = IPVT(K) IF (L.EQ. K) GO TO 70 T = B(L) B(L) = B(K) B(K) = T 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN 101

END C C C****************** *****t*** ********** **************** **************C C C SUBROUTINE CAXPY(N,CA,CX,INCX,CY,INCY) C C C** *********** ********** ******************************** ******* *c C C C NAASA 1.1.014 CAXPY FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CX(1),CY(1), CA INTEGER I,INCX,INCY,IX,IY,N C IF(N.LE.O)RETURN IF (ABS(REAL(CA)) + ABS(AIMAG(CA)).EQ. 0.0 ) RETURN IF(INCX.EQ.1.AND.IINCY.EQ.1)GOTO 20 C CCC Code for unequal increments or equal increments CCC Not equal to 1 C IX = 1 IY = 1 IF(INCX.LT.O)IX = (-N+1)*INCX + 1 IF(INCY.LT.O)IY = (-N+1)*INCY + 1 DO 10 I = 1,N CY(IY) = CY(IY) + CA*CX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C CCC Code for both increments equal to 1 C 20 DO 30 I = 1,N CY(I) = CY(I) + CA*CX(I) 30 CONTINUE RETURN END C C COMPLEX FUNCTION CDOTC(N,CX,INCX,CY,INCY) C C C NAASA 1.1.012 CDOTC FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C C FORMS THE DOT PRODUCT OF TWO VECTORS, CONJUGATING THE FIRST C VECTOR. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CX(1),CY(1),CTEMP INTEGER I,INCX,INCY,IX,IY,N C CTEMP = (0.0,0.0) CDOTC = (0.0,0.0) IF(N.LE.O)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GOTO 20 C 102

CCC Code for unequal increments or equal increments CCC Not equal to 1 C IX = 1 IY = 1 IF(INCX.LT.O)IX = (-N+1)*INCX + 1 IF(INCY.LT.O)IY = (-N+1)*INCY + 1 DO 10 I = 1,N CTEMP = CTEMP + CONJG(CX(IX))*CY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE CDOTC = CTEMP RETURN C CCC Code for both increments equal to 1 C 20 DO 30 I = 1,N CTEMP = CTEMP + CONJG(CX(I))*CY(I) 30 CONTINUE CDOTC = CTEMP RETURN END C C C** * ** **4**** ** ** ******** **** * ***** *********************************C C C SUBROUTINE CSCAL(J[,CA,CX, INCX) C C C*e****************************************************************C C C C NAASA 1.1.019 CSCAL FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C C SCALES A VECTOR BY A CONSTANT. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CA,CX(1) INTEGER I,INCX,N,NINCX C IF(N.LE.O)RETURN IF(INCX.EQ.1)GOTO 20 C CCC Code for increment not equal to 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX CX(I) = CA*CX(I) 10 CONTINUE RETURN C CCC Code for increment equal to 1 C 20 DO 30 I = 1,N CX(I) = CA*CX(I) 30 CONTINUE RETURN END C C C************************ ******************************************C C C SUBROUTINE CSSCAL (:, SA,CX, INCX) 103

c c C***********************N4*************************************** ***c C C C NAASA 1.1.018 CSSCAL FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C C SCALES A COMPLEX VECTOR BY A REAL CONSTANT. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CX(1) REAL SA INTEGER I,INCX,N,NINCX C IF(N.LE.O)RETURN IF(INCX.EQ.1)GOTO 20 C CCC Code for increment not equal to 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX CX(I) = CMPLX(SA*REAL(CX(I)),SA*AIMAG(CX(I))) 10 CONTINUE RETURN C CCC Code for increment equal to 1 C 20 DO 30 I = 1,N CX(I) = CMPLX(SA*REAL(CX(I)),SA*AIMAG(CX(I))) 30 CONTINUE RETURN END C C C**$**^*****************^***********************+**$$$***$$***$$***c INTEGER FUNCTION ICAMAX(N,CX,INCX) C0**i4(4E*********^~***4^**4i**************$$*$$****$$**$*$$**$$*$** C C C NAASA 1.1.021 ICAMAX FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CX(1) REAL SMAX INTEGER I,INCX,IX,N COMPLEX ZDUM REAL CABS1 CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM)) C ICAMAX = 1 IF(N.LE.1)RETURN IF(INCX.EQ.1)GOTO 20 C CCC Code for increment not equal to 1 C IX = 1 SMAX = CABS1(CX(1)) IX = IX + INCX DO 10 I = 2,N IF(CABS1(CX(IX)).LE.SMAX) GO TO 5 ICAMAX = I SMAX = CABS1(CX(IX)) 104

5 IX = IX + INCX 10 CONTINUE RETURN C CCC Code for increment equal to 1 C 20 SMAX = CABS (CX(1)) DO 30 I = 2,N IF(CABS1(CX(I)).LE.SMAX) GO TO 30 ICAMAX = I SMAX = CABS1(CX(I)) 30 CONTINUE RETURN END C C C** *** *********** ******** ** *************************************** *C REAL FUNCTION SCASUM(N,CX,INCX) C*************************** *********************** ** *** * ****C C C C NAASA 1.1.010 SCASUM FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C C TAKES THE SUM OF THE ABSOLUTE VALUES OF A COMPLEX VECTOR AND C RETURNS A SINGLE PRECISION RESULT. C JACK DONGARRA, LINPACK, 6/17/77. C COMPLEX CX(1) REAL STEMP INTEGER I,INCX,N,NINCX C SCASUM = O.OEO STEMP = O.OEO IF(N.LE. O)RETURN IF(INCX.EQ.1)GOTO 20 C CCC Code for increment not equal to 1 C NINCX = N*INCX DO 10 I = 1,NINCX,INCX STEMP = STEMP + ABS(REAL(CX(I))) + ABS(AIMAG(CX(I))) 10 CONTINUE SCASUM = STEMP RETURN C CCC Code for increment equal to 1 C 20 DO 30 I = 1,N STEMP = STEMP + ABS(REAL(CX(I))) + ABS(AIMAG(CX(I))) 30 CONTINUE SCASUM = STEMP RETURN END C C C The following subroutines are standard LINPACK routines C C to perform L-U decomposition and back substitution on a C C single precision complex matrix. See CC-Memo 407 sec 2.1 C C for documentation on these routines. C C C C C************ * * **** ************** C C 105

SUBROUTINE CGECO(A,LDA,N,IPVT,RCOND,Z) C C C C C NAASA 2.1.042 CGECO FTN-A 05-02-78 THE UNIV OF MICH COMP CTR C C INTEGER LDA,N,IPVT(1) COMPLEX A(LDA,1),Z(1) REAL RCOND C C CGECO FACTORS A COMPLEX MATRIX BY GAUSSIAN ELIMINATION C AND ESTIMATES THE CONDITION OF THE MATRIX. C C IF RCOND IS NOT NEEDED, CGEFA IS SLIGHTLY FASTER. C TO SOLVE A*X = B, FOLLOW CGECO BY CGESL. C TO COMPUTE INVERSE(A)*C, FOLLOW CGECO BY CGESL. C TO COMPUTE DETERMINANT(A), FOLLOW CGECO BY CGEDI. C TO COMPUTE INVERSE(A), FOLLOW CGECO BY CGEDI. C C ON ENTRY C C A COMPLEX(LDA, N) C THE MATRIX TO BE FACTORED. C C LDA INTEGER C THE LEADING DIMENSION OF THE ARRAY A. C C N INTEGER C THE ORDER OF THE MATRIX A C C ON RETURN C C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS C WHICH WERE USED TO OBTAIN IT. C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. C C IPVT INTEGER(N) C AN INTEGER VECTOR OF PIVOT INDICES. C C RCOND REAL C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A C FOR THE SYSTEM A*X = B, RELATIVE PERTURBATIONS C IN A AND B OF SIZE EPSILON MAY CAUSE C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION C 1.0 + RCOND.EQ. 1.0 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING C PRECISION. IN PARTICULAR, RCOND IS ZERO IF C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE C UNDERFLOWS. C C Z COMPLEX(N) C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) C C LINPACK. THIS VERSION DATED 07/14/77 106

C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. C C SUBROUTINES AND FUNCTIONS C C LINPACK CGEFA C BLAS CAXPY,CDOTC,CSSCAL,SCASUM C FORTRAN ABS,AIMAG,AMAX1,CMPLX,CONJG,REAL C C INTERNAL VARIABLES C COMPLEX CDOTC,EK,T,WK,WKM REAL ANORM,S,SCASUM,SM,YNORM INTEGER INFO,J,K,KB,KP1,L C COMPLEX ZDUM,ZDUM1,ZDUM2,CSIGN1 REAL CABS1 CABS1(ZDUM) = ABS(REAL(ZDUM)) + ABS(AIMAG(ZDUM)) CSIGN1(ZDUM1,ZDUM2) = CABS1(ZDUM1)*(ZDUM2/CABS1(ZDUM2)) C CCC Compute 1-NORM of A C ANORM = O.OEO DO 10 J = 1, N ANORM = AMAX1(ANOR[,SCASUM(N,A(1,J),)) 10 CONTINUE C CCC Factor C CALL CGEFA(A,LDA,N,IPVT,INFO) C C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))). C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND CTRANS(A)*Y = E C CTRANS(A) IS THE CONJUGATE TRANSPOSE OF A. C THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL C GROWTH IN THE ELEMENTS OF W WHERE CTRANS(U)*W = E. C THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW. C C SOLVE CTRANS(U)*W = E C EK = CMPLX(1.OEO,O.OEO) DO 20 J = 1, N Z(J) = CMPLX(O.OEO,O.OEO) 20 CONTINUE DO 100 K = 1, N IF (CABS1(Z(K)).NE. O.OEO) EK = CSIGN1(EK,-Z(K)) IF (CABS1(EK-Z(K)).LE. CABS1(A(K,K))) GO TO 30 S = CABS1(A(K,K))/CABS1(EK-Z(K)) CALL CSSCAL(N,S,Z,1) EK = CMPLX(S,O.OEO)*EK 30 CONTINUE WK = EK - Z(K) WKM = -EK - Z(K) S = CABS1(WK) SM = CABS1(WKM) IF (CABS1(A(K,K)).EQ. O.OEO) GO TO 40 WK WK/CONJG(A(K,K)) WKM = WKM/CONJG(A(K,K)) GO TO 50 40 CONTINUE WK = CMPLX(1.OEO,O.OEO) 107

WIN = CMPLX(:L.0E0,0.0E0) 50 CONTINUE KP1 = K + 1 IF (KP1.GT. N) GO TO 90 DO 60 3 = KPI, N SM = SN + CABS1(Z(J)+WKN*CONJG(A(K,J))) Z(3 = Z(J) + WK*CONJG(A(K,J)) S = S + CAkBSI(Z(J)) 60 CONTINUE IF (S.GE. SN) GO TO 80 T = UKM! - K UK = UKM! DO TO J KP1, N Z(JM Z(3 + T*CONJGCA(K,J)) 70 CONTINUE 80 CONTINUE 90 CONTINUE Z(K) = WK 100 CONTINUE S = 1.OEO/SCASUN(N,Z,1) CALL CSSCAL(N,S,Z,I.) C CCC Solve CTRANS(L)*Y = V C DO 120 KB = 1, N K = N + 1 - KB IF (K.LT. N) Z(K) = Z(K) + CDOTC(N-K,A(K+1,K),l,Z(K+1), 1) IF (CABS1(z(K)).LE. 1.OEO) GO TO 110 S = 1.OEO/CABSl(Z(K)) CALL CSSCAL(I',S,Z,l) 110 CONTINUE L = IPVT(K) T = ZCL Z(L) = Z(K) Z(K) = T 120 CONTINUE S = 1.OEO/SCASUMN(N,Z,1) CALL CSSCAL(N,S,Z,l) C YNORM = 1.OEO C CCC Solve L*V =Y C DO 140 K= 1, N L = IPVT(K) T = Z(L) Z(L = Z(K) Z(K) = T IF (K.LT. N) CALL CAXPY(N-K,T,A(K+1,K),l,Z(K+1),l) IF (CABS1(Z(K)).LE. 1.OEO) GO TO 130 S =1.OEO/CABS1(ZCK)) CALL CSSCAL(N,S,Z,l) YNORM = S*YNORIM 130 CONTINUE 140 CONTINUE S = 1.OEO/SCASUN(N,Z,1) CALL CSSCALCN,S,Z,l) YNORN = S*YNORM C CCC Solve U*Z = V 108

C DO 160 KB = 1, N K = N + 1 - KB IF (CABS1CZ(K)).LE. CABS1CA(K,K))) GO TO 150 S =CABS1CACK,K))/CABS1(Z(K)) CALL CSSCALCN,S,Z,l) YNORMI = S*YNORM 150 CONTINUE IF CcABS1(A(K,K:)) -NE. 0.OEO) Z(K = Z(K)/A(K,K) IF CCABS1(A(K,K:)).EQ. 0.0E0) Z(K)= CMPLX(1.OEO,0.OEO) T = -Z(K) CALL CAXPY(K-1,Tr,A(1,K),l,ZC1),l) 160 CONTINUE C MAKE ZNORM = 1.0 S = 1.0E0/SCASUM(N,.Z,l) CALL CSSCALCN,S,Z,IL) YNORM = S*YNORMN IF CANORM NE. 0.OEO) MCOND = YNORM/ANORM IF (ANORM EQ. 0.0E0) RCOND = 0.0E0 RETURN END 109

DESCRIPTION AND LIST OF CODE VSIE-ISO Jian-Ming Jin Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, Michigan 48109 ABSTRACT A FORTRAN program named as VSIE-ISO and developed by the author is described briefly and listed thereafter. 1. Objective of VSIE-ISO VSIE-ISO is a computer code written in FORTRAN language and developed for computing electromagnetic scattering from an inhomogeneous penetrable cylinder for both TE and TMN incidence. VSIE-ISO uses quadrilateral isoparametric elements and point-matching technique to solve the volume-surface integral equation (VSIE). For the formulation of VSIE, please read the first paper in this report. For the description of the numerical implementation, please read the third paper. VSI'E-ISO was used to generate the data presented in the third paper. 110

Due to the success of VSIE-ISO, we will extend it to include the treatment of conducting surfaces and discontinous interfaces. The resulting program should be able to treat complex problems as VSIEM does, but with a much higher accuracy. 2. Input Data All input data are provided by file INPUT except for relative permittivity and permeability and their derivatives which are specified by subroutines EMFUN and DRFUN. A detailed description of the input data is given in the main program. A sample file of INPUT for the geometry shown in Figure 1 is listed below. 3 5 1 7 Figure 1: A five quadrilateral element model with 20 nodes for a circular region. DIELECTRIC CIRCULAR CYLINDER 0 20.0 1 1 0. 180. 30. 180. 111

20 8 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 1 2 3 4 5 1 2 3 4 5 4 1. O. 0.70711 0.70711 0. 1. -0.70711 0.70711 -1. 0. -0.70711 -0.70711 0. -1. 0.70711 -0.70711 0.65 0. 0. 0.65 -0.65 0. 0. -0.65 0.3 0. 0.21213 0.21213 0. 0.3 -0.21213 0.21213 -0.3 0. -0.21213 -0.21213 0. -0.3 0.21213 -0.21213 15 13 1 3 14 9 17 15 3 5 16 10 19 17 5 7 18 11 13 19 7 1 20 12 13 15 17 19 14 16 1 3 2 3 5 4 5 7 6 7 1 8 2 10 4 11 6 12 8 9 18 20 3. Output Data The output data are written in file OUTPUT. The data contain the geometry and discretization information and the computed result. Following is a sample output. DIELECTRIC CIRCULAR CYLINDER TOTAL NUMBER OF' DISCRETE POINTS I X(I) Y(I) 20 1 2 3 4 5 6 7 8 9 1.00000 0. 70711 0.00000 -0.70711 -1.00000 -0.70711 0.00000 0.70711 0. 65000 0.00000 0.70711 1.00000 0.70711 0.00000 -0.70711 -1.00000 -0.70711 0.00000 112

10 11 12 13 14 15 16 17 18 19 20 0.00000 -0.65000 0.00000 0.30000 0.21213 0.00000 -0.21213 -0.30000 -0.21213 0.00000 0.221213 0.65000 0.00000 -0.65000 0.00000 0.21213 0.30000 0.21213 0.00000 -0.21213 -0.30000 -0.21213 TOTAL NUMBER OF AREA ELEMENTS 5 L N(1,L) N(2,L) N(3,L) N(4,L) N(5,L) N(6,L) N(7,L) N(8,L) 1 2 3 4 5 15 17 19 13 13 13 15 17 19 15 1 3 5 7 17 3 5 7 1 19 14 16 18 20 14 9 10 11 12 16 2 4 6 8 18 10 11 12 9 20 TOTAL NUMBER OF LINE ELEMENTS L N(1,L) N(2,L) N(3,L) 4 1 1 3 2 3 5 3 5 7 4 7 1 2 4 6 8 KEY PARAMETERS INCIDENT POLARIZATION NUMBER OF INCIDENT FIELD DIRECTIONS NUMBER OF BISTATIC DIRECTIONS WAVELENGTH FIELD VALUES AT NODES REAL IMAGINARY NODE 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 0.1071E-01 0.3847E-01 0.1417E+00 0.1945E+00 0.1983E+00 0.1945E+00 0.1417E+00 0.3847E-01 -0.5831E-01 0.1395E-01 0.7612E-01 0.1395E-01 -0.6544E-01 -0.6291E-01 -0.5218E-01 -0.4194E-01 -0.3658E-01 -0.8549E-01 -0.7804E-01 -0.1951E-01 0.2023E-01 0.9312E-01 0.2023E-01 -0.1951E-01 -0.7804E-01 -0.9881E-02 -0.6613E-01 -0.8084E-01 -0.6613E-01 0.2954E-01 0.1555E-01 -0.1739E-01 -0.4392E-01 -0.5415E-01 E 1 7 20.00000 MAGNITUDE 0.0862 0.0870 0.1431 0.1956 0.2191 0.1956 0.1431 0.0870 0.0591 0.0676 0.1110 0.0676 0.0718 0.0648 0.0550 0.0607 0.0653 113

18 19 20 -0.4194E-01 -0.4392E-01 -0.5218E-01 -0.1739E-01 -0.6291E-01 0.1555E-01 0.0607 0.0550 0.0648 BISTATIC SCATTERING CROSS SECTION FOR INCIDENT FIELD DIRECTION= 180.00 THETA 10*LOG(SIGMA/LAMBDA) PHASE 0.00 -4.03 -41.2 30.00 -4.19 -40.8 60.00 -4.63 -39.7 90.00 -5.25 -37.9 120.00 -5.89 -35.8 150.00 -6.38 -34.0 180.00 -6.56 -33.3 4. Program Description Following is a brief description of the main program and subroutines developed for VSIE-ISO. All others not described here are standard subroutines and since they are already listed in VSIEM they will not be listed here. Note that the equations referenced below are referred to those in the third paper. VSIE-ISO - EMFUM - DRFUN - MATRIXi MATRIX2 FRAF1 - FRAF2 - CALAE - CALAE1 - Main program acting as a driver for subroutines. Specifies Or and [r at a given point (x, y). Specifies Ver and V/yr at a given point (xy). Computes the contribution of the area integral to the system matrix. Computes the contribution of the line integral to the system matrix. Computes the contribution of the area integral to the far field. Computes the contribution of the line integral to the far field. Computes A(rlr') in (5). Computes A(rlr') in (5) for r at infinity. 114

CALBE - CALBE1 - CALNE - CALNEXI - CALNETA - CALNE1 - CALNEXI1 - Computes B(rlr') in (5). Computes B(rlr') in (5) for r at infinity. Calculates Ne(,r/7) in (14). Calculates ONf/O9. Calculates ONfe/lT. Calculates LP() in (9). Calculates OLP/OD. 115

C******************** *** ******** **************************************c C PROGRAM VSIE-ISO C C*********************** ***** *****************************************c C C C THIS PROGRAM USES QUADRILATERAL ISOPARAMETRIC ELEMENTS AND POINT- C C MATCHING TECHNIQUE TO SOLVE VOLUME-SURFACE INTEGRAL EQUATION FOR C C ELECTROMAGNETIC SCATTERING BY A DIELECTRIC CYLINDER. C C C C*** ** ** * *** ************ ******** ***** *********** **************** ****** **C C***** *************** **** ******* ************ ************************ *******C C INPUT FORMAT FOR PROGRAM VSIE-ISO ----VERSION OF AUG. 20, 1988 C C***** ************************ ********************** ******* ****** ** *****C C CARD 1 FORMAT (18A4) TITLE CARD; USE UP TO 72 COLUMNS C C***** *** ***************************************************************C C CARD 2 FORMAT (I3,F10.5) KODE,WAVE C C KODE=O COMPUTES BISTATIC SCATTERING PATTERN C C KODE=1 COMPUTES BACKSCATTERING PATTERN C C WAVE WAVELENGTH C C*********************************************************************C C CARD 3 FORMAT (I3,I3:4F10.5) IPP,IOPT,FIRST,LAST,INK,CANG C C IPP=1 E-POLARIZATION C C IPP=2 1I-POLARIZATION C C IOPT=O INTERIOR FIELDS NOT PRINTED C C IOPT=i INTERIOR FIELDS PRINTED C C FIRST INITIAL SCATTERING AND INCIDENCE ANGLE C C LAST FINAL ANGLE C C INK ANGULAR INCREMENT C C CANG ANGLE FOR BISTATIC COMPUTATION C C** * ** ******** ******* *** * **** ********* ***** ******** **** * ***** ***** * ****C C FORMAT FOR INPUT NUMBER OF NODES AND ELEMENTS C C*********************** ********************************** *************C C CARD 4 FORMAT (213) N,M C C N NUMBER OF NODES C C M NUMBER OF QUADRILATERAL ELEMENTS C C*****************************************************************+*****C C FORMAT FOR INPUT NUMBER OF SURFACE NODES AND SEGMENTS C C********************************************************* ************C C CARD 5 FORMAT (213) Nl,Ml C C N1 NUMBER OF NODES ON THE SURFACE C C M1 NUMBER OF QUADRATIC SEGMENTS ON THE SURFACE C C************************:**************************************** * ** **C C DO I=1,N TO READ THE NODE POSITION C C**********************************************************************C C CARD 6 FORMAT (I3,2F10.5) I,X(I),Y(I) C C I NODE NUMBER C C X(I) NODE X-COORDINATE C C Y(I) NODE Y-COORDINATE C C REPEAT CARD 6 FOR EACH NODE C C DO I=1,M TO READ THE NODE-ELEMENT RELATION C C CARD 7 FORMAT (I3,2X,8I3) I,(NUM(J,I),J=1,8) C C I ELEMENT NUMBER C C NUM(J,I) GLOBAL NUMBER OF NODE J IN ELEMENT I C C REPEAT CARD 7 FOR EACH ELEMENT C C********************** DO I=,M TO READ THE NODE-SEGMENT RE**********************************************TION C C DO 1=1,Mi TO READ THE NODE-SEGMENT RELATION C 116

C CARD 8 FORMAT (I3,21,3I3) I,(NUM1(J,I),J=1,3) C C I SEGMENT NUMBER C C NUM1(J,I) GLOBAL NUMBER OF NODE J AT SEGMENT I C C REPEAT CARD 8 FOR EACH SEGMENT C COMPLEX A(100,100),FIIC(100),FAR(91) REAL X(l00),Y(l00),LAST,INK INTEGER NUMC8,50),NUM1((3,50),KPVT(100) DIMENSION ID(18),IPOL(2) DATA RED,IPOL/0.OIL745329,4HEEEE,4HHHHH/ OPEN(UNIT=5,FILE=)INPUT') C OPEN(UNIT=6,FILE=:'OUTPUT') READ(5,1000) ID WRITE(6,1005) ID READ(5,1010) KODEWAVE READ(5,1020) IPP,1OPT,FIRST,LAST,INKCANG READ(5,1030) N,M READ(5,1030) N1,M1 DO I=1,1 READC5,1040) X(I),Y(I) END DO DO I=1,M READ(5,1050) (NUM(J,I),J=l,8) END DO DO I=1,Ml READ(5,1050) (NUMi(J,I),J1l,3) END DO C......PRINT INPUT DATA WRITE(6,1060) N WRITE(6,10) CI,X(I),Y(I),I=1,N) 10 FORMAT(9X,I3,5X,FIO.5,5XF10.5) WRITE(6,1070) M DO 20 L=1,M 20 WRITE(6,30) L,(NUMII,L),11,8) 30 FORMAT(1OX,I2,3X,8C3X,I3,2X)) WRITE(6,1080) Ml DO 40 L=1,Ml 40 WRITE(6,50) L,(NU]M1(I,L),I=1,3) 50 FORMAT(10X,I2,3X,3(3X,13,2X)) C.....CALL SUBROUTINES TO FORM SYSTEM MATRIX CALL MATRIX1(A,N,,X,Y, NUM,WAVE,IPP) CALL MATRIX2(A,N,11,M1,X,Y,NUM1,WAVE,IPP) C......CALL SUBROUTINE TO SOLVE THE MATRIX CALL CGECO(A,100,N,KPVT,RC,FINC) XK=2.*3.141592654/WAVE IF CKODE.NE.0) GO TO 70 NINC=1 NBIT=1+IFIX((LAST-FIRST)/INK) GO TO 80 70 NBIT=l NINC=l+IFIX((LAST-FIRST)/INK) 80 WRITE (6,400) IPOL(IPP),NINC,NBIT,WAVE IF(KODE.EQ.1) WRITE(6,600) IF(KODE.EQ.1) CANG=FIRST 90 TETA=RED*CANG CT=COS(TETA) ST=SIN(TETA) DO 100 I=1,N 117

HOLD=XK*(CT*X(I)+ST*Y(I)) 100 FINC(I)=CMPLX(COS(HOLD),SIN(HOLD)) C..... COMPUTE NEAR FIELD (NODAL FIELD VALUES) CALL CGESL(A,l00,'ff,KPVT,FINC,0) IF(IOPT.EQ.0) GOTO 101 WRITE(6,700) DO 98 I=1,1 98 WRITE(6,99) I,FINC(I),CABS(FINC(I)) 99 FORMAT(5X,13,5X,2E1E4.4,5X,F12.4) 101 CONTINUE CALL FARF1(FAR,FINC,N,M,X,YN,WAVE,FIRST,INKNBITIPP) CALL FARF2(FAR,FINIC,N,N1,Ml,X,Y,NUM1,WAVE,FIRST,INKNBIT,IPP) PI=3.141592654 IF(KODE.EQ.0) WRITE(6,610) CANG DO 102 I=1,NBIT ANGLE=FIRST+INK*FL.OAT (I-1) PMAG=CABS(FAR(I)) SIGMA=1O.*LOG10(2.*PMAG*PMAG/PI) PPHS=ATAN2(AIMAG(FAR(I)),REAL(FAR(I)))*180./PI WRITE(6,650) ANGLE,SIGMA,PPHS 102 CONTINUE IF(KODE.EQ.1) THEN CANG=CANG+INK FIRST=CANG IF(CANG.LE.LAST) GOTO 90 ELSE END IF 1000 FORMAT (18A4) 1005 FORMAT (IX,18A4) 1010 FORMAT (13,F10.5) 1020 FORMAT (13,13,4F10.5) 1030 FORMAT (213) 1040 FORMAT (3X,2F10.5) 1050 FORMAT (5X,8I3) 1060 FORMAT(10X,'TOTAL NUMBER OF DISCRETE POINTS',SX,I3/ & lOX,' I',10X,'X(I)',1OX,'YCI)'/) 1070 FORMAT(//1OX,'TOTAL NUMBER OF AREA ELEMENTS',5X,I3/ & lOX,' L',5X,'N(1,L)',2X,'N(2,L) ',2X,'N(3,L)',2X,'N(4,L), & 2X,'N(5,L)',2X,'N(6,L)',2X,'N(7,L)',2X,'N(8,L)'/) 1080 FORMAT(//10X,'TOTAL NUMBER OF LINE ELEMENTS',5X,I3/ & lox,' L',5X,'N(1,L)',2X,'N(2,L)',2X,'N(3,L)'/) 400 FORMAT (//25X,14HKEY PARAMETERS/I &1OX,21HINCIDENT POLARIZATION,22X,lAl/ &10X,35HNUMBER OF INCIDENT FIELD DIRECTIONS,I9/ &10X,29HNUMBER OF BISTATIC DIRECTIONSI15/ &1OX,1OHWAVELENGTH,F34.5) 600 FORMAT (///,20X,28HBACKSCATTERING CROSS SECTION/17X, &36HTHETA 10*LOG(SIGMA/LAMBDA) PHASE/) 610 FORMAT (///,19X,33HBISTATIC SCATTERING CROSS SECTION/18X, &29HFOR INCIDENT FIELD DIRECTION=,F7.2/17X, &36HTHETA 10*LOG(SIGMA/LAMBDA) PHASE/) 650 FORMAT (9X,F13.2,F15.2,F16.1) 700 FORMAT(//20X,'FIELD VALUES AT NODES',! &5X,'NODE REAL IMAGINARY MAGNITUDE'!) STOP END 118

SUBROUTINE EMFUN(:X,Y,EPS,MU,IPP) COMPLEX EPS,MU,EXCH C..... SPECIFY PERMITTIVITY AND PERMEABILITY EPS=(72.,-162.) MU=(1.0,0.O) IF(IPP.EQ.2) RETURN EXCH=EPS EPS=MU MU=EXCH RETURN END SUBROUTINE DRFUN(X,Y,DXEPS,DYEPS,IPP) COMPLEX DXEPS,DYEPS,DXMU,DYMU C..... SPECIFY THE DERIVATIVES OF PERMITTIVITY AND PERMEABILITY DXEPS=(O.,O.) DYEPS=(O.,O.) DXMU=(O.,O.) DYMU=(0.,O.) IF(IPP.EQ.2) RETURN DXEPS=DXMU DYEPS=DYMU RETURN END 119

SUBROUTINE MATRIXI (A,N,M,X,Y,NUI,WAVE, IPP) COMPLEX A(100,100),AE,AIJ,EPS,MU,B(8) REAL X(N),Y(N),W(,3),T(3),NE(8),NEXI(8),NETA(8) REAL XlC8),YlC8) INTEGER NUJN(8,N) DATA W/0.55555555,56,0.8888888889,0.5555555556/ DATA T/-0.77459666592,*0.,0.7745966692/ DO 100 I=1,N DO 100 3=1,N 100 A(I,J)=CMPLX(O. 10.) DO 120 I=1,N CALL EMFUN(X(I),Y(I),EPS,MU,IPP) 120 A(I,I)=l./EPS DO 900 L=1,M DO 150 Il=1,8 Xl1(Il) =X(NUM (I1 L;)) 150 Yl(Il)=Y(NUJM(Il,L;)) DO 700 3=1,3 WJ=W(J) XI=T(J) DO 600 K=1,3 WK=W(K) ETA=T (K) CALL CALNE(NE,XI,E —TA) Xx=0. DO 200 I1=1,8 200 YY=YY+NE(Il)*Yl(I1.) CALL CALNEXI (NEXI, XI, ETA) CALL CALNETA (NETA, XI, ETA) xxI=0. YXI=0. XETA=O. YETA=O. DO 300 I1=1,8 YXI=YXI+Y1l(Il)*NEXI(Il) XETA=XETA+X1l(Il)*NETA(Il) 300 YETA=YETA+Y1l(Il)*NETA(Il) DETJ=ABS (XXI*YETA-XETA*YXI) DO 1000 II=1,N CALL CALAE(AE,X(II),Y(II),XX,YY,WAVE,IPP) DO 400 1=1,8 400 B(I)=WJ*WK*NE(I)*AE*DETJ DO 1000 1=1,8 JJ=NUM(I,L) A(II,JJ)=A(II,JJ)+B(I) 1000 CONTINUE 600 CONTINUE 700 CONTINUE 800 CONTINUE 900 CONTINUE RETURN END 120

SUBROUTINE MATRIX:2(A,N,N1,N1,X,Y,NUM1,WAVE,IPP) COMPLEX A(100,100.),BE,AIJ,EPS,MU,B(3) REAL X(N),Y(N),W(4I),T(4),NE(3),NEXI(3),X1(3),Y1(3) INTEGER NUMi1(3,H1.) DATA W/0.34785484!51,O.6521451549,0.6521451549,0.3478548451/ DATA T/-0.8611363:116,-0.3399810436,O.3399810436,0.8611363116/ DO 100 I=1,Nl CALL ENFUNCICI),Y(I),EPS,MU,IPP) 100 A(I,I)=A(I,I)+O.5*(l.-1./EPS) DO 900 L=1,Ml DO 150 Il=1,3 Xl (Il)=X(NUM1 (Ii,L')) 150 Y1(I1)=YCNTJH1(I1L')) DO 700 3=1,4 WJ=W (3) XI=T(J) CALL CALNE1(NE,XI:) xx=0. yy=0. DO 200 I1=1,3 XX=XX+NE(Il) *Xl (Il) 200 YY=YY+NE(I1)*Yl(IIL) CALL CALNEXIl1(NEXT,XI) XxI=0. YXI=0. DO 300 I1=1,3 XXI=XXI-X1l(Il)*NEXkI(I1) 300 YXI=YXI+Y1l(Il) *NEI:ICIl) DETJ=SQRT (XXI*XXI+-YXI*YXI) XXI=XXI/DETJ YXI=YXI/DETJ DO 1000 II=1,N CALL CAL`BE(`BE,l(I-),Y (Ul),II,NYY,IU 1,NI,lJkJE,TPP) DO 400 1=1,3 400 B (I) =WJ *HE (I) * BE*LDETJ DO 1000 1=1,3 JJ=NUM1 (I,L) ACII,JJ)=A(II,JJ)+B(I) 1000 CONTINUE 700 CONTINUE 800 CONTINUE 900 CONTINUE RETURN END 121

SUBROUTINE FARF1(FAR,F,N,M,X,Y,NUMN,WAVE,FIRST,DANG,NBIT,IPP) COMPLEX F(100),AE,AIJ,EPS,IJ,B(8),FAR(l) REAL X(N),Y(N), WC3),T(3),NE(8),NEXI(8),NETA(8) REAL XlC8),Y1(8) INTEGER NUM(8,M) DATA W/0.5555555556,0.8888888889,0.5555555556/ DATA T/-O.77459666392,0.,0.7745966692/ DO 170 I=1,NBIT 170 FAR(I)=CMPLXCO.,0.) DO 900 L=1,N DO 150 II=1,8 Xl (Il)=X(NUN(I1,L:)) 150 Yl(Il)=Y(NUM(Il,L:)) DO 700 3=1,3 WJ=W(J) XI=T(J) DO 600 K=1,3 WK=W(K) ETA=T(K) CALL CALNE(NE,XI,ETA) xX=O. YY=O. DO 200 11=1,8 XX=XX+NE(I1)*X1(II) 200 YY=YY+NE(Il) *Y(II) CALL CALNEXI(NEXI, XI,ETA) CALL CALNETA(NETA,XI,ETA) XXI=O. YXI=0. XETA=O. YETA=O. DO 300 11=1,8 XXI=XXI+X1(I11)*NEXICI1) YXI=YxI+Yl(Il)*NEX.ICIl) XETA=XETA+Xl CIl)*BETA (Il) 300 YETA=YETA+Yl(Il)*N'ETA(Il) DETJ=ABS(XXI*YETA-xETA*YXI) DO 1000 II=1,NBIT ANG=(FIRST+FLOAT(II-1)*DANG)*3.141592654/180. CALL CALAE1(AE,ANG,XX,YY,WAVE,IPP) DO 400 1=1,8 JJ=NUM(I,L) 400 FAR(II)=FAR(II)+WJ*WK*NE(I)*AE*DETJ*F(JJ) 1000 CONTINUE 600 CONTINUE 700 CONTINUE 800 CONTINUE 900 CONTINUE RETURN END 122

SUBROUTINE FARF2(FAR,F,N,N1,Ml,X,Y,NUM1,WAVE,FIRST,DANG,NBIT,IPP) COMPLEX F(100),BE,AIJ,EPS,MU,FAR(1) REAL X(N),Y(N),W(4),T(4),NE(3),NEXI(3),X1(3),Y1(3) INTEGER NUN1(3,M1) DATA W/0.3478548451,0.6521451549,0.6521451549,0.3478548451/ DATA T/-0.8611363116,-0.3399810436,0.3399810436,0.8611363116/ DO 900 L=1,M1 DO 150 I1=1,3 X (I1)=X(NUM1 (Il,L)) 150 Yl(I1)=Y(NUM1(I1,L)) DO 700 J=1,4 WJ=W(J) XI=T(J) CALL CALNE (NE,XI) XX=O. YY=O. DO 200 I1=1,3 XX=XX+NE(I1) *X1 (I1) 200 YY=YY+NE(I1)*Y1(I1) CALL CALNEXI (NEXI,XI) XXI=O. YXI=O. DO 300 I1=1,3 XXI=XXI-X1(I1)*NEXI(I1) 300 YXI=YXI+Y1(Il)*NEXI(I1) DETJ=SQRT(XXI*XXI+YXI*YXI) XXI=XXI/DETJ YXI=YXI/DETJ DO 1000 II=1,NBIT ANG=(FIRST+FLOAT(II-1)*DANG)*3.141592654/180. CALL CALBE1(BE,ANG, XX,YY,XXI,YXI,WAVE,IPP) DO 400 I=1,3 JJ=NUM1(I,L) 400 FAR(II)=FAR(II) +WJ*NE(I)*BE*DETJ*f(jj) 1000 CONTINUE 700 CONTINUE 900 CONTINUE RETURN END 123

SUBROUTINE CALAE(.&E,X,Y,XP,YP,WAVE,IPP) COMPLEX AE,EPS,MU,HZERO,HONE,DXEPS,DYEPS PI=3.141592654 XK=2.*PI/WAVE XJI=X-XP YJI=Y-YP RJI=SQRT(XJI*XJI+YJI*YJI) RKO=XK*RJI CALL HANKZ2(RKO,2,HZERO,HONE) CALL EMFUN(XP,YP,EiPS,MU,IPP) CALL DRFUN(XP,YP,I)XEPS,DYEPS,IPP) AE=XK*XK*(MU-1./EPS)*HZERO*CMPLX(.0., O.25)-(DXEPS*XJI/RJI+ & DYEPS*YJI/RJI)/(EPS*EPS)*X*HONE*CMPLX(0.,O. 25) RETURN END SUBROUTINE CALAE1(AE, ANG,XP,YP,WAVE, IPP) COMPLEX AE,EPS,MU,HZERO,HONE,DXEPS,DYEPS PI=3.141592654 XK=2.*PI/WAVE RP=SQRT C(XP*XP+YP*YP) IF(RP.EQ.O.) GOTO 1 COP=XP/RP SIP=YP/RP COT=COS (ANG) SIT=SIN(ANG) HOLD=XK*RP*(COT*COP+SIT*SIP) HZERO=CMPLX(COS(HCILD),SIN(HOLD)) HONE=CMPLX(O.,l.) *HZERO CALL EMFUN(XP,YP,E;PS,MU,IPP) CALL DRFUN(XP,YP,DXEPS,DYEPS,IPP) AE=XK*XK*(MU-1./EFS)*HZERO*CMPLX(0., 0.25) + (DXEPS*COT+ t DYEPS*SIT)/(EPS*EF'S)*XK*HONE*CMPLX(0.,0.25) RETURN END 124

SUBROUTINE CALBE(BE,X,Y,XP,YP,XXI,YXI,WAVE,IPP) COMPLEX BE,EPS,MU,HZERO,HONE PI=3.141592654 XK=2.*PI/ WAVE XJI=X-XP YJI=Y-YP RJI=SQRT(XJI*XJI+'YJI*YJI) RKO=XK*RJI CALL HANKZ2(RKO,1,HZERO,HONE) CALL EMFUN(XP,YP,IEPS,NIJ,IPP) BE=(l.-1./EPS)*(X.JI/RJI*YXI+YJI/RJI*XXI)*HONE*XK*CMPLX(O.,O.25) RETURN END SUBROUTINE CALBEl ('BE,ANG,XP,YP,XXI,YXI,WAVE, IPP) COMPLEX BE, EPS,MU HZERO,HONE PI=3.141592654 XK=2. *PI/ WAVE RP=SQRT (XP*XP+YP*YP) COP=XP/RP SIP=YP/RP COT=COS CAIG) SIT=SIN (ANG) HOLD=XK*RP*(COT*CCOP+SIT*SIP) HZERO=CMPLX(COS (BOLD),SIN (BOLD)) HONE=CMPLX(O.,1.)*HZERO CALL EMFUN(XP,YP,E:PS,NM,IPP) BE=(l.-1l./EPS)*(CCOT*YXI+SIT*XXI)*HONE*XK*CMPLX(O.,O.25) RETURN END 125

SUBROUTINE CALNE (NE, XI,EA REAL NEW8 NE(1)=-0.25*Cl.-X:O)*(l.-ETA)*CXI+ETA+1.) NEC2)=0.25*(l.+XI')*C1.-ETA)*(XI-ETA-1.) NEC3)=0.25*(l.+XI,')*C1.+ETA)*(XI+ETA-1.) NE(4W=0.25*(l.-XI')*Cl.+ETA)*(-XI+ETA-1.) NE(5)=0.5*(l.-XI*]KI)*C1.-ETA) NE(6)=0.5*(1.+XI)*$(1.-ETA*ETA) NE (7)=0.5* (1.-XI*XKI)*(1l.+ETA) NE(8)=0.5*(l.-XI)*'(l.-ETA*ETA) RETURN END SUBROUTINE CALNEX-I(NEXI,XI,ETA) REAL NEXI(8 NEXI(1)=0.25*C1.-ETA)*(XI+ETA+1.)-0.25*(l.-XI)*(l.-ETA) NEXI(2)=0.25*(l.-ETA)*(XI-ETA-1.)+0.25*(l.+XI)*(l.-ETA) NEXI(3)=0.25*(l.+ETA)*(XI+ETA-1.)+0.25*Cl.+XI)*(l.+ETA) NEXI(4)=-0.25*(l.+~ETA)*(-XI+ETA-1.)-0.25*(l.-XI)*(l.+ETA) NEXI (5)=-XI* (1 -ETA) NEXI(6)=0. 5*(l.-ET'A*ETA) NEXI (7) =-XI* (1. +ETA) NEXI(8)=-O.5*(1l.-E,'TA*ETA) RETURN END SUBROUTINE CALNETA (NETA,XI,ETA) REAL NETA(8) NETAC1)=0.25*(l.-X.I)*(XI+ETA+1.)-0.25*(l.-XI)*(l.-ETA) NETA(2)=-0.25*(l.+,XI)*(XI-ETA-1.)-0.25*(l.+XI)*(l.-ETA) NETAC3)=0.25*(l.+XI)*(XI+ETA-1.)+O.25*(l.+XI)*(l.+ETA) NETA(4)=0.25*C1.-XI)*(-XI+ETA-1.)+0.25*(l.-XI)*(l.+ETA) NETA(5)=-0.5*(l.-XI*XI) NETAC6)=-ETA*(1l.+XI) NETA(7)=0.5*C1.-XI*XI) NETAC8)=-ETA*(1l.-XI) RETURN END 126

SUBROUTINE CALNEl (NE,XI) REAL NE(3) NE(1)=-0.5*ul.-XI:)*XI NEC(2)=0.5* (l. +XI)*kX NE(3)=l.-XI*XI RETURN END SUBROUTINE CALNEX.-l1(NEXI,XI) REAL NEXI(3) NEXI(1)=0.5*XI-O.!-*(l.-XI) NEXI(2)=0.5*XI+O.51*Cl.+XI) NEXI(3)=-2.*XT RETURN END 127