ELECTROMAGNETIC SCATTERING FROM LARGE PLANAR PLATES OF ARBITRARY SHAPE by Kasra Barkeshli, Arindam Chatterjee and John L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, Michigan 48109-2122 November 1989 389604-8-T = RL-2568

ABSTRACT ELECTROMAGNETIC SCATTERING FROM LARGE PLANAR PLATES OF ARBITRARY SHAPE by Kasra Barkeshli, Arindam Chatterjee and John L. Volakis The problem of electromagnetic scattering from large plana plates of arbitrary shape and resistivity is formulated using an improved conjugate gradient solution method. The method employs the discrete convolution theorem and Fast Fourier transform to achieve higher efficiency and accuracy over previously reported results. The accuracy of the technique is confirmed by comparison with measured data and patterns available based on alternative formulations for a variety of targets.

i TABLE OF CONTENTS LIST OF FIGURES............................ LIST OF TABLES................................ LIST OF APPENDICES............................ CHAPTER I. Introduction............................. II. Formulation............................... 2.1 Integral Equations........................ 2.2 Continuous Convolution....................... 2.3 Discrete Convolution....................... III. R esults.................................. 3.1 Comparison with measured data............... 3.2 Efficiency................................ IV. Conclusions and Future Work.................. * * 111 iv V 1 4 4 5 8 14 14 16 29 V. Users' M anual.................. 5.1 Input variables............... 5.1.1 Target geometry generation... 5.1.2 Resistivity taper specification. 5.1.3 Scattering pattern computations 5.2 Output................... 5.3 Sample problems............... 30....30.... 30.... 33....34.... 36.... 36 APPENDICES............. 40

ii BIBLIOGRAPHY............................... 63

iii LIST OF FIGURES Figure 2.1 Geometry of a polygonal plate illuminated by a plane wave....... 6 2.2 Evaluation of the self-cell element using approximate integration, four-term Taylor series expansion and the circular disk approximation 13 3.1 Backscattering pattern of a 2Ax2A square plate........... 18 3.2 E-polarization edge-on RCS for a rectangular plate......... 19 3.3 Convergence of the CGFFT solution for a square conducting plate. 20 3.4 Backscattering pattern of a circular plate (ka=5)........... 21 3.5 Backscatter pattern for an equilateral triangular plate of side 2A.. 22 3.6 Target geometry of the test plate................... 23 3.7 Azimuthal backscattering pattern of the test plate.......... 24 3.8 Effect of resistive tapering on the monostatic scattering of a 2x2 square plate............................... 25 3.9 3-Dimensional profile of a parabolic taper............... 26

iv LIST OF TABLES Table 3.1 Comparison between measured and calculated values at edge-on incidence (E-pol.) for a circular disk................... 27 3.2 Relative speedup between algorithms employing the Analytical and the Discrete transform of the Green's function............ 28 3.3 Speedup between vectorized and non-vectorized code......... 28

v LIST OF APPENDICES Appendix A. Spectral Domain Considerations..................... 41 A.1 Fourier Transform of the Current Density........... 41 A.2 Discrete Differentiation and its Transform........... 42 B. Listing.................................... 43

CHAPTER I Introduction The problem of electromagnetic scattering from large conducting and resistive plates is of considerable interest in electromagnetics as they constitute simple but nevertheless important components in man-made structures. Although they have been studied at a wide frequency range, experience with various numerical and asymptotic methods of solution as well as comparison with measured data reveals that there is a serious difficulty in predicting the scattering behavior of plates at grazing incidences where the edge currents and corner diffraction mechanisms are significant. Also, in the case of numerical methods, computer memory requirements play a major role in choosing the efficient method of solution. This report presents a technique to solve the problem of scattering from large planar plates iteratively using the conjugate gradient (CG) method. The conjugate gradient method is an iterative solution technique and has recently been found useful in electromagnetic applications ([1],[2],[3], [4],[5] and [6]). It is shown that the CG method in conjunction with the Fast Fourier Transform (FFT) is an efficient and in many instances the viable alternative to direct methods where the memory requirement for large targets becomes prohibitive. The CGFFT method was initially implemented by employing the analytical ex 1

2 pression for the Fourier transform of the free space Green's function to carry out the convolution integrals appearing in the pertinent integral equations [4],[7]. The elements of the dyadic Green's function were subsequently formed by carrying out the differential operations in the transform domain. This approach assumes an infinite domain for the Green's function as the Fourier transform is defined over the whole space (-oo, +oo). Thus-as far as the Green's function is concerned-the finiteness of the target's physical extent is not accounted for. As a result, unless a large FFT pad with extended zero elements is used, the method suffers from aliasing errors. This problem is more serious for the TM cases where the surface current density exhibits high spectral content due to edge singularity effects. It is the experience of the authors that a pad of at least three times the size of the target is often needed to get acceptable results for TM cases at oblique and close to grazing incidences. This makes the method inefficient in terms of both computer memory and speed particularly for monostatic computations since the solution must be repeated for each incidence. To alleviate this difficulty an alternative approach in handling the Green's function is taken resulting in improved accuracy and a substantial increase in speed. Here the problem is explicitly formulated in terms of discrete convolution of the current density with a discrete integral function of the Green's dyadic. In this case, the Fourier transform properly accounts for the finite domain of the Green's function, thus, avoiding aliasing. The piecewise constant expansion functions were used with point-matching although a provision for incorporating piecewise linear or sinusoidal (roof-top) basis functions has also been kept. An interactive Fortran code was developed based on this study which can handle a wide variety of targets with irregular symmetries. This is done by approximating

3 the target as a polygon with the coordinates of its corners specified by the user. This report is divided into two parts. The first part describes the formulation of the problem in terms of both continuous and discrete Green's functions. This is followed by an evaluation of the performance of the two formulations. The second part is a user's manual for the computer program including sample test cases for reference purposes.

CHAPTER II Formulation 2.1 Integral Equations Consider a thin inhomogeneous non-magnetic plate of resistivity Tq illuminated by a plane wave Ei as shown in Figure 2.1. The integral equation for the unknown surface current density J is given by El(r) = 7(r)J(r) + jkZ J(r') * r(r; r')ds' (2.1) where Z denotes the free space impedance, k = (2ir/A) is the free space wave number and r denotes the electric dyadic Green's function in unbounded space given by[4] r(r; r') = {I + k VV}G(r; r') (2.2) with e-jklr-r'l G(r; r') 4rlr- r' (2.3) in which r and r' denote the observation and integration points, both on the plate's surface. The explicit form of r is 1 02 1 02 ( k202&) k2XzOy r1 (1 — 2 G(r; r') (2.4) k2 ay0x k2 y2 4

5 Thus equation ( 2.1) represents a set of coupled integral equations: 1 a2 1 a2 E(x, y) = 7(x, y)J(x, y) - [(1 + + k )l + k y(2.5) k2 dI2 )', V2 aXay E 1 a2 1 a2 EY(, 7y) = 77(, Y)Jy(, Y)-[ n + (1 + +( )n,] (2.6) where II is the electric vector Hertz potential given by H(r) = -jkZ j J(r')G(r; r')ds'. (2.7) The integral in ( 2.1) is a convolution and can therefore be evaluated by invoking the convolution theorem. Introducing the spectral frequencies fx and fy and Fourier transform pair g= '{( } 8~, (2.8) with g(fx, fy) = L g(x, y)e-j27(f + fy)dxdy, (2.9) -,oo -00 (xy) = L (f=, f)ej2r(fx + f)Y)dfdf. (2.10) and employing the convolution theorem, the above system of equations can be formally written as Ei = r7J + jkZ.Fl{r J} (2.11) The above equation can now be efficiently solved via the conjugate gradient method. The description of the CG method and the algorithm used in this study are available in [1],[7]. 2.2 Continuous Convolution In a traditional implementation of the CGFFT method, the current density is expanded in terms of a subsectional surface basis function ' as M-1 N-1 J(x,y) = E E Jmn,mn(xy) (2.12) m=O n=O

z A / I x E. ll a O EH-pol. x ^ E. e Fig. 2.1. Geometry of a polygonal plate illuminated by a plane wave.

7 where mn(x, y) = ( -x (-x, y - yn), (2.13) -( (? Y) ~ ) (x,y)= ( (x y) 0, (2.14) 0 O~(x,y) and 4_, and y, are the expansion functions in the x and y directions, respectively. The Fourier transform of the current density vector is thus given by (Appendix I) J=J 3 (2.15) where J is the discrete Fourier transform of the sampled current. Similarly, the Fourier transform of the dyadic Green's function is expressed as (1 ) (-f G(f f) (2),16) f(-fzf,) (1 - f2) with e-jko/x' + y2 1 G(x,y)= 4rvx2 G(ff)= f f-1' (2.17) in which fc and fe are the spectral frequencies and use has been made of the following relationship for the partial derivative OX () j2,rfxf(f). (2.18) Equations (2.16) and ( 2.17) constitute analytical expressions for the free space dyadic Green's function. Substituting these into ( 2.11) and testing the resulting equation at discrete points (point-matching), the following expression is obtained Ei, = tiJ, + jkZ."l1{(r * ) j* } * (2.19) where the subscript ij denotes the value of the quantity at the point (xiyj) on the plate. It should be noted that to perform the Fourier transformation implied

8 by ( 2.11), an FFT pad twice the size of the plate must be employed. In general, however, a much larger pad is required when the analytical transform of the Green's function is used. To comply with a standard size FFT pad, we must choose N' = 2-, where - is any integer and N' > 2N - 1 > NNyquist with N reperesenting the number of elements/unknowns. Thus, 7 > log2(2N - 1)+ p (2.20) in which p is referred to as the order of the FFT pad. 2.3 Discrete Convolution The use of analytical expressions for the Fourier transform of the Green's function developed in the previous section assumes the entire space as the domain of the Green's function. Since the finiteness of the structure is not accounted for in this approach, aliasing errors are inevitable. For TM scattering, this necessitates the use of prohibitively large FFT pads in order to achieve reasonable accuracy. This approach also suffers from difficulties associated with the numerical handling of the known ring singularity of the transform of G1. By employing the discrete Fourier transform of the Green's function defined over an area twice the domain of the scatterer, aliasing is avoided (provided the required sampling criterion is satisfied). This leads to substantially more accurate results and faster convergence. More importantly, the results are insensitive to the FFT pad provided it is at least twice the plate size. In the new method, the integral equation is put into a discrete convolutional form for direct application of the discrete convolution theorem. Substituting for the 'The singularity corresponds to the circular ring f2 + f2 = 1 in equation ( 2.17).

9 current expansion in the integral on the right hand side of ( 2.1) yields [ mn * mn,,( y)] * f(, y; ', y')ds' (2.21) mn which upon interchanging the order of summation and integration, may be written as E Jmn * mn(x', y') r(x, y; x', y')ds' (2.22) Introducing ( 2.22) into ( 2.1) and applying point-matching yields the system of equations Eiij = rijJij + jkoZo Z Sij. Jmn (2.23) mn The dyadic function S = (2.24) is the new discrete kernel whose elements are given by = (1 I+y2 G(x, y; ', y')l, )(x' - xm,' - y )ds' =- - — 1 82 f r k3 = J2, G(x, y; x', y')y,(x' -zx, y -yY)ds' S? k2 9xay..,(,n = (i + - j G( x, y; x', y')( x' - xm, y '-y)ds' P ( 9X2 ) - Yn)ds| (2.25) where mnn is the incremental surface element corresponding to the mnth cell on the plate (x, -- a < x < xm +-,yn - -' < y < yn +- ). Obviously, the convolutional nature of the operation is preserved once the above functions are evaluated at the appropriate field points. Applying the discrete convolution theorem in ( 2.23) now yields Ei- =,iJij + jkoZoDFT-1{. J} (2.26)

10 In the above, 2 denotes the discrete Fourier transform of E and its accurate evaluation is a key to the success of the formulation. To calculate these elements, it is noted that in the discrete sense, the partial derivatives may be carried out by finite difference formulae. In particular, using a 3-point difference formula, we obtain the relation (see Appendix I) a O-g(x) 2 j2.rfSinc(rfAx)(f (2.27) which may be considered as the discrete counterpart to ( 2.18). Using ( 2.27) - may be expressed as ^ (1 - [2rfSinc(rfAzx)]2) - (2r)2ffySinc(rf)Sinc(fyAy) A -(2r)2ffyfSinc(irfyAy)Sinc(irf,4x) (1 - [27fySinc((rfy)y)]2) (2.28) where ~ is the discrete Fourier transform of the sequence e-j kc V7+' T mn n 4= 2 ds -Jk/TY W.m, 4./,-jx2 c + y2d = Lm+ / - 4/ ds- (2.29) J,,-qi $,-4l 4rV7/ + ya Note that in ( 2.27) and ( 2.28), Sinc(x) = sin(x)/x (2.30) It is noted that ( has an integrable singularity when ym = Yn = 0, which corresponds to the self-cell interaction. This term can be evaluated analytically using various approximations as stated below. Approximate integration: An approximate integration of the self-cell term yields[8] lI1 1 Ay oo {Axr lnltan( + tan-1 ) + a 2 A31 1 a AxAy +Aylnltan(4 + -tan )-jk 2 ' _. 4 2 Ay 2

11 Taylor series expansion: expanding the integrand in a four-term Taylor series expansion e-JkR 1( (j k r)2 (j kR)3 1 -(jR)3 (2.32) R -2 6 -k The integral can then be expressed as 1 k2 k3 AA M o (II - jk2l- I3 +i4 ) ' (2.33) 4w 2 b -2 2 where I = -ds = x ln(y + R) + yln(x + R) I2 = ds = AxAy IRd xyR x3 3 = j s = + ln(y + R) + Yln( + R) 3 6 4 = R2dS. Circular disk approximation: Using a circular disk approximation j2r fro -jkr 0oo --- rdrdo Jo Jo 47rr ro -jkro/2Sc( kro(2.34) nc(- (2.34) o = Figure 2.2 shows a comparison of these for the square cells (A =Ay = A) of different sizes. As seen, they all give values that are essentially indistinguishable for A < 0.1A.

12 The remaining terms 'mn are evaluated numerically using Guassian quadrature integration. The solution of ( 2.26) via the CG method will be referred to in the report as the CGDFT method.

Analytical Approximations to the Self-Cell Element 0.06 oS 0.04 0 X 0.02 0.00 Soemple Siz, A/A 0 v SJ 0 0 I4 Somple Size, A/A Fig.2.2 Evaluation of the self-cell element using approximate integration(-,) fourterm Taylor series expansion(. ), and circular disk approximationd-4.

CHAPTER III Results The above formulations were employed to compute the scattering by a number of plates. In this chapter, a sample of numerical results some of which are compared with available measured data or results based on alternative formulations are given. 3.1 Comparison with measured data Figure 3.1 shows the radar cross section for a rectangular conducting plate illuminated by E and H polarized plane waves. The edge-on behavior of a plate of constant width (b = 2A) and varying length (2A < a < 2.5A) is given in figure 3.2 with the electric field aligned with the shorter side. The results obtained are compared with measurement data reported by Hey and Senior[9]. It should be noted that the CGFFT method employing the continuous transform of the Green's function was found inadequate for an accurate prediction of the edge-on scattering behavior[4]. This is clearly demonstrated in Figure 3.3 where we compare the RCS of a square plate (of side 2A) as computed by the CGDFT method with p = 1. It is found to be in good agreement with the corresponding measured data. Since plates of arbitrary shape may be modeled as polygons, circular and triangular plates can also be treated just as easily. Figure 3.4 gives the radar cross 14

15 section of a circular plate of radius a with ka = 5. These results agree well with the Sommerfeld-Macdonald approximation for close-to-normal incidence angles. However at oblique incidences, this approximation breaks down. Table 3.1 provides a comparison with measured data[10] of the edge-on backscatter as computed with the CGDFT as a further testament of the method's accuracy at close-to-normal incidences. Also, in figure 3.5 backscatter RCS data for an equilateral triangular plate of side 2A are compared with the corresponding measured data. The agreement is clearly good considering a possible small inaccuracy in the measurements due to alignment difficulties. The scattering characteristics of geometrically complex targets may also be simulated by approximating the target by a polygon of n sides. This is illustrated in Figure 3.6 where the plate has been modeled as a polygon of 8 corners. The neargrazing (conical) azimuthal scattering pattern of this target is shown in Figure 3.7. Resistive plates are considered next. In practice conducting surfaces are replaced with resistive cards for cross section reduction purposes. By defining q as a function of position, different resistivity tapers may be treated. As an example, the resistivity can be expressed by a nonlinear function (x77(,y)=7c + (77 -c) [( -X2~ )I (y y (3.1) xi/2 + y '/2 (3. where tc and 7 may be considered as the resistivities at the center and the edges of a rectangular plate, respectively and Tr and r, denote the tapers in the corresponding directions. Figure 3.8 shows the effect of uniform and non-uniform (parabolic) resistive tapers on the monostatic cross section of a 2A x 2A plate. The parabolic resistivity profile is shown in figure 3.9.

16 3.2 Efficiency It was mentioned in the introduction that applying the traditional approach (using the analytical expression for the transform of the Green's function) requires a large FFT pad to reduce aliasing effects and characterize the weakly singular behavior of G in the vicinity of the ring f2 + fA = 1. Figure 3.3 illustrates the convergence of the far zone scattered field for a square plate as the size of the FFT pad is progressively increased. Also shown is the improved result using the method presented here. It is noted that a pad of at least 3 times the number of unknowns is needed to obtain acceptable results with the traditional method and although the general behavior of the backscattering cross section approaches that of the correct result, the convergence to the measured value is not clear near grazing incidence even with higher order pads. On the other hand, the new method gives a reasonable prediction using a FFT pad of order 1. In an earlier work [11], the vectorizable nature of the CGFFT algorithm was explored by identifying the major processes involved in a given iteration. Since in the CG algorithm a considerable amount of computation time is spent in operations which exhibit no data recurrence at a particular point in the iteration, these operations may be vectorized to increase the speed of calculations. This is particularly the case for the FFT which is a highly vectorizable algorithm and plays a major role in the speed and efficiency of the optimized code. In this study both scalar and optimized FFT routines (available on the IBM 3090's ESSL library) were used. Sample timing information are given in Tables 3.2 and 3.3. When the data from both of these tables are combined, one observes an overall speedup of more than 60 times over the scalar implementation using the continuous transform of the Green's

17 function. It should be noted that a vectorization of 97% is observed when employing the vector FFT routine.

18 Satterin from a erfectly Conducing Square Plate 30, E-POL Col*ated 0 E-POL Meured.20 - H-POL Caclite- O H-POL Mesuraed i Ai,,, 10: \/. 2,........, —... ~ 0 -10 Jo 0 0 10 20 30 40 50 60 70 80 90 Ang, of Icideme, dg. Figure 3.1: Backscattering pattern of a 2Ax2A square plate

19 E-Polarization Edge On RCS for a Rectangular Plate 10 a me *O 5 0 -5 -10 2 2.1 2.2 2.3 2.4 2.5 Side Length, a/X Figure 3.2: E-polarization edge-on RCS for a rectangular plate

20 Convergence of the CGFFT solution for a square conducting plate 2X x 2K Po de -O,< 0 so Lu I0 30 25 20 15 10 5 0 -5 -10 -15 -20 0 30 60 Angle of Incidece o,, dEg. 90 Figure 3.3: Convergence of the CGFFT solution for a square conducting plate

21 Backscattering pattern of a circular plate (ka=5) 20 -A H-POL. Sommerfeld-MacDonald Approx. ' E-POL. Sommerfeld-MacDonald Approx. 10 -0 '0 " Sommerfeld-MacDonald ' Approximation ^ ( m %u 40 I 0 30 60 90 -2 0 30Breakdown of60 90 Angle of Incidence (o, deg.

22 Backscattering pattern of a 2X2 triangular plate 30 25 20 CO 0 Vt,< 0 *r* () co 0 CO co 15 10 5 0 -5 -10 -15 -20 -25 -30 -90 -60 -30 0 30 60 Angle of Incidence too, deg. 90 Figure 3.5: Backscatter pattern for an equilateral triangular plate of side 2X

2 %. y< -- ' 2.5X Fig 3 T g o t t plate *::...*: *....... *...*..*....*..*....*: ** ** ** *~~ ~ ~ *....***.....**** **** **** *~~~~~~~~ *..... *.... * * * * * * * * * * * * * * * * ******:* *::* **:..... ******* ******* *.*:.:::*::*.**,* * ** *.... ****** ****** *~~~~~~~~ *..........**.:.:*** *** **** *** *** * * * ** * * ** * * **.:* * *: *....***** ** *** *..::..... ****** *** ********* *~~~:..................::::* *** *** *** *** *~~~~~~~~~~~~~~ ~ ~ ~ ~ ~ ~ ~ * *o *~~~~~~~~~~~~~ ~ ~ ~ ~ ~ ~ ~ * * *o * ~ ~ ~ ~ ~ ir 3.6 T*g *~et *f *h *~ *plate ******** ********

24 P3 a.o 0n 10 — 0 --10 — -20 — -30 --90 --- CGDFT (0.033k per element) I - DFTM/Moment Method (provided by Northrop Corp) ASPECT ANGLE (DEGREES) Figure 3.7: E-Polarization Conical (00 = 800) backscattering pattern of the test plate

25 Effect of Resistive Tapering on Monostatic Scattering of a Plate 2X x 2k 30 20 m 0 0 Uo C) cn uPo A: 10 0 -10 -20 -30 0 30 60 90 Angle of Incidence (pO, deg. Effect of resistive tapering on the monostaic scattering of a 2x2 X2 square plate Figure 3.8:

26 22 I I I X / / / / Figure 3.9: 3-Dimensional profile of a parabolic taper (Max R==-) 2

27 ka radius(a) Measured values Calculated values 1.0 0.1592 -7.696 -8.653 1.5 0.2387 -6.383 -6.274 2.7 0.4297 -5.229 -3.099 Table 3.1: Comparison between measured and calculated values at edge-on incidence (E-pol.) for a circular disk.

28 Green's function FFT Pad order Time (in ms) per iteration Relative speedup Discrete 1 65.49 20 1 65.62 19.96 Analytical 2 291.21 4.5 3 1310.62 1 Table 3.2: Relative speedup between algorithms employing the Analytical and the Discrete transform of the Green's function. Green's FFT Pad Time(in ms) per iteration Vector Function order Vectorized Non-vectorized Speedup Discrete 1 65.49 216.26 3.30 1 65.62 215.75 3.28 Analytical 2 291.21 840.34 2.89 3 1310.62 3597.12 2.75 Table 3.3: Speedup between vectorized and non-vectorized code.

CHAPTER IV Conclusions and Future Work The formulation discussed here may be easily extended to cylindrical structures with a minor modification of the dyadic Green's function and the incident field expressions. It should be noted that the problem is simplified for TM case where the longitudinal current component is the only unknown and only a single integral equation is to be solved. 29

CHAPTER V Users' Manual This chapter describes the significance of the various input variables in the CGFFT code and explains the procedure for running the program. It is the intention of the authors that the reader be fully acquainted with the CGFFT code after reading this manual. The program can be run in either interactive or non-interactive mode. The interactive mode is designed to help make the task of the inexperienced user easier. The interactive mode of running the program will be described first. 5.1 Input variables 5.1.1 Target geometry generation The CGFFT code is capable of performing scattering computations for targets of irregular symmetries (this will be shown later). The user can generate a target of his/her own choice by specifying the co-ordinates of the corners of the target. The program prompts for the target choice on the screen by printing the following message: **SCATTERING FROM THIN CONDUCTING PLATES** SPECIFY THE TARGET 30

31 1: RECTANGULAR PLATE 2: CIRCULAR DISC 3: TRIANGULAR PLATE 4: POLYGONAL PLATE 5: TEST CASE ENTER THE TARGET NUMBER: Option 5 includes a few special targets which were generated for testing purposes. It should be emphasized here that these targets could just as well be generated using the 'TARGET # 4: POLYGONAL PLATE' choice. On entering the target number, the program assigns the inputted value to the variable 'NTARG'. Inputting a number less than 1 or greater than 5 evokes the response 'WRONG TARGET! TRY AGAIN...'. If the target number specified was between 1 and 4, the program asks the user about the TARGET CENTER AND ENCLOSURE LIMITS: XO,YO,XL,YL The target center should be the point of symmetry for symmetrical targets and a convenient point (from which it is easier to calculate the corner co-ordinates) near the center point for asymmetrical ones. 'XOYO' are the input variables for the target center. The enclosure limits are the length and breadth of the target specified by the variables 'XL,YL'. For non-rectangular targets, the enclosure limits are specified by the dimensions of the largest rectangle that would just fit the geometry. The next message printed by the program is SAMPLE DENSITIES:UPLX,UPLY The number of samples per wavelength desired by the user can thus be specified. If the sampling is low, aliasing errors occur whereas oversampling results in high

32 memory requirement. The optimum number of samples per wavelength can thus only be decided by trial and error but the authors have found that 20-25 samples per wavelength (for large targets of sides 2 wavelengths long) give reliable results. The variables 'UPLX,UPLY' specify the sample densities. The target geometry specification is now complete except for the order of the FFT pad. FFT PAD OF ORDER: The discrete Fourier transforms are carried out by a radix 2 FFT routine. The sampling interval and the size of the FFT pads are chosen so that the Nyquist criterion in both the spatial and frequency domains as well as the requirements of linearity in the convolutions are satisfied. Thus, the period N' of the array to be transformed is chosen so that N'= 2: N' > NNyquist, N' > 2 x N- 1 (5.1) where N is the number of unknowns and 7 is an integer. In practice 7 is chosen according to the rule 7> log2(2N- ) + p (5.2) where p is an integer setting the order of the FFT pad dimension to ensure adequate frequency sampling in the spectral domain when performing the inverse transform operation. The array elements beyond the scatterer's physical range are set to zero. For computations using the discrete Green's Function, the pad should be of order 1 whereas for calculations with the analytical Green's Function, the pad order is limited only by the size of the memory and time considerations (the authors have found that a pad order exceeding 3 is too time-consuming without yielding a significant improvement in the results).

33 For target numbers 1, 3 and 4, the program asks for the coordinates of the cornerpoints: ENTER THE COORDINATES OF EACH CORNER SEQUENTIALLY: The only thing to remember here is that the co-ordinates of the points should be entered such that the points are arranged in either clockwise or counter-clockwise fashion. For target number 2 (a circle), the program asks the user for the radius of the circular target. For option number 5 (test cases), the geometry is automatically generated since the parameters are specified within the code. GEOMETRY DISPLAY? 1)YES 2)NO In order to make sure that the geometry generated is correct, it is recommended that the user check the geometry before plunging into the CGFFT calculations. Entering '1' in response to the above message, the geometry appears in the form of asterisks in a grid of dots. The entire target is actually modelled as a matrix where the asterisks indicate the target points('1') and the dots signify zeroes. The user may continue with the calculations or modify the geometry or exit from the program at this point depending on the response to the message: 1) CONTINUE 2) MODIFY GEOMETRY 3) EXIT 5.1.2 Resistivity taper specification The plate's resistivity may be specified externally or internally by setting the variable 'TPR' to a real value (TPR> -2.0). If a value of -1.0 is selected, the plate is assumed to be perfectly conducting and calculations involving resistivity are

34 omitted. A non-negative value for TPR denotes the taper in the x and y directions as indicated in equation ( 3.1). RESISTIVITY TAPER(REAL): -2.) RESISTIVITY VECTOR READ EXTERNALLY -1.) PERFECTLY CONDUCTING PLATE 0.) UNIFORM 1.) LINEAR 2.) PARABOLIC The variables 'ETC' and 'ETE' read in the resistivity values at the center and the edge of the plate, respectively. For non-symmetrical targets, it is advisable to read in the resistivity values externally. 5.1.3 Scattering pattern computations This section enables the user to specify the range of the 0 and < angles, the polarization angle b, the type of scattering calculations( monostatic or bistatic), the basis function, the Green's function (Analytical or Discrete) and the tolerance desired. TETA & PHI The program asks the user the initial 0 and f angles for which the scattering computations are to be done. The code accepts both positive and negative values for 0 and b. POLARIZATION ANGLE: 0.)H [HZ0O], 90.)E [EZ,0] If H-polarization scattering calculations are desired, 0. should be entered. This sets i to 0. For E-polarization computations, 4 should be set to 90.

35 1)BISTATIC 2)MONOSTATIC(BACKSCATTERING) COMPUTATIONS For bistatic calculations, 1 is to be inputted. Entering 2 causes the program to carry out monostatic (backscattering) computations. 1)ELEVATION 2)AZIMUTH CUT In order to find out the elevation pattern, the user should input 1. For azimuthal pattern calculation, 2 should be entered. INITIAL & FINAL ANGLES AND INCREMENT The user can substantially decrease the computation time by specifying the initial and final angles so as to exploit the symmetry inherent in the chosen target. A note of caution for the inexperienced user:even though the initial and final angles are real, the input for the angle increment expects an integer. BASIS: 1)CONVENTIONAL DFT 2)SURFACE PULSE 3)ROOF-TOP The user can employ either of the three basis functions for scattering computations. However, in this code the surface pulse basis has been only incorporated. 1)DISCRETE OR 2)ANALYTICAL TRANSFORM OF THE DYADIC G The user can employ either the Discrete or Analytical transform of the Dyadic Green's function. As discussed in Chapter 2, the analytical transform was employed in the original version of the code. It is present as an option but its use is not recommended due to its associated inefficiency. MAX & MIN NO. OF ITERATIONS:ITMAX,ITMIN & TOLERANCE

36 The user can specify the maximum and minimum number of iterations. However, it is recommended by the authors that the maximum number of iterations be kept to a reasonable value (say, 100) and the minimum number kept to 2. The tolerance value specified by the authors during running the program was.035 (3.5%). It is left to the user to tailor the values of 'ITMAX,ITMIN & TOLERANCE' according to his needs. 5.2 Output The information about the number of iterations and the tolerance achieved for each angle is displayed on the screen. In this way, the angles for which convergence was not achieved can be identified. The values for the scattering computations are stored in a data file in such a way that the data can be easily plotted. 5.3 Sample problems A sample problem is presented in this section to allow the user a better insight into the running of the CGFFT code. In the problem solved here, the backscattering pattern of a polygonal plate is computed using the discrete transform of the dyadic Green's function. **SCATTERING FROM THIN CONDUCTING PLATES** SPECIFY THE TARGET 1: RECTANGULAR PLATE 2: CIRCULAR DISC 3: TRIANGULAR PLATE 4: POLYGONAL PLATE

37 5: TEST CASE ENTER THE TARGET NUMBER: 4 TARGET CENTER AND ENCLOSURE LIMITS: XO,YO,XL,YL 0. 0. 2. 2. SAMPLE DENSITIES: UPLX,UPLY 25.5 25.5 FFT PAD OF ORDER: ENTER THE COORDINATES OF EACH CORNER SEQUENTIALLY: CORNER #1: -2.232 0.0 CORNER #2: -0.5 -1.0 CORNER #3: 0.5 -1.0 CORNER #4: 1.088 -0.809 CORNER #5: 1.451 -0.309

38 CORNER #6: 1.451 0.309 CORNER #7: 1.088 0.809 CORNER #8: 0.5 1.0 CORNER #9: -0.5 1.0 GEOMETRY DISPLAY? 1)YES 2)NO 2 1) CONTINUE 2) MODIFY GEOMETRY 3) EXIT 1 TETA & PHI 90. 90. POLARIZATION ANGLE: 0.)H [HZ0], 90.)E [EZ=0] 90. 1)BISTATIC 2)MONOSTATIC(BACKSCATTERING) COMPUTATIONS 2 1)ELEVATION 2)AZIMUTH CUT

39 1 INITIAL & FINAL ANGLES AND INCREMENT 90. 270. 3 BASIS: 1)CONVENTIONAL DFT 2)SURFACE PULSE 3)ROOF-TOP 2 1)DISCRETE OR 2)ANALYTICAL TRANSFORM OF THE DYADIC G 1 MAX & MIN NO. OF ITERATIONS:ITMAX,ITMIN & TOLERANCE 150 2.035 In the non-interactive mode of running the above program, the input values are stored in a file and the program reads from it. This mode is faster than the interactive mode.

APPENDICES 40

41 APPENDIX A Spectral Domain Considerations A.1 Fourier Transform of the Current Density The Fourier transform of the current density vector is defined as J(fx, f) = J(xy)e-j27r(fxx + fY)dxdy. (A.1) -00 Substituting for the current expansion ( 2.12)in the above, = Jmn mn(x,y)-j2(f + fY)dxdy J-oo 0- mn Jmn r r" mn(x,y)e-j2'(fzx + fyY)dxdy (A.2) mn J- J-oo and noting that?mn = J (x - x yx - yn)e-j2(fx + fyY)dxdy -00 -00 = e-j2r(fxm + fyYn) (A.3) equation (A.1) reduces to J = '(f, fy). j Emne-j2r(fsx + f y,) (A.4) mn with the summation on the right hand side defined as the discrete Fourier transform of J: m= n,,e-j2(fx + f. (A.5 mn

42 Thus, the desired relationship between the continuous and discrete Fourier transforms of the current density vector is J= J. (A.6) A.2 Discrete Differentiation and its Transform Numerical differentiation of a discrete function can be carried out by calculating finite differences at a number of points. A 3-point finite difference formula is given by - g - ) (A.7) with its transform ag,-Jf " - e 2i!S ax gA(fx) (A.8) which can be further reduced to Og -j2 -j2rfSinc('rfZAx)g. (A.9) More accurate expressions may be derived by using higher order difference formulae. For the 5-point case, ag 8g(x, + A) - g(x, + ) + g(x - ) -8g(xi - ) (A X(Xi) =2 (A.10) x 12Ax and the corresponding result is g -icic3 1 Oxg -ji2rf3[2Sinc(rf.Ax) - 1Sinc(3irf3Ax)]' (A.11) ax 3 4 9 A

43 APPENDIX B Listing This appendix provides a listing of the computer program PLTCGF developed and used in this study.

Print fil "pltcgf" C =========-= C DESCRIPTION OF THE PROGRAM C THIS PROGRAM COMPUTES THE CURRENT DENSITY AND THE RADAR SCATTERING C CROSS SECTION OF PLANAR PLATES OF ARBITRARY SHAPES MODELED BY C POLYGONS. ALSO, SEVERAL SPECIFIC GEOMETRIES ARE TREATED SEPARATELY C INCLUDING RECTANGLES, TRIANGLES AND CIRCLES. C C KASRA BARKESHLI C RADIATION LABORATORY C EECS DEPARTMENT C UNIVERSITY OF MICHIGAN C JULY 1989 C C COMPILATION PROCEDURE: C C....MTS (WITH VECTORIZATION) C C NOTICE: FOR VECTORIZATION ON THE IBM-3090, THE FOLLOWING STATEMENT C AND SIMILAR ONES IN THE SUBSEQUENT SUBROUTINES MUST BE UNCOMENTED C WITH THE '@' SYMBOL POSITIONED IN COLUMN 1. C ******************************************************************* CQPROCESS DIRECTIVE ('*VDIR' ) C PROGRAM PLTCGF C C $ RUN NEW:FORTRANVS SCARDS-PLTCGF PRINT-LIST SPUNCH-LOAD C PAR-OPT(3) VECTOR(LEVEL(2),REPORT) C $ RUN -LOAD+NAAS:NEW.ESSL+*FORTRANLIB C C....APOLLO WORKSTATIONS C C $ ftn pltcgf C $ bind pltcgf.bin -b pltcgf.obj C $ pltcgf.obj C C C PARAMETER DESCRIPTION C C C.. GEOMETRY C C NTARG....................... TARGET NUMBER C XO,YO...................... TARGET CENTER COORDINATES C XL,YL......................TARGET ENCLOSURE C UPLX,UPLY..................NUMBER OF UNKNOWNS PER WAVELENGTH IN C X AND Y DIRECTIONS C MX,MY.......................NUMBER OF UNKNOWNS IN X AND Y DIRECTI( C NX,NY......................FFT PAD SIZE IN X AND Y DIRECTIONS C IPAD.......................ORDER OF THE FFT PAD C ITAG......................INTEGER ARRAY CONTAINING THE ADDRESSE' C OF THE ACTUAL POINTS ON THE TARGET C... ELECTRICAL PROPERTY C C ETA................. RESISTIVITY C ITAP...................... TAPERING PARAMETER: C 0) UNIFORM RESISTIVITY C 1) LINEAR TAPER C 2) PARABOLIC TAPER C 3) RESISTIVITY VECTOR READ EXTERNALLY C C...PATTERN COMPUTATION C C TETA, PHI.................. SPHERICAL ANGLES OF INCIDENCE Page 1 ONS S

Print file "pltcgf" C PSI............ C C C ISCAT........... C C ICUT............ C ANGIANGFINC... C C... GREEN'S FUNCTION C C IG.............. C C C...EXPANSION FUNCTIONS C Page 2 *.. *ee *.. *........... POLARIZATION ANGLE: 0 )H POLARIZATION HZ=0 90)E POLARIZATION EZ=0.........PATTERN COMPUTATION FLAG: 1)BISTATIC 2)MONOSTATIC (BACKSCATTERING)......... NUMBER OF BISTATIC PATTERN CUTS......... INITIAL AND FINAL ANGLES AND INCREMENT............ METHOD OF GREENS FUNCTION EVALUATION: 1)DISCRETE TRANSFORM 2)ANALYTICAL TRANSFORM C IBASE...... C C C C C... ITERATION C C ITMAX,ITMINJ C C C...CGFFT VARIABLES C C E........... C JS.......... C R........... C P........... C A,AT....... C C ZT,W....... C C...TIMING C C BGTIMEINTIb C C................ TYPE OF BASIS FUNCTION FOR CURRENT EXPANSION: 1)CONVENTIONAL DFT 2)SURFACE PULSE 3)ROOF-TOP (NOT AVAILABLE), TOL............ MAXIMIM AND MINIMUM NUMBER OF ITERATIONS AND ITERATION TOLERANCET................. EXCITATION (INCIDENCE) FIELD.............. UNKNOWN CURRENT DENSITY.............. RESIDUAL VECTOR.............. SEARCH VECTOR..............THE INTEGRO-DIFFERENTIAL OPERATOR AND ITS FOURIER TRANSFORM.............. S P E C T R A L D O M A I N W O R K I N G A R E A ME,ENTIME.......BEGIN, INTERMEDIATE AND END TIMES PARAMETER(MXT-128,MYT-128,NXF-256,NYF"256) PARAMETER (MT"MXT*MYT, MTI=2*MXT*MYT, NTF=NXF*NYF, NTFI=2*NXF*NYF) INTEGER ITAG(MT) INTEGER*4 IT(3) REAL X(MT),Y(MT),RES(MT) COMPLEX GXX(NXF,NYF),GXY(NXF,NYF),GYX(NXF,NYF),GYY(NXF,NYF) COMPLEX JS(MTI),E(MTI),A(MTI),R(MTI),P(MTI) COMPLEX ZT(NTFI),AT(NTFI),W(NTF) COMPLEX XJ CHARACTER*2 FILE(2) COMMON /DIM/MX,MY,NX,NY,NT,MG,MGI,DX,DY,DS DATA TPI/6.28318530717959/,ZO/376.991118/,XJ/(O.,1.)/ DATA FILE/' jx','jy'/ ------------------------------------------------------------------.TARGET GEOMETRY SPECIFICATION WRITE (6,*) WRITE(6,*)'**SCATTERING FROM THIN CONDUCTING PLATES**' WRITE (6,*) WRITE(6,*)'SPECIFY THE TARGET:' WRITE (6,*) WRITE(6,*)'1: RECTANGULAR PLATE' WRITE(6,*)'2: CIRCULAR DISK' C C C C.. C 1

Print file "pltcgf" Page 3 WRITE(6,*)'3: TRIANGULAR PLATE' WRITE(6,*)'4: POLYGONAL PLATE' WRITE(6,*)'5: TEST TARGETS' WRITE(6,*) READ (5, *, ERR-1) NTARG IF ((NTARG.GE. 1).AND.(NTARG.LE.4))THEN WRITE(6,*)'TARGET CENTER AND ENCLOSURE LIMITS: XO,YO,XL,YL' READ (5, *, ERR-1) XO,YO, XL, YL ELSEIF (NTARG.EQ. 5)THEN WRITE(6, *)'TEST CASES' WRITE(6,*)'6: TARGET # 1' WRITE(6,*)'7: TARGET # 2' WRITE(6,*)'8: TARGET # 3' WRITE(6,*)'9: TARGET # 4' READ (5, *, ERR-1) NTARG END IF IF((NTARG.LT. 0).OR. (NTARG.GT.9))THEN WRITE(6,*) ' WRONG TARGET! TRY AGAIN...' GO TO 1 ENDIF WRITE (6, *) 'SAMPLE DENSITIES:UPLX, UPLY' READ (5, * )UPLX, UPLY WRITE(6,*)'FFT PAD OF ORDER:' READ (5, *) IPAD CALL GEOMTR(NTARG,XO, Y, XL, YL, UPLX, UPLY, IPAD, ITAG, X, Y) WRITE(6, *)'1) CONTINUE 2) MODIFY GEOMETRY 3) EXIT' READ(5, *) IRES IF(IRES.EQ.2)THEN GO TO 1 ELSEIF(IRES.EQ.3)THEN GO TO 600 END IF IF((MX.GT.MXT).OR. (MY.GT.MYT).OR. (NX.GT.NXF).OR. (NY.GT.NYF))THEN WRITE(6,2) WRITE (6, 3)MX, MY, NX, NY WRITE (6, 4)MXT,MYT,NXF,NYF GO TO 600 ENDIF 2 FORMAT(/,'**REQUIRED ARRAY DIMENSIONS EXCEED THOSE SPECIFIED.') 3 FORMAT(' REQUIRED: MXT-',1 3,' MYT-', 3,' NXF-', I3,' NYF-',I3) 4 FORMAT(' SPECIFIED:',4(5X, 3),/,'MODIFY DIMENSIONS & RERUN.') NT-NX*NY C C ---------------------------------- C C...RESISTIVITY PARAMETERS C WRITE(6, *)'RESISTIVITY TAPER(REAL):' WRITE(6,*)' -2.) RESISTIVITY VECTOR READ EXTERNALLY' WRITE(6,*)' -1.) PERFECTLY CONDUCTING PLATE' WRITE(6,*)' 0.) UNIFORM' WRITE(6, *) ' 1.) LINEAR' WRITE(6, *) ' 2.) PARABOLIC' READ (5, *)TPR INDX-0 IF(TPR.EQ.-2. ) THEN C C OPEN (11, FILE' rest. in' ) DO 2000 J-1,MY DO 1000 I-1,MX INDX-INDX+1 READ (11, *) RES (INDX) 1000 CONTINUE 2000 CONTINUE

Print file "pltcgf" CLOSE (11) ELSEIF(TPR.EQ. 0.)THEN C C...GET THE RESISTIVITY VALUE FOR THE UNIFORM CASE C WRITE(6,*)' ETA-' READ(5,*)ETA DO 4000 J-1,MY DO 3000 I-1,MX INDX-INDX+1 RES(INDX)-ETA 3000 CONTINUE 4000 CONTINUE ELSEIF(TPR.GT. 0.)THEN C C...GET THE RESISTIVITY TAPER AND SAVE C OPEN(10, FILE-' rest.out' ) WRITE(10, *) 'RESISTIVITY TAPER' WRITE (10, *)MXMY WRITE(10,*) (I, I-1,MX) WRITE(10,*) (J J-1,MY) WRITE(6, *) 'ETA(CENTER), ETA(EDGE)' READ (5, *) ETC,ETE WIDTHX-XL/2. WIDTHY-YL/2. DO 6000 J-1,MY DO 5000 I-1,MX INDX-INDX+1 AC1-ABS (X(INDX) -XO) AC2-ABS (Y(INDX) -YO) RES (INDX)+ ETC+ (ETE-ETC) * ((AC1/WIDTHX) **TPR+ (AC2/WIDTHY) **TPR) WRITE(10, *)RES(INDX) 5000 CONTINUE 6000 CONTINUE CLOSE ( 10) ENDIF C C --------------------------------------------- C C...PARAMETERS FOR SCATTERING PATTERN COMPUTATIONS C WRITE(6,*)'TETA & PHI' READ (5, *) TETA, PHI WRITE(6,*)'POLARIZATION ANGLE: 0.)H [HZ-0], 90.)E [EZ-0]' READ(5, *)PSI WRITE (6, *) ' 1) BISTATIC 2) MONOSTATIC (BACKSCATTERING) COMPUTATIONS' READ(5,*)ISCAT WRITE(6,*) '1)ELEVATION 2)AZIMUTH CUT' READ(5, *)ICUT WRITE(6,*)'INITIAL & FINAL ANGLES AND INCREMENT' READ(5, *)ANGI,ANGF, INC C C ---- - -------------------------------------------------- C C...THE SURFACE BASIS FUNCTIONS TO BE INCORPORATED C WRITE(6,*)'BASIS: 1)CONVENTIONAL DFT 2)SURFACE PULSE 3)ROOF-TOP' READ (5,*)IBASE C C --- —----------------------------- C C...THE METHOD OF GREENS FUNCTION EVALUATION C WRITE(6,*) '1)DISCRETE OR 2)ANALYTICAL TRANSFORM OF THE DYADIC G' Page 4

Print file "pltcgf " Page 5 READ (5, *) IG CALL GREENS(GXX, GXY,GYX, GYY, IG, IBASE) C C ------------------------------------------ C C.. ITERATION PARAMETERS C WRITE(6,*)'MAX & MIN NO. OF ITERATIONS:ITMAX,ITMIN & TOLERANCE' READ (5,) ITMAX, ITMIN, TOL C C ------------------------------------------ C C...INITIAL GUESS (OF ZERO) C DO 6 I - 1,MGI JS(I) - (0.,0.) 6 CONTINUE OPEN (UNIT-3, FILE-' sigma0' ) WRITE(3,7) 7 FORMAT (' INPUT DATA.',, ' "SIGMAO" ') 10 CONTINUE WRITE (6, 15) TETA, PHI 15 FORMAT('ANGLES OF INCIDENCE: TETA-',F5.1,', PHI-',F5.1) ITER - 0 C C ----- --------- -------------------------------------- C C...INCIDENT FIELD AND ITS NORM C CALL INCIDE(E,TETA,PHI, PSI, X,Y) ENORM - VNORMX(E) C C ------------------- -------------- a --- —------ C C...CGFFT ZEROTH ORDER RESIDUAL CALCULATION C CALL OPERAT(A, -1, JS,TPR RES ZT,GXX, GXY,GYX,GYY, AT, W,ITAG) DO 20 I - 1,MGI R(I) - A(I)-E(I) P(I) - (.,0.) 20 CONTINUE 30 CONTINUE C C --- -------------------------------------------------------- C C RESET TIME C CALL TIME(0,PR,RES,RC4,RC8) C C...CGFFT MAIN ITERATION LOOP C CALL OPERAT (A, 1,R,TPR, RES, ZT, GXX, GXY, GYX, GYY, AT,W, ITAG) B - 1./VNORMX(A) DO 40 I - 1,MGI P(I) - P(I)-B*A(I) 40 CONTINUE CALL OPERAT(A, -1,P,TPR, RES, ZT, GXX,GXY, GYX, GYY, AT, W, ITAG) T - 1./VNORMX(A) DO 50 I - 1,MGI R(I) - R(I)+T*A(I) JS(I) - JS(I)+T*P(I) 50 CONTINUE C C -------- --------------------------------------------------- C C... TOLERANCE CHECK

Print file "pltcgf" Page 6 C ITER - ITER+1 C COMPUTE RELATIVE RESIDUAL ERROR RSS - SQRT(VNORMX(R)/ENORM) IF((RSS.LE.TOL).AND.(ITER.GT.ITMIN)) THEN WRITE (6,400) ELSE IF(ITER.GE.ITMAX) THEN WRITE(6,500) ELSE GO TO 30 END IF 60 FORMAT (5 (1X, G12.5)) C C CPU TIME FOR THIS ANGLE OF INCIDENCE C CALL TIME(26,PR,RES,RC4,RC8) C C ------------- ---------------------------------------- C C... SCATTERING CALCULATIONS C IF(ISCAT.EQ.2)THEN CALL SCATER(SIG,JS,TETA,PHI,X,Y) IF(ICUT.EQ.1)THEN WRITE (6,180) TETA, ITER,RSS, SIG WRITE (3,200) TETA, SIG TETA-TETA+INC IF(TETA.LE.ANGF) THEN GO TO 10 END IF ELSEIF (ICUT.EQ.2) THEN WRITE (6, 180) PHI, ITER, RSS, SIG WRITE (3,200) PHI, SIG PHI-PHI+INC IF(PHI.LE.ANGF) THEN GO TO 10 END IF END IF ELSEIF(ISCAT.EQ.1)THEN WRITE (6, *)ITER,RSS C C...FOR BISTATIC CASE, SAVE THE CURRENT DENSITY C DO 110 I-1,2 L1-(I-1)*MG OPEN(I,FILE-FILE(I)) WRITE (I,*) 'TITLE' WRITE(I, *)MX,MY WRITE(I,*) (J,J-1,MX) WRITE(I,*)(K,K-1,MY) DO 70 L-1,NT W(L)-(0., 0.) 70 CONTINUE DO 80 L - 1,MG W(ITAG(L)) - JS(L+L1) 80 CONTINUE DO 100 K-1,MY L2-(K-1)*NX DO 90 J-1,MX WRITE (I, *) CABS (W(L2+J)) 90 CONTINUE 100 CONTINUE CLOSE(I) 110 CONTINUE IANGI-ANGI IANGF-ANGF

Print file "pltcgf" Page 7 IF(ICUT.EQ.1)THEN DO 150 I-IANGI,IANGF,INC TETA-I CALL SCATER(SIG,JS,TETA,PHI,X,Y) WRITE (6, 180) TETA, SIG WRITE (3,200) TETA, SIG 150 CONTINUE ELSEIF(ICUT.EQ.2)THEN DO 160 I-IANGI,IANGF,INC PHI=I CALL SCATER(SIG,JS,TETA,PHI, X,Y) WRITE (6,180)PHI, SIG WRITE(3,200)PHI, SIG 160 CONTINUE END IF END IF WRITE(3,*)'END OF DATA.' CLOSE(3) CLOSE (9) WRITE(6,*)'SCATTERING COMPUTATIONS COMPLETED' WRITE(6,*)'ANOTHER TARGET? 1) YES 2)NO' READ(5, *)IRES IE(IRES.EQ.1)GO TO 1 WRITE (6, *) 'NORMAL TERMINATION. ' 180 FORMAT(' ANGLE-', 1X, F5.1,2X, 'SIGMAO',F9.3) 200 FORMAT (F5.1,2X, F9.3) 300 FORMAT(I4, 1X,G13.4) 400 FORMAT(' CONVERGENCE ACHIEVED.') 500 FORMAT(' ITMAX EXCEEDED; CONVERGENCE CRITERION NOT MET.') 600 STOP END C * TARGET GEOMETRY SPECIFICATIONS C ininini B"-IC "=55Z SUBROUTINE GEOMTR(NTARG,XO,YO,XL,YL,UPLX,UPLY, IPAD, ITAG,XT,YT) INTEGER ITAG(*) REAL XT(*),YT(*),XC(10),YC(10) CHARACTER*1 ITARG(256,256) LOGICAL TAGTRITAG COMMON /DIM/MX,MY, NX,NY,NT,MG,MGI,DX,DY,DS TAG-. FALSE. IF (NTARG.EQ.1)THEN C C... RECTANGULAR PLATE C ELSEIF (NTARG.EQ.2)THEN C C... CIRCULAR DISK C WRITE(6,*)' ENTER THE RADIUS:' READ(5, *)RAD ELSEIF ((NTARG.EQ. 3).OR. (NTARG. EQ.4) ) THEN IF (NTARG. EQ. 3) THEN C C...TRIANGULAR PLATE C NC-3 ELSEIF (NTARG.EQ.4)THEN C C...POLYGONAL PLATE C WRITE(6,*)'NUMBER OF CORNERS:' READ *,NC END IF NT-NC-2 WRITE(6,*)'ENTER THE COORDINATES OF EACH CORNER SEQUENTIALLY:'

Print file "'pltcgf " Page 8 DO 30 I-1,NC WRITE(6,*)'CORNER #',I READ (5, *) XC (I), YC (I) 30 CONTINUE ELSEIF (NTARG.EQ. 6) THEN XL=3. YL=2. X0-l.5 Y0-0. C1-SQRT (3.) ELSEIF (NTARG.EQ.7) THEN XL=4. YL=2. XO=2. Y0-0. Cl=SQRT (3.) C2-1.+C1 ELSEIF (NTARG. EQ. 8) THEN XL-3.5 YL-2. XO-1.75 YO-0. C1-2.5 ELSEIF (NTARG.EQ. 9) THEN XL-2.5 YL-2. XO-.25 Y0-0. CI-SQRT (3.) C2=C1/2. END IF XL2-XL/2. YL2-YL/2. XSTRT-XO-XL2 YSTRT"YO0-YL2 MX-XL*UPLX MY"YL*UPLY NX-.2** (INT (ALOG (2. *MX) /ALOG (2.) ) +IPAD) NY="2** (INT (ALOG (2. *MY) /ALOG (2.) ) +IPAD) DX=1./UPLX DY-1./UPLY DO 20 J-",NY DO 10 I-1,NX ITARG(I, J)-'.' 10 CONTINUE 20 CONTINUE INDX-0 DO 60 J=",MY Y-YSTRT+ (J-0.5) *DY Ll""(J-l) *NX DO 50 I-",MX X-XSTRT+ (I-0.5) *DX L-Li+I C C C C...TARGET # 1 RECTANGLE C IF(NTARG.EQ.1)THEN TAG —.TRUE. C C... TARGET # 2 CIRCULAR DISK C ELSEIF (NTARG.EQ.2) THEN IF( (SQRT((X-X0)**2+(Y-Y0)**2))'.LT.RAD )TAG=.TRUE.

Print file "pltcgf" Page 9 C C ----------------------------------------------- C C...TARGETS # 3 AND 4 TRIANGLE AND POLYGON C ELSEIF( (NTARG.EQ.3).OR. (NTARG.EQ.4))THEN DO 40 K=1,NT IF( + TRITAG (X, Y, + XC(1),YC(1),XC(K+1),YC(K+1),XC(K+2),YC(K+2)) + )TAG-.TRUE. 40 CONTINUE C C ------------------------------------------------------------------ C C...TARGET # 6 C ELSEIF (NTARG.EQ. 6) THEN IF( + (X.LE.C1).AND. (ABS (Y).LE.X/C1) +.OR. + (X.GE.C1).AND. ( ((X-C1) **2+Y**2).LE.1.) + )TAG-.TRUE. C C ------------------------------------------- C C...TARGET # 7 C ELSEIF (NTARG. EQ. 7)THEN IF( (X.LE.C1).AND. (ABS (Y).LE.X/C1) +.OR. + (X.GE.C1).AND.(X.LE.C2).AND.(ABS(Y).LE.1.) +.OR. + (X.GE.C2).AND.(((X-C2)**2+Y**2).LE.1.) )TAG-.TRUE. C C ------------------------------------------------- C C...TARGET # 8 C ELSEIF (NTARG.EQ. 8) THEN IF( (X.LE.C1).AND. (ABS(Y).LE.1.) +.OR. + (X.GE.C1).AND. ( ((X-C1) **2+Y**2).LE.1.) )TAG=.TRUE. C C --------------------------------------------------- C C...TARGET # 9 C ELSEIF (NTARG.EQ. 9)THEN IF( ((ABS(Y).LE.C2).AND. (Y.GE.C1* (ABS(X)-.5))) +.OR. + ((X.GE. 0.).AND.(Y.LE.C2) +.AND.(((X-.5)**2+Y**2).LE.1.)) )THEN TAG-. TRUE. IF ( (X. LE. O.).AND. (ABS(Y).LE.DY/2. ))THEN TAG-. FALSE. ENDIF ENDIF END IF C C.. TAG ASSIGNMENT C IF (TAG) THEN ITARG(I,J)-' *' INDX-INDX+1 ITAG(INDX) -L

Print file "Pltcgf" ag 1 Page I 0 XT (INDX) -X YT (INDX) =Y TAG-.FALSE. END IF CONT INUE CONT INUE 50 60 C C C C..C 70 80 C C C C C C C.TOTAL NUMBER OF 'NONZERO' ELEMENTS WRITE (6, * WRITE(6,*)'XI= 'f,XSTRT+0.5*DX,' XF= ',X,' MX= ',MX,' WRITE(6,*)'YI= ',YSTRT+0.5*DYrf YF= ',Y,' MY= ',MY,' WRITE(6, *) MG-INDX MGI=2*MG WRITE (6,f *) 'GEOMETRY DISPLAY? 1) YES 2) NO' READ (5, *)IRES IF (IRES. EQ. 1) THEN DO 70 J=NYl,-1 WRITE(6,80) (ITARG(I,,J),,I=1,NX) CONTINUE FORMAT (128 (1XAl)) END IF WRITE (6, *) WRITE (6, *)'TARGET GEOMETRY SPECIFICATION COMPLETED.' RETURN END NX= ',rNX NY= I',NY * IDENTIFYING A POINT IN THE INTERIOR OF A TRIANGLE SPECIFIED BY * * ITS CORNERS* LOGICAL FUNCTION TRITAG(X,YX1,,YlX2,Y2,X3,Y3) TRITAG-. FALSE. Zl = ( X - Xl) * ( Y2-Yl ) - ( X2-Xl ) * ( Y Z2 - ( X3 - Xl) * ( Y2-Yl ) - ( X2-X1 ) * ( Y3 Z3 - ( X - Xl) * ( Y3-Yl ) - ( X3-Xl ) * ( Y Z4 - ( X2 -Xl) * ( Y3-Yl ) - ( X3-Xl ) * ( Y2 Z5 - ( X - X2) * ( Y3-Y2 ) - ( X3-X2 ) * ( Y Z6 - ( Xl - X2) * ( Y3-Y 2) - ( X3-X2 ) * ( YI. IF( -Y1 -Y1 -Y1 -Y1 -Y2 -Y2 ) + (zl*z2.GE. 0.0).AND.(Z3*Z4.GE. 0.0).AND.(Z5*Z6.GE. 0.0) + )TRITAG-.TRUE. RETURN END min m n m mm min in mnw= * EVALUATION OF DYADIC GREEN'S FUNCTION AND ITS FOURIER TRANSFORM * SUBROUTINE GREENS (GXXGXYGYXGYY, IG, IBASE) REAL*8 AUXi1(20000),AUX2 (20000) REAL*8 DX2DDY2DDFXDFYFXFYDERXDERYCXXCXYCYYDBSINC REAL* 8 TP IDrTPI2D, S COMPLEX XJ, XJK, MYCSQR, FCT COMPLEX GI, GOGX, GY, GOSELF, GONSLF COMPLEX GXX(NXNY),GXY(NXNY),GYX(NXNY),GYY(NXNY) INTEGER N(2) COMMON /DIM/MXMYNXNYNTMGMGIDXDYDS DATA PI/3.141592653589793/,TPI/6.28318530717959/ DATA TPI2/39.47841760435747/,TPI2D/39.47841760435747D0/ DATA FPI/12.56637061435917/,ZO/376.991118/ DATA TPID/6.28318530717959D0/ DATA XJ/(0.rl.)/,XJK/(0.6.28318530717959)/ FCT-XJ*ZO fTP I/NT DS=DX*DY DX2-DX/2.

Print file F'pltcgfv DY2-DY/2. DX2D-DBLE (DX2) DY2D-DBLE(DY2) DFX1.DO/DBLE (DX*NX) DFYil.DO/DBLE(DY*NY) NXH = NX/2+1 NYH = NY/2+1 C DO 2 J=1,NY DO 1 I-1,NX GXX(IJ)-(0., 0.) GXY(IJ)-(O, 0.) GYX(I,J)-(O.,O.) GYY(IJ)=(O,0.) 1 CONTINUE 2 CONTINUE C C C C... DISCRETE GREEN'S FUNCTION C IF(IG.EQ. 1)THEN N(1) =NX N(2)=NY RHO-SQRT (DS/PI) RBOUND-RHO DO 20 J=-MY+1,MY-1 IF(J.LT.0O)THEN JQ=NY+J ELSE JQ-J END IF ETA-J*DY ETAP-ETA+DY2 ETAM-ETA-DY2 DO 10 Ii-MX+1,MX-1 IF(I.LT O)THEN IP-NX+I ELSE IP=I END IF PSI-I*DX PSIP-PSI+DX2 PSIM-PSI-DX2 R-SQRT (PSI**2+ETA**2) IF(R.GT.RBOUND) THEN GO - GONSLF(R) ELSE GO m GOSELF(PSIPETAP) + GOSELF(PSIMETAM) + - ( GOSELF(PSIPETAM) + GOSELF(PSIMETAP) END IF GXX(IP+1,JQ+1) - GO I FPI GXX(IP+1, JQ+1) -FCT*GXX(IP+1, JQ+1) 10 CONTINUE 20 CONTINUE C CALL SCFT2(1,GXX,1,NXGXX,1,NXNXNY,1,1.O,AUXl,20000,AUX2,f20COO) C CALL SCFT2(OGXX,1,NXGXX,1,NXNXf,NY,1,1.O,AUX1,20000,AUX2,2000) CALL FFTNDM(GXXN,2, -1) DO 40 J-1,NY IF(J.LE.NYH)FY-(J-1) *DFY IF(J.GT.NYH)FY-(J-NY-1) *DFY DERY-FY*DBS INC (TPID*FY*DY2D) CYY-1. DOIDERY*DERY DO 30 I-1lNX IF(I. LE. NXH)FX=(I-1) *DFX IF(I.GT.NXH)FX-(I-NX-1) *DFX Page 11

Print fleo "'pltcgf" -Page 12 DERX-FX*DBS INC (TP ID *FX*DX2D) CXX-1. DO DERX*DERX CXY=-DERX*DERY GYY(IJ)=GXX(IJ) *SNGL(CYY) GXY(IrJ) -GXX(IrJ) *SNGL (CXY) GYX(I, J) -GXY(I, J) GXX(I, J) -GXX(I, J) *SNGL (CXX) 30 CONT INUE 40 CONTINUE C C C...ANALYTICAI TRANSFORM OF THE GREEN'S FUNCTION C ELSEIF(IG.EQ.2)THEN DO 60 J-1,NY IF(J. LE.NYH) FYI(J-1) *DFY IF(J.GT.NYH)FY-(J-NY-1) *DFY CYY-TPI2D*(1.DO-FY*FY) DO 50 I-1,NX IF(I. LE. NXH)FX-(I-1) *DFX IF(I.GT.NXH)FXu.(I-NX 1)*DFX CXX-TPI2D* (1.DO-FX*FX) CXY-TP I2D*FX*FY SIFX*FX+FY*FY-l.DO GI-2. *TPI*MYCSQR(S) IF(CABS(GI).LT.1.E-20)GI-l.E-10 GXX(I, J) =FCT*SNGL (CXX) /GI GXY(I, J) =FCT*SNGL (CXY) /GI GYX(I,J)=GXY(IJ) GYY(I, J) FCT*SNGL (CYY) /GI 50 CONTINUE 60 CONTINTUE END IF C C C... BASIS FUNCTIONS C DO 80 J-1,NY IF(J. LE. NYH) FY-(J-1) *DFY IF(J.GT.NYH)FY-(J-NY-1) *DFY DO 70 I-1,NX IF (I. LE. NXH)FX- (I1-4) *DFX IF (I.GT.NXH) FX- (I-NX-1) *DFX IF(IBASE.EQ.1)THEN ELSEIF(IBASE.EQ.2)THEN GXX (I, J) -GXX (I, J) *PWCNST (DX, FX) *PWCNST (DY, FY) GXY(IJ)-GXY(IJ) *PWCNST (DXFX) *PWCNST (DYFY) GYX(I, J) -GYX(I, J) *PWCNST (DX, FX) *PWCNST (DY, FY) GYY(I, J) -GYY(IrJ) *PWCNST (DX, FX) *PWCNST (DY, FY) ELSEIF(IBASE.EQ.3)THEN GXX (IfJ) -GXX (II, J) *PWSINE (DX, FX) *PWCNST (DY, FY) GXY(ItJ) -GXY(ItJ) *PWCNST (DX, FX) *PWSINE (DY, FY) GYX (I, J) -GYX (I, J) *PWS INE (DX, FX) *PWCNST (DY, FY) GYY (I,tJ) -GYY (ItJ) *PWCNST (DX, FX) *PWS INE (DY, FY) END IF 70 CONTINUE 80 CONTINUE RETURN END C C C GONSLF COMPLEX FUNCTION GONSLF(R) COMPLEX GFREESXJK

Print file "pltcgf" COMMON /DIM/MX,MY,NXNYNTMGMGIDXDYDS DATA TPI2/39.47841760435747/,XJK/(O.,6.28318530717959)/ GFREES (R) CEXP (-XJK*R) /R GONSLF-TPI2*DS*GFREES (R) RETURN END Page 13 C C C *C GOSELF * COMPLEX FUNCTION GOSELF(XY) COMPLEX XJK COMMON /DIM/MXMYNXNYNTMGMGIDXDYDS DATA TPI2/39.47841760435747/,XJK/(O.,6.28318530717959)/ X2-X*X Y2-Y*Y X3-X*X2 Y3 -Y* L~Y2 XY-X *Y R2-X2 +Y2 R-SQRT(R2) ALX-ALOG (X+R) ALY-ALOG (Y+R) Gl=X *ALY+Y*ALX G2-DS G3-XY*R/3. +(X3*~ALY+Y3*'ALX) /6. G4-XY*R2 /3. GO SELF-TP 12*+(Gl-XJK*JG2 -TPI2/2.*.G3+XJTK*TPI2/6.*. G4) RE TURN END TRANSFORM OF THE SURFACE PIECEWISE CONSTANT EXPANSION FUNCTION " REAL FUNCTION- PWCNST(DjF) DATA TPI/6 d.2883185307O 717959 / C-TPI*F" PWCNST-SNS INC (C*D)/2.) RET~URN END C C C C C C = TRANSFORM OF THE OVERLAPPING SINUSOIDAL EXPANSION FUNCTION REAL FUNCTION PWSINE(DF) DATA TPI/6.28318530717959/ C=TPI*F C1=TPI*TPI*(1.-F*F) IF(C.EQ. 0.)THEN PWSINE-1. ELSE PWSINE-2.*eTPI (COS (C*D) -COS (TPI*D) ) / (D*C1*SIN (TPI*D)) END IF RETURN END *t INCIDENT FIELD C C C C(PROCESS DIRECTIVE ('*VDIR') SUBROUTINE INCIDE(ETETAPHI, PSIXY) REAL X(*),Y(*) COMPLEX E(*),XJPHASE COMMON /DIM/MXMYNXNYNTMGMGIDXDYDS DATA PI/3.141592653589793/,TPI/6.28318530717959/,zO/376.991118/ DATA XJ/(O.,1.)/ STE - SIN(TETA*PI/18O.) CTE - COS(TETA*PI/18O.) SPH - SIN(PHI*PI/180.) CPH - COS(PHI*PI/18O.)

Print file "pltcgf " SPS = SIN(PSI*PI/180.) CPS = COS(PSI*PI/180.) RX = TPI*STE*CPH RY = TPI*STE*SPH EO=Z0 EX = EO * ( CPS*CTE*CPH - SPS*SPH ) EY = EO * ( CPS*CTE*SPH + SPS*CPH ) C*VDIR IGNORE RECRDEPS DO 10 I-1,MG PHASE = CEXP ( XJ * ( RX*X(I) + RY*Y(I)) ) E(I) - EX * PHASE E(I+MG) - EY * PHASE 10 CONTINUE RETURN END Page 14 C C C * EUCLEDEAN NORM CALCULATION * REAL FUNCTION VNORMX(W) REAL*8 SUM,AMAG COMPLEX W(*) COMMON /DIM/MX, MY, NX, NY, NT,MG, MGI, DX, DY, DS SUM - 0.DO DO 10 I - 1,MGI AMAG - DBLE ( CABS ( W(I) ) ) SUM - SUM+AMAG*AMAG 10 CONTINUE VNORMX - SNGL(SUM) RETURN END C C * INTEGRO-DIFFERENTIAL OPERATOR INVOLVING CONVOLUTIONS CARRIED OUT * C * IN THE TRANSFORM DOMAIN * C C@PROCESS DIRECTIVE (' *VDIR' ) SUBROUTINE + OPERAT (A, ISIGN, Z,TPR, RES, ZT,GXX, GXY, GYX, GYY, AT,W, ITAG) INTEGER ITAG(*) REAL RES(*) COMPLEX XJ,A(*),Z(*),ZT(*),AT(*),W(*) COMPLEX GXX(NX,NY),GXY(NX,NY),GYX(NX,NY),GYY(NX,NY) COMMON /DIM/MX,MY,NX,NY,NT,MG,MGI,DX,DY,DS,MT DATA PI/3.141592653589793/,TPI/6.2831853071796DO/,ZO/376.991118/ DATA XJ/(0.,1.)/ C C ------------------------------------- C C...FORWARD FOURIER TRANSFORM C CALL SPECTR (Z, ZT, -1, W, ITAG) C I 0 C --------------------------------------------------------- C C... REGULAR OPERATOR C IF(ISIGN.EQ.-1)THEN C*VDIR IGNORE RECRDEPS DO 20 J-1,NY L=(J-1) *NX DO 10 I-1,NX K-L+I AT (K) -GXX (I, J) *ZT (K) +GXY (I, J) *ZT (K+NT) 10 CONTINUE 20 CONTINUE C*VDIR IGNORE RECRDEPS DO 40 J=1,NY

Print file "pltcgf" Page 15 L=NT+(J-1)*NX DO 30 I-1,NX K-L+I AT (K) =GYX (I, J) *ZT (K-NT) +GYY (I,J) *ZT (K) 30 CONTINUE 40 CONTINUE C C...... --- —------—. --- —---------------------- C C..ADJOINT OPERATOR C ELSEIF(ISIGN.EQ. 1)THEN C*VDIR IGNORE RECRDEPS DO 60 J-1,NY L=(J-1) *NX DO 50 I-1,NX K-L+I AT(K)-CONJG(GXX(I,J) ) *ZT(K)+CONJG(GXY(I,J)) *ZT(K+NT) 50 CONTINUE 60 CONTINUE C*VDIR IGNORE RECRDEPS DO 80 J-1,NY L-NT+(J-1) *NX DO 70 I-1,NX K-L+I AT(K)-CONJG(GYX(I,J) ) *ZT (K-NT) +CONJG (GYY (I,J)) *ZT (K) 70 CONTINUE 80 CONTINUE END IF C C ------------------------- ----------------------- C C.. INVERSE FOURIER TRANSFORM C CALL SPECTR(A,AT, 1,W, ITAG) C O I C ----------- ------------------------ ------------— o --- C C...CALCULATIONS FOR RESISTIVE PLATES C IF(TPR.GT.-1.)THEN IF(ISIGN.EQ.-1) THEN C*VDIR IGNORE RECRDEPS DO 100 J=1,MY L=(J-1) *MX DO 90 I-1,MX K-L+I A (K) -RES (K) *Z (K) +A (K) 90 CONTINUE 100 CONTINUE C*VDIR IGNORE RECRDEPS DO 120 J-1,MY L-MT+ (J-1)*MX DO 110 I-1,MX K-L+I A(K)-RES (K-MT) *Z (K-MT) +A (K) 110 CONTINUE 120 CONTINUE ELSEIF(ISIGN.EQ.1) THEN C*VDIR IGNORE RECRDEPS DO 140 J-1,MY L-(J-1) *MX DO 130 I-1,MX K-L+I A (K) -RES (K) *Z (K) +A (K) 130 CONTINUE

Print file "'plt~cgf " 140 CONTINUE DO 160 J=lMY L-MT+(J-1) *MX DO 150 I-1,MX K-L+I A (K) =R.ES (K-MT) * Z (K-MT) +A (K) 150 CONTINUE 160 CONTINUE ENDIF END IF RETURN END Page 1 6 C C C * FOURIER TRANSFORM VIA FFT SUBROUTINE SPECTR (ZZT, ISIGN, WITAG) INTEGER N (2),fITAG (*) REAL*8 AUX1 (20000)rAUX2 (20000) COMPLEX Z (*),ZT(*),W(*) COMMON /DIM/MXMYNXNYNTMGMGIDXDYDS N (1) =NX N (2) -NY * C C....FORWARD FOURIER TRANSFORM C IF(ISIGN.EQ.-1)THEN DO 10 I - 1,2 C PRECONDITIONING IN SPATIAL DOMAIN DO 15 IIinlNT 15 CONTINUE Li - (I-1)*MG DO 20 L - 1,MG W(ITAG(L)) - Z(Ll+L) 20 CONTINUE C CALL SCFT2 (1,W, 1,NXW, 1,NXNXNY,l1,l.0,AUXl,20000,AUX2,20000) C CALL SCFT2(0,W,1,NXW,1,NXNXNY,1,1.0,AUXl,20000,AUX2,20000) CALL FFTNDM (W, N, 2 -1) Li = (I-1)*NT DO 40 L = 1,NT ZT (L1+L) = W (L) 40 CONTINUE 10 CONTINUE C C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C C... INVERSE FOURIER TRANSFORM C ELSE IF(ISIGN.EQ.1)THEN DO 50 I - 1,2 Li - (I-1)*NT DO 60 L - 1,NT W(L) - ZT(L1+L) 60 CONTINUE C CALL SCFT2(lW,1,NXW,1,NXNXNY,-1,1.0,AUXl,20000,AUX2,20000) C CALL SCFT2(0,W,1,NXWri,NXNXNY,-1,1.0,AUX1,20000,AUX2,?20000) CALL FFTNDM (W, Nr2,r1) Li - (I-1)*MG DO 70 L - 1,MG Z(L1+L) - W(ITAG(L)) 70 CONTINUE 50 CONTINUE END IF RETURN

Print file "pltcgf " Page 17 END C ======-=-5 =========SZ========5======================-===== C * ROUTINE FOR COMPUTATION OF FAR ZONE SCATTERING COEFFICIENT SUBROUTINE SCATER(SIG,JS,TETA,PHI,X,Y) REAL X(*),Y(*) COMPLEX JS(*),XJ,PHASE,SX,SY,NTETA,NPHI COMMON /DIM/MX,MY,NX,NY,NT,MG,MGI,DX,DY,DS DATA PI/3.141592653589793/,TPI/6.2831853071796DO/,XJ/(0.,1.)/ C ----- ---------------------------------------- STE = SIN(TETA*PI/180.) CTE = COS(TETA*PI/180.) SPH - SIN(PHI*PI/180.) CPH - COS(PHI*PI/180.) RX - TPI*STE*CPH RY - TPI*STE*SPH SX - (0.,0.) SY = (0.,0.) DO 10 I-1,MG PHASE - CEXP(XJ*(RX*X(I)+RY*Y(I))) SX = SX + JS(I) * PHASE SY - SY + JS(I+MG) * PHASE 10 CONTINUE SX - DS*SNSINC(DX*RX/2.) *SNSINC (DY*RY/2.) *SX SY - DS*SNSINC(DX*RX/2.)*SNSINC(DY*RY/2.)*SY NTETA - CTE*(CPH*SX+SPH*SY) NPHI =-SPH*SX+CPH*SY SIG = PI*(CABS(NTETA)**2+CABS(NPHI)**2) IF(SIG.LE..E-10)SIG-1.E-10 SIG - 10.*ALOG10(SIG) RETURN END C C C C C C C C C * SNSINC REAL FUNCTION SNSINC(A) IF(A.EQ.0)THEN SNSINC - 1. ELSE SNSINC - SIN(A) / A ENDIF RETURN END * DBSINC DOUBLE PRECISION FUNCTION DBSINC(A) REAL* 8 A IF(A.EQ. 0.DO)THEN DBSINC - 1.DO ELSE DBSINC - DSIN(A) / A END IF RETURN END * SQUARE ROOT OF A SIGNED REAL NUMBER * COMPLEX FUNCTION MYCSQR(A) REAL*8 A COMPLEX XJ DATA XJ/(O.,1.)/ IF (A. GE. ODO) MYCSQR-SNGL (DSQRT (A)) IF (A. LT.ODO)MYCSQR-XJ*SNGL (DSQRT(-A)). RETURN END

Prin C C C C C C C C C C C C C C C C t I* mm file "Pltcqf ', Page 18 -ile 'plicgf" Page 18 FFT ROUTINE FOR COMPUTATION OF N DIMENSIONAL FOURIER TRANSFORM TASKS: 1)BIT REVERSAL 2)TRIGONOMETRIC RECURRENCE 3)TRANSFORM ALL DONE IN THIS PROGRAM * * * = PARAMETER DESCRIPTION DATA..............A REAL ARRAY IN WHICH DATA ARE STORED AS IN A MULTIDIMENSIONAL COMPLEX FORTRAN ARRAY NDIM.................DIMENSION OF DATA AND THE FFT NN................... INTEGER ARRAY OF LENGTH NDIM ISIGN................DIRECTION OF THE TRANSFORM: -1 -FORWARD FFT 1 -INVERSE FFT TIMES THE PRODUCT OF LENGTHS OF ALL DIMENSIONS - - - s - a s - = - = - = = = = = = = = = = = = = = = = ==1~~~~ SUBROUTINE FFTNDM(DATA,NN,NDIM,ISIGN) REAL*8. WR,WI,WPR,WPI,WTEMP,THETA DIMENSION NN (NDIM), DATA ( *) NTOT-1 DO 10 IDIM-1,NDIM NTOT-NTOT*NN(IDIM) 10 CONTINUE NPREV 1 DO 80 IDIM-1,NDIM N-NN(IDIM) NREM-NTOT/ (N*NPREV) IP 1-2*NPREV IP2-IP1*N IP3-IP2*NREM I2REV-1 C C - ----------------------------- C C...BIT REVERSAL SECTION C DO 40 I2-1,IP2,IP1 IF (I2.LT. I2REV) THEN DO 30 I1-I2,I2+IP1-2,2 DO 20 I3-I1,IP3,IP2 I3REV-I2REV+I3-I2 TEMPR-DATA(I3) TEMPI-DATA (3+1) DATA (3) -DATA(I3REV) DATA (I3+1) -DATA (I3REV+1) DATA (I3REV) -TEMPR DATA (I3REV+1) -TEMP I CONTINUE CONTINUE END IF IBIT-IP2/2 IF((IBIT.GE.IP1).AND.(I2REV.GT.IBIT))THEN I2REV-I2REV-IBIT IBIT-IBIT/2 GO TO 1 END IF I2REV-I2REV+IBIT CONTINUE 20 30 1 40 C c -------------------—. C C C...DANIELSON-LANCZOS FORMULA C 2 IFP1-IP1 IF(IFP1.LT. IP2) THEN IFP2-2*IFP1

Print file "pltcgff" Page 19 THETA=ISIGN*6.28318530717959D/ (IFP2/IP1) WPR=-2.DO*DSIN(0.5D0*THETA) **2 WPI-DSIN (THETA) WR-1.DO WI-0.DO DO 70 I3-1,IFP1,IP1 DO 60 I1-I3,I3+IP1-2,2 DO 50 I2I1, IP3, IFP2 K1-I2 K2-K1+IFP1 TEMPR-SNGL (WR) *DATA(K2) -SNGL (WI) *DATA(K2+1) TEMPI-SNGL (WR) *DATA(K2+1) +SNGL (WI) *DATA(K2) DATA (K2) -DATA (K1) -TEMPR DATA (K2+1) -DATA (K1+1) -TEMP I DATA (K1) -DATA(K1) +TEMPR DATA (K1+1) -DATA (K1+1) +TEMPI CONTINUE CONTINUE WTEMP=WR 50 60 C C -------- - -------------------------------------------- C C...TRIGONOMETRIC RECURRENCE C WR=WR*WPR-WI*WP I +WR WI=WI *WPR+WTEMP *WP I+WI 70 CONTINUE IFP1=IFP2 GO TO 2 END IF NPREV-N*NPREV 80 CONTINUE RETURN END

BIBLIOGRAPHY 63

64 BIBLIOGRAPHY [1] T. K. Sarkar, A. Ercument, and M. Rao, "Application of FFT and the conjugate gradient method for the solution of electromagnetic radiation from electrically large and small conducting bodies," IEEE Trans. Antenna Propagat., AP-34(5):635-640, May 1985. [2] N. N. Bojarski, k-space formulation of the electromagnetic scattering problem, Technical Report AFAL-TR-71-5, March 1971. [3] P. M. Vanden Berg, "Iterative computational techniques in scattering based upon the integrated square error criterion," IEEE Trans. Antenna Propagat., AP-32:1063-1071, October 1984. [4] T. J. Peters and J. L. Volakis, "The application of a conjugate gradient fft method to scattering from thin planar material plates," IEEE Trans. Antenna Propagat., AP-36(4):518-526, April 1988. [5] M. F. Catedra, J. G. Cuevas, and L. Nuno, "A scheme to analyze conducting plates of resonant size using the conjugate gradient method and the fast fourier transform," IEEE Trans. Antenna Propagat., AP-36:1744-1752, December 1988. [6] C. Y. Shen, K. J. Glover, M. I. Sander, and A. D. Varvatsis, "The discrete fourier transform method of solving differential-integral equations in scattering theory," IEEE Trans. Antenna Propagat., AP-37:1032-1041, August 1989. [7] K. Barkeshli and J. L. Volakis, "Improving the convergence rate of the conjugate gradient fft method using subdomain basis functions," IEEE Trans. Antenna Propagat., AP-37(4), April 1989. [8] Y. Rahmat-Samii and R. Mittra, "Electromagnetic coupling through small apertures in a conducting screen," IEEE Trans. Antenna Propagat., AP-25(2), March 1977. [9] J. S. Hey and T. B. A. Senior, "Electromagnetic scattering by thin conducting plates at glancing incidence," Proceedings of the Physical Society, 72:981, 1958. [10] J. S. Hey, G. S. Stewart, J. T. Pinson and P. E. V. Prince, "The scattering of electromagnetic waves by conducting spheres and disks," Proceedings of the Physical Society, 69:1038, 1956.

65 [11] K. Barkeshli and J. L. Volakis, "A vector concurrent application of a conjugate gradient FFT algorithm to electromagnetic radiation and scattering problems," IEEE Transactions on Magnetics, 25(4), July 1989.