RL- 780 SCATTERING BY DISTRIBUTIONS OF SMALL THIN PARTICLES by David Allen Ksienski A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1984 Doctoral Committee: Professor Professor Professor Professor Professor Professor Thomas B.A. Senior, Chairman Chiao-Min Chu Ward D. Getty Wilfred Kaplan Chen-To Tai Herschel Weil

ABSTRACT SCATTERING BY DISTRIBUTIONS OF SMALL THIN PARTICLES by David Allen Ksienski Chairman: Thomas B.A. Senior The scattering of electromagnetic radiation by distributions of particles occurs in a variety of circumstances. At radio wave frequencies the operation of radar units is affected by rain and ice crystals suspended within clouds, while at higher, optical frequencies the amount of solar radiation reaching the earth's surface can be affected by pollutants in the upper atmosphere. This study is restricted to particles which are small compared to the wavelength of the illuminating electromagnetic radiation. The particles are assumed to be composed of a homogeneous lossy dielectric, with conductors characterized by large permittivities and lossy materials described by permittivities with relatively large imaginary components. In the investigation of scattering by small particles, it is often convenient and useful to solve for the scattered field in terms of a low frequency expansion (a Taylor series in powers of the maximum dimension of the particle over the wavelength of the incident electromagnetic field). Unfortunately, the usual techniques for computing the low frequency expansion fail if the particle is collapsed to a plate

with vanishing thickness. The difficulty arises from an unanswered problem in classical physics, the construction of a vector potential. A solution to this problem is obtained and presented in the context of low frequency scattering. In many cases of low frequency scattering, only the first term of the expansion is necessary to adequately characterize the scattered field. This is obtained by solving a static field scattering problem. Unfortunately, for particles which are very thin (but of finite thickness), existing numerical codes become highly unstable. An algorithm is developed expressly for the thin plate scattering problem which appears accurate over a wide class of thin plates, permitting arbitrarily shaped plates with complex permittivities. The solution is obtained using a finite element method with linear basis functions over triangular elements. The results of the program include the calculation of the dipole moments associated with the plates, and the manner in which these dipole moments affect the electrical properties of the entire distribution is discussed.

ACKNOWLEDGEMENTS The author wishes to express his gratitude for the guidance and expertise provided by the chairman, Professor Thomas B.A. Senior. The author also wishes to acknowledge the technical and emotional support of his colleagues at The University of Michigan. Special thanks are due Ms. Wanita Rasey and Ms. Jean Ringe for their skillful and expeditious typing of the manuscript. - i ii -

TABLE OF CONTENTS Page DEDICATION 1 ACKNOWLEDGMENT iii LIST OF FIGURES V LIST OF TABLES vi CHAPTER I. INTRODUCTION 1, II. LOW FREQUENCY SCATTERING 7 2.1 The Low Frequency Expansion 7 2.2 Scattering from a Plate 10 III. STATIC SCATTERING FROM A DIELECTRIC PLATE 18 3.1 The Dipole Moment 18 3.2 Static Scattering from a Thin Dielectric Plate 18 3.3 An Overview of Numerical Methods 24 3.4 Functional Description of the Program 27 IV. THE DIPOLE MOMENT: NUMERICAL RESULTS 48 4.1 Accuracy of the Results 48 4.2 A Linear Predictor Model 53 4.3 Shape Effects 69 V. SCATTERING BY DISTRIBUTIONS OF PARTICLES 73 5.1 The Clausius-Mosotti-Lorentz-Lorenz Equation 73 5.2 Cubical Distributions of Particles 80 5.3 A Cubical Exclusion Volume 82 5.4 Restriction of the Clausius-Mosotti Equation to Distributions of Plates 87 VI. CONCLUSIONS 91 APPENDIX 94 BIBLIOGRAPHY 172 -1iv

LIST OF FIGURES Figure Page 3.1 Diagram of a Plate. 20 3.2(a) Definition of Elements for a Circular Plate. 28 3.2(b) Input to Computer Program to Define a Circular Plate. 29 3.3(a) Definition of Elements for a Square Plate. 30 3.3(b) Input to Computer Program to Define a Square Plate. 31 3.4(a) Definition of Elements for a Rectangular Plate. 32 3.4(b) Input to Computer Program to Define a Rectangular Plate. 33 3.5 Dipole Moment and Calculated Potentials for the Square Plate Defined in Fig. 3.3. Under x Excitation, T = 2, k/t = 10. 42 3.6 Dipole Moment and Calculated Potentials for the Square Plate Defined in Fig. 3.3. Under x Excitation, T = 101, V/t = 10. 43 3.7 Perspective Plot of Potential on Top Surface of Square Plate Shown Partially Hidden by 0 = 0 Reference Square: x Excitation, z/t = 10, T = 2. 44 3.8 Perspective Plot of Potential on Top Surface of Square Plate Imposed on 0 = 0 Reference Square: z Excitation, Vt = 10, T = 2. 45 3.9 Perspective Plot of Potential on Top Surface of Square Plate Shown Partially Hidden by q = 0 Reference Square: x Excitation, k/t = 10, T = 101. 46 3.10 Perspective Plot of Potential on Top Surface of Square Plate Imposed on ) = 0 Reference Square: z Excitation, Z/t = 10, T = 101. 47 4.1 Diagram of a Plate Used by Linear Predictor. 55 5.1 Parallel Plate Capacitor with Spherical Exclusion Volume. 75

LIST OF TABLES Table 4.1 4.2 4.3 P /V 11 P /V 11 P /V 33 for for for Plates Plates P1 ates with Circular Cross-Section. with Square Cross-Section. with Square Cross-Section Page 50 51 52 4.4 P /V for Several Convex 4.5 P /V for Several Convex 3 3 4.6 P /V for Several Convex 1:1 4.7 P /V for Several Convex 2 2 4.8 P /V for Several Convex 3 3 4.9 P /V for Several Convex and complex permittivity 4.10 P /V for Several Convex 22 and complex permittivity 4.11 P33/V for Several Convex and complex permittivity 4.12 Eigenvalues Extracted by Plates with Length/Width = 1 Plates Plates Plates Plates Plates with with with with with Length/Width Length/Width Length/Width Length/Width Length/Width 1 2 2 2 2 59 60 61 62 63 64 65 Plates with Length/Width = 2 Plates with Length/Width = 2 Linear Predictor 66 68 -vi -

CHAPTER I. INTRODUCTION The interaction of electromagnetic radiation with distributions of particles is significant in many areas. For example, the manner in which radio waves scatter from water droplets produced by clouds and the ice crystals which may be suspended within these same clouds is important for the operation of radar units as well as the reception capabilities of radios and television sets. The interaction of higher frequency electromagnetic radiation with particle distributions is also important, and the scattering of visible light from contaminants produced by factories is a phenomenon which is often all too visible. For the purpose of analysis, the individual particles may be considered as composed of a lossy dielectric. Conductors are then characterized by large permittivities while lossy materials are characterized by permittivities with relatively large imaginary components. One class of interaction or scattering of electromagnetic waves from particles is low frequency scattering. Despite the name, a particular band of frequencies is not implied, but rather the ratio of the size of the particle to the wavelength of the electromagnetic radiation is restricted to be small. Although an exact range of ratios is not dictated, the degree to which the ratio is small determines the extent to which the low frequency approximations are valid. When low frequency scattering is analyzed, the scattered electromagnetic radiation from a particle may be modeled by expressing -l

-2 - the scattered field as a Taylor series expansion in powers of the maximum dimension of the particle over the wavelength of the incident electromagnetic field (Stevenson, 1954). In many cases, the ratio is sufficiently small to permit characterization of the scattered field by the first term in the series which is a zeroth order term and corresponds to solving a static field problem (infinite wavelength). The solution of the static field problem is invariably much easier than the solution of the general dynamic problem, and when the size of the particle is sufficiently small compared to the wavelength, it permits a compact description of the scattered field. Particles which are thin (one dimension smaller than the others) are important both as approximations which permit various simplifications in the analysis and as a region where existing numerical codes exhibit instabilities. Thin particles are also a common constituent of aerosols, with the relatively large surface area helping the particle remain suspended in the atmosphere. The numerical solution of the static field problem for bodies with axial symmetry has been presented in Senior and Willis (1982), and involves the evaluation of an integral over the surface of the body. Unfortunately, as the thickness of the body is decreased, numerical inaccuracies increase (Willis, 1982). The reformulation of this integral into a form which is extremely stable and accurate for thin plates is one of the contributions of this investigation. Thin plates have been studied by Harrington and Mautz (1975) using the electric field integral equation and by Inspektorov (1982) using the magnetic field integral equation. Harrington and Mautz (1975) approximate the plates using a resistive sheet and thus ignore

-3 - any effects resulting from a normal polarization of the sheet. Although Inspektorov (1982) considers finite thickness plates, his formulation becomes increasingly unstable as the plate thickness decreases. Both formulations consider the dynamic scattering problem and require the solution of a pair of coupled integral equations. For the problem of scattering by electrically small particles (i.e., the Rayleigh region), a simpler approach is possible. As shown by Keller et al (1972), the scattered field may be characterized by a polarization tensor. Further, if the object contains at least one axis of symmetry, at most three of the tensor elements are independent, corresponding to polarization along three perpendicular axes. The calculation of the dipole moments necessitates the previous calculation of the potential which may be obtained from a single scalar integral equation. This generally requires the solution of a potential problem (Senior, 1982), however this may be circumvented by using zero-degree harmonics (Senior and Ksienski, 1984). Unfortunately, when the scattering object is collapsed to a plate of infinitesimal thickness both of the techniques fail. Further, for thin plates (small but finite thickness), the solution of the potential problem develops numerical difficulties (Willis, 1982). Distributions of particles are important for examining composite materials (e.g., Polder and Van Santen, 1946; and Bergman, 1978), as well as clouds of particles, and regular arrays of elements designed to provide an artificial dielectric for the purpose of microwave lenses (e.g., Kock, 1948). For particles and wavelengths such that the particle interaction is far field and the particle size and inter-particle distance are small compared to a wavelength, the

-4 - distribution may be analyzed using the Clausius-Mosotti-LorentzLorenz formulation (e.g., von Hippel, 1954). If the particles are spherical, the Clausius-Mosotti-Lorentz-Lorenz formulation may be reduced to a formulation due to Maxwell Garnett (1904), which has been experimentally verified for sparse distributions of particles (for which the formulation is rigorously correct) as well as for dense distributions of particles (Bohren and Battan, 1980). This investigation is concerned with low frequency scattering from distributions of particles. The particles, as well as the inter-particle distance, are assumed small relative to the wavelength of the exciting field. Additionally, one dimension of the particle is assumed smaller than the other two dimensions, and in Chapter II the particles considered have zero thickness and are nothing more than a boundary condition on an open surface. Chapter II discusses the problem of low frequency scattering from a single thin particle. Although low frequency scattering has been previously discussed for a solid body, the problem of scattering from an open surface defeats the standard techniques (eg., Stevenson, 1954; Kleinman, 1965; and Senior, 1982). The solution of the problem involves the construction of a vector potential F, given v x F = f which is a problem from classical physics. Solving problems which involve scattering from particles with zero thickness is not merely an academic exercise. Particles which are very thin may be modeled by open surfaces, and this generally permits some simplification in the analysis. For example, scattering from resistive sheets was considered by Harrington and Mautz (1975) and shown to be an effective model for thin (finite thickness) dielectric shells.

-5 - Chapters III and IV consider the problem of static scattering from a thin dielectric plate. The solution of this problem is of course necessary in the low frequency expansion. In contrast to Chapter II, the problem considered in Chapters III and IV is primarily numerical. Previous studies in static scattering from dielectric particles have encountered numerical difficulties when the objects became thin (e.g., Senior, 1975; Herrick, 1976; and Willis, 1982). This is believed to result from the singularity associated with the surface integral formulation. In Chapter III the surface integral formulation of the problem is recast using a volume integral which has a less singular kernel. The integro-differential equation is solved for the potential which is constrained to vary linearly along its smallest dimension. The finite element method (e.g., Zienkowicz, 1982) is employed using triangular elements with linear basis functions, which guarantee Co continuity. The result of the approach is a highly efficient, very stable matrix problem with all of the matrix elements evaluated accurately using an analytic evaluation of the volume integral. The results are consistent with expectations and the program is able to operate over a wide range of thicknesses. The dipole moments for several shapes are presented and these are compared to the results of a simple linear predictor model which is derived in Chapter IV. Chapter V is concerned with sparse distributions of particles. These distributions of particles are analyzed in terms of the density and dipole moments of the constituent particles to arrive at an artificial dielectric description of the distribution. The primary

-6 - focus of this chapter is the rigorous derivation of the ClausiusMosotti-Lorentz-Lorenz equation (e.g., von Hippel, 1954), which describes the effective permittivity of the artificial dielectric. The chapter concludes with a comment on distributions of plates, and the accuracy of the method is verified for a dense cubical distribution of thin plates. Finally, it must be noted that this dissertation reflects the author's preference of electric quantities over magnetic quantities. Thus, the use of an electric exciting field, electric potentials, dielectric constants, and electric dipole moments. In fact, exactly analogous magnetic problems can be solved by merely reading magnetic for electric, permeability for permittivity, etc., which in Chapter IV would then result in the calculation of magnetic dipole moments.

CHAPTER II. LOW FREQUENCY SCATTERING 2.1 The Low Frequency Expansion The low frequency expansion involves an expansion of the scattered field in powers of the frequency. If the dimensions of the particle are small compared to the wavelength of the illuminating electromagnetic radiation then only a few terms of the expansion may be necessary to accurately characterize the scattered field. The numerical problems encountered with the low frequency expansion are generally simpler than those associated with the general dynamic problem. An additional benefit is that by expanding the solution explicitly in powers of frequency the scattered field becomes known for a range of frequencies, in contrast to the dynamic case where the problem must be resolved for each frequency desired. The method of the low frequency expansion has been developed in Stevenson (1953), Kleinman (1965), and Senior (1982). However, if the body is collapsed to a surface with zero thickness (perhaps as a model of a thin metal sheet), mathematical difficulties are encountered. There are several reasons for using the zero thickness model, not the least of which is the fact that the mathematical problems encountered are interesting. For solid bodies which are thin, numerical difficulties often arise in attempting to evaluate the associated integrals over the opposing surfaces which are close -7 -

-8 - together. Finally, by reducing the problem to a zero thickness plate it is hoped that simplifications would result in the analysis. In performing a low frequency expansion, the first step is to expand all quantities in powers of frequency. Specifically, the electric and magnetic fields, E and H, are expanded in powers of (ik), where k is related to the circular frequency w by k = ( 1/p) U and c are the permeability and permittivity of free space and the time convention used is e. The field quantities E and H may be expressed in terms of sources. Charge p and current J are sources which, since they arise from the scattering problem, are assumed to be constrained to the surface of the plate. These quantities are also expanded in powers of frequency, 00 00 E Em(ik)m,= 2 H (ik)m m=o m=o 00 00 PM(ik)m m m P = 2 (ik)", J = *(ik) m=o m=o When the expansions are inserted into Maxwell's equations V x H = J + at v x E

-9 - V. J = where B = AH and D = eE, the following relation is obtained inter alia: Z v x = -E (2.1 0 1 0 The constant Z0 is the free space impedance, H is the first order magnetic field, and Eo is the zeroth order, or static, electric field. Equation (2.1) also constitutes a problem which must be solved in the low frequency expansion. The problem is to determine R given E0 Since Eo is assumed known, it may be defined in terms of a scalar potential, q, which is also assumed known. The precise context in which (2.1) arises varies with the different treatments, and the problem is not generally stated in terms of a first order magnetic field and a zeroth order electric field. However, the formulation in (2.1) does give the problem some physical significance. The body is assumed to have zero net charge, which is mathematically stated by requiring Jn * v ds = 0 B The standard solution (Stevenson, 1954) for a solid body is then given by ZoR = -v x n(H - )G ds' B where G is the free space static Green's function and ~ is the solution to the interior Neumann problem

-1 0 - an an The boundary value of ap/an is known since E0 and hence vf is known everywhere. The solution of the interior Neumann problem does involve the solution of an integral equation, which from a numerical perspective is time consuming. This problem may be circumvented through the use of zero-degree harmonics, as shown in Senior and Ksienski (1984). Unfortunately, both the standard solution via the Neumann problem and the solution via zero-degree harmonics break down if the body is collapsed to a zero thickness plate. The solution via zero-degree harmonics breaks down if the volume of the body vanishes because the evaluation of the kernel associated with the integral becomes ambiguous. The method of solving the problem by first solving the Neumann problem becomes impossible since the Neumann problem is a three dimensional problem which must be solved in the interior of the body which is only a two dimensional surface of discontinuity. Further, the symmetry of the problem would appear to indicate that both p and D are continuous across the surface of the plate which forces H to be identically zero. As the solution of (2.1) -1 for H is necessary for the continued development of a low frequency 1 expansion, it is apparent that the expansion cannot be obtained for a zero thickness plate. A solution which is valid for the case of a zero thickness plate is the subject of the remainder of this chapter. 2.2 Scattering from a Plate In obtaining an H which satisfies Eq. (2.1), an additional constraint which is applied is that H must be physically reasonable.

-11 - If H is expressed in terms of a current distribution, this requirement 1 may be stated precisely by constraining the current to the scattering object and, in the case of the plate, by requiring that the current does not flow off the edge of the plate (i.e., the normal component of J on the edge of the plate must be zero). Obtaining a solution for H in terms of the current may be facilitated by solving I V *3 = cp (2.2) 11 S 0.2, for J, where 11 l+ Po c 0 " o z and c is the speed of light. Since the current J and charge p0 are both surface distributions, the divergence of 3 is specified with the surface differential operator v S The plate is assumed to lie in the x-y plane, and the notation ap/az|j denotes the discontinuity of aj/9z across the surface of the plate. The solution of (2.2) for J1 produces a solution of (2.1) if H1 is defined in terms of the current distribution as H = v x f J G ds' i ~J B This may be shown by substituting this definition of H into (2.1) and taking the field point away from the plate, taking the field point away from the plate,

-12 - E = -Z v x = -Z x x G ds') B = -Z vv G ds' B Bringing the differential operator inside the integral and converting the operator to the primed coordinate system yields o o Z 1 = Z f v'G ds' B = -Z v Gvs J ds' + Zov (3 G) T d 0 1 B CB where T is the outward normal to the plate on the edge of the plate and lies in the plane of the plate and CB denotes the contour surrounding the plate. The contribution of the second integral is zero if r * J is constrained to be zero as discussed above. 1 Finally, using equation (2.2) reduces the above equation to Eo v GpO ds' B which is the definition of Eo in terms of the charge distribution. It should be noted that both the solution (2.1) for H1 and (2.2) for J1 allow considerable freedom. Since only the curl of H1 is specified, H is only unique to within the gradient of a scalar. Similarly, in (2.2), only the divergence of J1 is specified and thus J is known only to within the curl of a vector. As only a particular 1

-13 - solution is required for (2.1) or (2.2), this ambiguity can be exploited to yield a solution which is easy to implement. Three solutions are proposed for (2.2), the first being the rather natural restriction of representing J as the gradient of a scalar. If J is represented as 1 J = v then from (2.2) V2 = cp, which is a two dimensional Poisson problem for b. Since the forcing function is cpo and p0 tends to infinity near the edges of the plate, the Poisson problem will in general be difficult to solve. However, if the geometry of the plate is simple, for example a circular plate, then the charge distribution resulting from a uniform incident electrostatic field is known analytically. The Poisson problem may then be solved analytically, and the solution for a circular plate is given in Senior and Ksienski (1984). Unfortunately, for a general plate geometry the solution of a Poisson problem with an unbounded forcing function is not a problem which is well suited for numerical methods, and an alternate solution is needed. -inc If J is now restricted to J = xJ assuming Eo = -vx then aJ X = Cp 2X0

-14 - and J = c/P dx By choosing J to lie in the same direction as the incident electric field, a solution is obtained by simply integrating the charge distribution. The integral specified above proceeds along the surface of the plate in the x direction. If the plate is symmetric and convex, then the normal component of 3 will vanish along the edge. The requirement that the plate be symmetric and convex is a loose description of a geometrical constraint which can be stated precisely. Specifically, the plate must have two non-collinear axes of symmetry such that the intersection of the plate with any line parallel to either axis is simply connected. The satisfaction of this requirement will then permit the determination of the field scattered by the plate when it is illuminated by an arbitrary uniform electrostatic field. However, if the plate is not symmetric and convex (in the sense defined above), it is impossible to force the normal component of the current to vanish along the edge, even with an arbitrary function of y which may be added to the integral. The third solution is intended as a correction to the preceding solution when the plate is not both symmetric and complex. The current is broken into two components, d = J(1) + j(2)

-15 - with J(1) chosen as above, and J(2) represented as the gradient of a scalar. Specifically, J(1) J(2) V2' 5 = X x = ViS = 0 with * J = O. The validity of J(1) is thus preserved from the 1 -(2) preceding solution, while J(2) as a homogeneous solution to (2.2) is added to satisfy the boundary condition. The solution for J(2) is obtained by solving a two dimensional Neumann problem, with boundary conditions which are finite since they arise from the normal component of J(1)along the edge of the plate. The normal component of J(1) along the edge of the plate is finite since it is defined as J i) = xc dx, -1/2 and the singularities in p0 are of the order x. Finally, the two dimensional Neumann problem does have a solution since J vs T' dk' = - CB j(1) T' dz CB (cont.)

= - fvs j ( ) ds B = - f cP0 ds B by the zero net charge condition. This last solution is valid for any flat plate, and thus the low frequency expansion is now repaired. The solution which has been obtained was intended to be complementary to a solution for solid perfectly conducting bodies obtained by Senior (1982). This goal has been achieved (Senior, 1983), and the solution has also proven useful for the problem of scattering from a resistive plate (Senior and Naor, 1984), as well as providing a solution for a problem from classical physics, the construction of a vector potential (Senior and Ksienski, 1984). The above formulation is necessary in analyzing the scattering by a thin plate via the low frequency expansion, and first and higher order terms in the expansion may now be obtained through the methods presented in this chapter. The low frequency expansion is valid when the dimensions of the particle are small compared to the wavelength of the illuminating radiation. To the extent that this requirement is satisfied, only a few terms in the low frequency expansion may be necessary to accurately characterize the scattered field. For a sufficiently small particle only the zeroth order term is necessary

-17 -to accurately characterize the scattered field, and this is considered in the remainder of the dissertation.

CHAPTER III. STATIC SCATTERING FROM A DIELECTRIC PLATE 3.1 The Dipole Moment In Chapter II, the problem considered was that of finding a first order magnetic field given a zeroth order electric field. The zeroth order fields are important not only as the first term in a low frequency expansion upon which the higher order terms depend, but also in their own right. The zeroth order fields determine the electric dipole moment p (and magnetic dipole moment m), which often is all the information that is needed. For example, in sparse distributions of particles, discussed in Chapter V, the distribution will be characterized by the density of the distribution and the dipole moments of the individual particles. Although the zeroth order term and hence the dipole moment are most easily visualized in a static field, the concept of a low frequency analysis permits the use of a dielectric with complex permittivity which has physical significance in a dynamic field. Thus the dipole moment calculated for complex permittivities can provide much information about the far field scattered from a lossy dielectric particle, and for the analysis contained in Chapter V it is all the information that is required. 3.2 Static Scattering from a Thin Dielectric Plate In the analysis of scattering from a thin flat dielectric plate the usual surface integral formulation (Senior, 1976) t = 2 1 - T 2 t G, dx'(31) 1 + T xj +1 + Tr dx -18 -

-19 - which is valid for any solid dielectric body with permittivity (permeability) - becomes highly unstable. This difficulty is believed to arise from the highly singular kernel and the problem of correctly evaluating contributions from the integration on the edge of the plate. In a previous numerical investigation (Willis, 1982), the instability was found to be worst for very thin plates and for large permittivities. As the current investigation is directed at computing the scattered fields specifically, though not exclusively, for the case of plate thickness approaching zero, coupled with permittivities approaching infinity, it was felt that an alternate formulation was needed. An integral equation was sought which involves integration only over a single flat surface. Although this was not obtained, an integrodifferential equation was found which involves integration and differentiation over a single flat surface. The associated kernel has a very mild singularity, and this new integro-differential equation resulted in greater numerical stability than Eq. 3.1. Without loss in generality, the plate is assumed to lie centered in the z = 0 plane. The plate is of thickness t, and the surface of the plate comprises the edge and the two parallel walls, as shown in Fig. 3.1. The general scattering problem may be solved by decomposing the incident electric field into x,y, and z components and then superposing the inc associated scattered fields. If E = x, then (3.1) becomes t 2 1 + - T 2 it aG ds' + ft aG d e w where e represents the edge and w represents the two parallel walls.

-20 - I length I "'thickness Fig. 3.1: Diagram of a Plate.

-21 - t/2 t 2 1 - X T t -- 1 + T 1 + - c -t/2 t/2 + - a s -t/2 (t ) dz' ds' where c is the intersection of the the intersection of the plate with of thickness t, and extends from z t/2 t = 2 +1 -T 2 p l +T 1 + J: -t/2 s edge with the the plane z = = -t/2 to z = vs ~ ( stvG. ) S. plane z = z' and s is z. The plate is t/2. + - t -( | ds' dz' 2 1 - T - + T l+T t/2 -t/2 Vs.V'G + t V2G 's s5 + azI z +,tz - ds' dz' t/2 = -1 T-2r F 1 +T ' +T t/2 - t/2 S -pt 6(R-Rl + V5t.5 at aG _ ds' dz' DZ' ~Z' where R is the field point and R' is the source point. Noting that R is on the boundary of the volume of integration,

-22 - t/2, J + t t 2 -T t - tI G d (1 - - + 2 (1 - V)2 J vsV SG + z X ds' dz' -t/2 S -t/ 2 s' 1 t = -x + (1 - T) t/2 t J |S. vsG + X ds' dz' -t/2(3.2) (3.2) Since no approximations have been made yet, Eqs. (3.1) and (3.2) have identical solutions. If ~t is now expanded in powers of (z/t), it is apparent from the symmetry of the problem that 00 = >2 t4 (z/t)2i i=o (3.3) i.e., odd powers of (z/t) are ot will be kept in which case t -x o -x + unneeded. To simplify Eq. (3.2) only the first equation to be solved is t/2 (1 - T) t J fv't v'G ds' dz'. (3.4) - t/2 S An alternate derivation of Eq. (3.4) can be obtained by using the divergence theorem. Starting with Eq. (3.1)

-23 - t 2 1 - T 2 S' (ty'G) dv' +T 1+T v =2 x + tI 2G - + -x + T r2 Vv'.v'G + T'2G dv' t/2 1+T x - 1 + 2 s fv v 'G ds dz' -t/2 S which may be reduced to Eq. (3.2). For Ein = z, Eq. (3.1) becomes 2t 2 +2 t Gds' + ft =G ds' 1+ 1 + T 2 --- ds' e w and following the same procedure as before we obtain t/2 t t = -z + (1 - T) t vt. vG + t aG ds' dz' (3.5) -- S S DZ aZI' -/2 s inc t without making any approximations. For Ei = z, t may be expanded in powers of (z/t), co t t2i+ 1 2i f = 2i +1 (/t) i=0 and this time the even powers of (z/t) are unneeded. Keeping the first term in the expansion, and neglecting the contribution of V,' V'G since it results in only quadrupole and higher order trms terms

-24 - t/2 z + - ') J ds' dz' -t/2 S t/2 -z + (1 - T) ft G ds' s -t/2 Restricting R to lie on the top surface of the plate (z = t/2) = -t+ G(z = t/21z' = t/2) 21 t - G(z = t/21z' =-t/2) ds' (3.6) The static scattering problem has now been formulated for a thin dielectric plate with the incident electric field either tangential (Eq. 3.4) or normal (Eq. 3.6) to the plane of the plate. The numerical solution of these integral equations occupies the remainder of this chapter. 3.3 An Overview of Numerical Methods Before discussing the numerical implementation of the integrodifferential equation developed in Section 3.2, a comment about numerical methods as encountered in electromagnetic problems is perhaps appropriate. First, the variable for which the solution is sought often constrains the accuracy which can be obtained. If solutions are expressed in terms of the variable p, the charge density, then problems may be encountered near edges, where p tends to infinity.

-25 - Many solutions are posited in terms of J, the current density, which as the surface integral of p is somewhat better behaved. The accuracy that a particular numerical method can achieve is dictated by the degree to which any approximations made in the discretization are justified. For example, in finite difference codes a derivative of a function is often approximated by taking the difference of the value of the function at two nearby points and normalizing by the distance between the two points. This is a reasonable approximation for a derivative if the function is nearly linear, however if the function is not well behaved, such as a Green's function near the point of singularity, the approximation degenerates. In finite element analysis, the domain of the problem is divided into elements, upon which basis functions are then individually imposed. The two most common basis functions are pulse functions, which are constant over individual elements and generate discontinuities between elements, and linear functions, which permit linear variation over the individual elements and can guarantee CO continuity over the entire domain. Thus, an integral equation which can be expressed in terms of a variable which is approximately linearly varying should be amenable to solutions using linear basis functions. The use of pulse functions is not an indication of belief that the answer contains step functions, but rather a concession to the numerical difficulties of using more complicated basis functions. Higher order basis functions exist (Zienkowicz, 1982), however only C continuity may be enforced in general, and numerical difficulties seem to have precluded their popularity in electromagnetics. On the other hand, instead of partitioning the object into separate domains, another technique is to

-26 - define basis functions over the entire object. Although this technique is useful (Harrington, 1982), it is predicated upon at least some a priori knowledge about the nature of the solution. After the problem has been discretized, there often remain problems of numerical integration, if an integral is involved or numerical differentiation if a derivative is needed. The numerical differentiation is usually handled as part of the discretization and any attendant errors can be analyzed in terms of the discretization. The numerical integration often involves a substantial amount of additional computation time as the integrals over each element must be evaluated individually. The program which was developed over the course of this investigation analytically evaluates the integrals which are defined over the individual elements so that no approximations are required after discretization. The basis functions used are linear and are defined over triangular elements. The variable which is solved for is the potential, l, which is even more stable than the current, 3 since J is the derivative of p Additionally, in most cases + does indeed appear to be effectively linear, thus justifying the use of linear basis functions and permitting accurate solutions with a minimal number of subdomain divisions. Further, this leads to a simplified model of the program which can produce fairly accurate predictions of the dipole moment associated with a particular particle given gross parameters such as the height, length, width and permittivity. The range of applicability of the model is any convex plate, and the model and its results are discussed in Chapter IV,

-27 - 3.4 Functional Description of the Problem The first operation required in the numerical solution of the static scattering problem must be performed by the user. The plate must be divided into a set of triangular elements. This process can also be performed by automatic mesh generating programs, but this refinement was not felt necessary for the present investigation. The triangular elements are defined in terms of the vertices which delimit the triangles, and the vertices are defined by their coordinates. Figures 3.2(a), 3.3(a), and 3.4(a) show subdivisions of the circular, square, and rectangular plates into triangles with the triangles and points numbered. Figures 3.2(b), 3.3(b), and 3.4(b) are the associated inputs to the program which then define these subdivided plates. The lines beginning with the letter "p" indicate the point number, followed by the x and y coordinates. The lines beginning with the letter "t' indicate triangle number, followed by the numbers of the three points which delimit the triangle. The s3 line at the bottom indicates that the definition is only for the first quadrant and that the actual plate consists of this definition mirrored about both the line x = 0 and the line y = 0. Other options are sO, indicating no mirroring; sl, indicating mirroring about the line x = 0; and s2, indicating mirrorinc about y = 0. These options are used with Figs. 3.1 and 3.2 to generate a half circle and a triangle which are used in the next chapter. There are effectively only two restrictions about how the plate can be defined. Point numbers and triangle numbers must be given in increasing order starting with zero, but the actual numbering of the plates and triangles is arbitrary. The plate need not be simply connected, can be made up of any number of triangles and

-28 - 4 p3 3 p2 y 2 pi 0 PO p4 p8 0 1 2 3 p1l 4 X Fig. 3.2(a): Definition of Elements for a Circular Plate.

-29 - hl6 sided approximation to circle, using 12 point definition pO 0 0 pl 0 1.5 p2 0 3 p3 0 4 p4 1.5 0 p5 2 2 p6 1 3 p7 1.5307337 3.69551813 p8 3 0 p9 3 1 plO 2.828427125 2.828427125 pll 4 0 p12 3.69551813 1.5307337 tO 0 4 5 tl 0 1 5 t2 1 5 6 t3 1 2 6 t4 2 6 3 t5 3 6 7 t6 6 7 10 t7 6 5 10 t8 5 9 10 t9 4 5 9 tlO 4 8 9 tll 8 9 11 t12 9 11 12 t13 9 10 12 s3 Fig. 3.2(b): Input to Computer Program to Define a Circular Plate.

-30 - 4 p4 3 p3 y 2 p2 1 P ' 0 PO 0 p5 p9 p12 1 2 3 X p14 4 Fig. 3.3(a): Definition of Elements for a Square Plate.

-31 - h14 point definition of square using 2 planes of symmetry pO 0 0 pl 0 1 p2 0 2 p3 0 3 p4 0 4 p5 1 0 p6 1 1 p7 1 2 p8 1 3 p9 2 0 plO 2 1 pll 2 2 p12 3 0 p13 3 1 p14 4 0 tO 0 1 6 tl 1 2 7 t2 2 3 8 t3 3 4 8 t4 0 5 6 t5 1 6 7 t6 2 7 8 t7 5 6 10 t8 6 7 11 t9 7 8 11 tlO 5 9 10 tll 6 10 11 t12 9 10 13 t13 10 11 13 t14 9 12 13 t15 12 13 14 s3 Fig. 3.3(b): Input to Computer Program to Define a Square Plate.

-32 - 2 p3 p2 19 p18 p17 y 1 pi 0 pO p4 p8 p12 p16 0 1 2 3 X 4 Fig. 3.4(a): Definition of Elements for a Rectangular Plate.

-33 - h20 point definition of rectangle pO 0 0 pl 0 1 p2 0 1.66666666 p3 0 2 p4 1.66666666 0 p5 1.66666666 1 p6 1.66666666 1.66666666 p7 1.66666666 2 p8 3 0 p9 3 1 plO 3 1.66666666 pll 3 2 p12 3.66666666 0 p13 3.66666666 1 p14 3.66666666 1.66666666 p15 3.66666666 2 p16 4 0 p17 4 1 p18 4 1.66666666 p19 4 2 tO 0 4 5 tl 0 1 5 t2 1 5 6 t3 1 2 6 t4 2 6 7 t5 2 3 7 t6 4 8 9 t7 4 5 9 t8 5 9 10 t9 5 6 10 tlO 6 10 11 tll 6 7 11 t12 8 12 13 t13 8 9 13 t14 9 13 14 t15 9 10 14 t16 10 14 15 t17 10 11 15 t18 12 16 17 t19 12 13 17 t20 13 17 18 t21 13 14 18 t22 14 18 19 t23 14 15 19 s3 Fig. 3.4(b): Input to Computer Program to Define a Rectangular Plate.

-34 - plates, and the triangles can be of any shape desired. Parameters such as thickness, permittivity, and the direction of the incident field are specified after the elements of the plate are defined to permit respecification of these quantities while retaining the definition of the elements. The program generates all linking information such as which triangles have which points in common, which points are connected to other points, as well as which points delimit the perimeter of the plate, the last being needed in calculating the dipole moment. The static scattering problem can be solved assuming an incident field in either the x or y directions (tangential to the plane of the plate) or in the z direction (normal to the plane of the plate). The evaluation of Eq. (3.4) is facilitated by first breaking S, the domain of integration, into triangular subdomains. In these regions, pt is restricted to be linearly varying, so that for each element t 4o = CxX + Cyy + CO Then define a u-v coordinate system, with u = c, where c = c x + c y t/2 t _G / S G c ds' dz' qo = -x + (1 - T) f cl c| -dsI dz'. -t/2 S Since the boundaries of S. are line segments, the integral in the above equation can be evaluated as

-35 - t/2 -t/2 V2 dz 'f 1 -1/2 dv' {(av' + b - U)2 + (v' - v)2 + (Z' - z)2 t/2 - t/2 Ii dz' First evaluate V2_ I,=f {(l + a2)v'2 + 2(a(b - u ) - v)v' + (b - U) + v2 + (Z I - Z 9. 1/2dv Sv From Gradshteyn and Ryzhik (1980, Sec. 2.261) V f{A + 1 v'c7 By'I + CV 12} -1/2 dv' V 2 1n(2/CR + 2Cv'I + B) V 1 and A= (b - u )2 + v2 + (z'I - z)2 B = 2(a(b - u ) - v) (1 + a2 )

-36 - Ii = n ( 2 vr'VT~+a{ (l + + (b - U)2 + V2 + (Z' - Z)2} 1/2 + a2)V'2 + 2(a(b - u) - v)v' V '=V2 2(1 + a2)V' + 2(a(b - u) - v)), Take the field point on the top surface, z = t/2. To find I, let pj =ZI - t/2. -t 1 ln(2'l az(l + a2)V'12 + 2(a(b - u) - V)v I VI =V2 - V)) + (b - U )2 + V2 +T2 12 + 2(1 + a2)V' + 2(a(b - u) then, change parameters A = (1 + a)'/ B = (1 + a2)v'2 + 2(a(b - u) - v)v' + C = 2 (1 + a2) VI + 2 (a(b - u ) - v) (b _ U)2 + V2 = (av'+b-u)2 -t I ln(2A{B A 1= Then, integrating by parts I1= 2A ln(2A{B + T2 / V 0 2 + C) V -t 1 (cont. )

-37 - 0 _-n ~ 2A 1 + Tl2]7-1/ 2 2nI A 2A{B + -n2) 1/2wC 2-{ V 2 V 1 dpj To evaluate the second term, callI 21 (B + n12) -1/2 2A(B + T,2)1/2 + c it I 2 V 2 dvi -I 0 22 (B + y,2 -1/2 I 2-22A(B + T1) +C2 V 2 V 1 Let ~ = (B + T2 / p T = (~ 2 - B)" 21dri = E r- ~ ~ d C I t 2(~ ) J 2AC + C 1 V 2 V 1 Let y = 2AE~ + C, C = (y - C)/2A, d~ = 1/2A dy 2A t2+B+C i-ff 2 2A WC 2 2A - - C) 2 - B4A72 y V2 1 2A V1 dy

-38 - 2 2 AA2Z7+ C V 2A2 2 V1 dy Changing parameters again, z = 2A/ + C 1 a = C - 4BA b = -2C z I 2 = 2A/Ft-2 + C z2 Ic' T ~b~ 2z V 2 V 1 dy 2/V 13 z2 I z 1 V2 vra- + by + y y dy ' V I Then, from Gradshteyn and Ryzhik (1980, sec. 2.267.1) (with R = a+by+y.2). z V Zq 2 13 A + a z V 1 1 2 1 V V 2 2 V V 1 1 This may be further reduced by using Gradshteyn and Ryzhik (1980, Secs.- 2.266 and 2.261);

-39 - 2 V V Z I, Z 2 2+a- 1 n 2a + by + 2 t 2 V z 1 1 1 V z 2 2 + b {ln(2i+ 2 + 2y + b)} V z 1 1 V 0 + 1/2 2 I = In (2A{B + 2} + C) - I A 3 V1 -t For the electric field polarized parallel to the plane of the plate the problem is described by Eq. (3.6). This equation, with its integral over the surface of a triangle, can be evaluated with the aid of the techniques described in Wilton et al (1984). In the numerical solutions of Eqs. (3.4) and (3.6) the integrals occasionally degenerate into simpler forms depending on the shape of the triangular subdomain and its position relative to the field point. The program tests for these problems and uses alternate expressions as appropriate. For x or z excitation the contributions of the individual elements are determined analytically. The problem associated with the electric field polarized in the y direction is entirely analogous to the problem associated with the electric field polarized in the x direction, which is discussed above. The matrix problem is then generated by relating the elemental contributions to area coordinates (Zienkowicz, 1982) which are in turn

-40 - defined by the vertices of the triangles. Thus, the problem is formulated in terms of potentials at vertices. Further, the program will set to zero the potentials of points lying on the axes, if warranted by the symmetry of the plate and the direction of excitation. The formulation yields a very stable matrix formulation for most values of T (at least for Real T > 0), since the self cell terms of the matrix formulation of the integro-differential equation are usually the largest elements of the matrix, resulting in extremely low condition numbers for the matrix. After the potentials are obtained, the dipole moment normalized by the volume is obtained using the formulae for a symmetric dielectric body given in Senior (1976). A sample output for the square plate subjected to an incident electric field polarized in the x direction is shown in Fig. 3.5 with T = 2 and length-to-height ratio equal to ten. Figure 3.6 shows the output obtained from the same problem but with T changed to 101. In the solution of any complex numerical problem, it is important to examine the extent to which the initial approximations are justified. Since the basis functions are linear the solution should be approximately linear for the solution to be considered accurate. To determine this, perspective plots were generated showing the variation of the potential across the top surface of the plate. For the case of a square disk with length-to-height ratio of ten, the potential has been plotted for T = 2 with x excitation in Fig. 3.7, for r = 2 with z excitation in Fig. 3.8, for T = 101 with x excitation in Fig. 3.9, and for T = 101 with z excitation in Fig. 3.10. For the electric field in the x direction the computed potential does indeed appear approximately linear. The electric

-41 - field in the z direction produces some "buckling" near the edges, particularly for T = 101, which causes the elemental divisions to become visible. This might indicate that perhaps another formulation might be appropriate, although the solution is quite adequate for the current investigation. The problem could be overcome by shrinking the element size toward the edges, as was done for the rectangular shape, or simply using more elements. In Figs. 3.7 through 3.10, the middle structure shows the variation of the potential across the top surface of the square plate. The f = 0 reference square is not part of the solution and is included only to lend perspective. The structure is described by quadrilaterals, and the variation in the shape and size of these quadrilaterals illustrates the variation in the gradient of the potential. In Fig. 3.7 the potential is approximately linear and this should be contrasted to Figs. 3.8 through 3.10 where the potential exhibits some nonlinearities.

-42 - condition number is 1.504066e+00 14 point defintion of square using 2 planes of symmetry t = 0.800000, tau = 2.000000 + i 0.000000 potential has been computed assuming x excitation plate has mirror symmetry about x=0 and y=0 The dipole moment is 0.913873 + i 0.000000 point # 0 is located at x = 0.000000, y = 0.000000 and has potential 0.00000e+00 point 0 is associated with points and triangles (tri,Pl,P2) (0,1,6), (4,5,6), point # I is located at x = 0.000000, y = 1.000000 and has potential 0.000000e+00 point 1 is associated with points and triangles (tri,Pl,P2) (0,6,0), (1,2,7), (5,6,7), point # 2 is located at x = 0.000000, y = 2.000000 and has potential 0.000000e+00 point 2 is associated with points and triangles (tri,P1,P2) (1,7,1), (2,3,8), (6,7,8), point # 3 is located at x = 0.000000, y = 3.000000 and has potential 0.000000e+00 point 3 is associated with points and triangles (tri,,lP2) (2,8,2), (3,4,8), point # 4 is located at x = 0.000000, y = 4.000000 and has potential 0.000000e+00 point 4 is associated with points and triangles (tri,P1,P2) (3,8,3), point # 5 is located at x = 1.000000, y = 0.000000 and has potential 9.372243e-01 point 5 is associated with points and triangles (tri,P1,P2) (4,6,0), (7,6,10), ( 10,9,10), point # 6 is located at x = 1.000000, y = 1.000000 and has potential 9.305898e-01 point 6 is associated with points and triangles (tri,P1,P2) (0,0,1), (4,0,5), (5,7,1), (7,10,5), (8,7,11), (11,10,11), point # 7 is located at x = 1.000000, y = 2.000000 and has potential 9.083354e-01 point 7 is associated with points and triangles (tri,P1,P2) (1,1,2), (5,1,6), (6,8,2), (8,11,6), (9,8,11), point # 8 is located at x = 1.000000, y = 3.000000 and has potential 8.580674e-01 point 8 is associated with points and triangles (tri,Pl,P2) (2,2,3), (3,3,4), (6,2,7), (9,11,7), point # 9 is located at x = 2.000000, y = 0.000000 and has potential 1.872970e+00 point 9 is associated with,..ts and triangles (tri,P1,P2) (10,10,5), (12,10,13), (14,12,13), point # 10 is located at x = 2.000000, y = 1.000000 and has potential 1.856799e+00 point 10 is associated with points and triangles (tri,Pl,P2) (7,5,6), (10,5,9), (11,11,6), (12,13,9), (13,11,13), point # 11 is located at x = 2.000000, y = 2.000000 and has potential 1.803127e+00 point 11 is associated with points and triangles (tri,Pl,P2) (8,6,7), (9,7,8), (11,6,10), (13,13,10), point # 12 is located at x = 3.000000, y = 0.000000 and has potential 2.805956e+00 point 12 is associated with points and triangles (tri,P1,P2) (14,13,9), (15,13,14), point # 13 is located at x = 3.000000, y = 1.000000 and has potential 2.767921e+00 point 13 is associated with points and triangles (tri,P1,P2) (12,9,10), (13,10,11), (14,9,12), (15,14,12), point # 14 is located at x = 4.000000, y = 0.000000 and has potential 3.763731e+00 point 14 is associated with points and triangles (tri,PI,P2) 15, 12, 13), Fig. 3.5: Dipole Moment and Calculated Potentials for the Square Plate Defined in Fig. 3.3. Under x Excitation, T = 2, Z/t = 10.

-43 - condition number is 9.026534e+00 14 point defintion of square using 2 planes of symmetry t = 0.800000, tau = 101.000000 + i 0.000000 potential has been computed assuming x excitation plate has mirror symmetry about x=0 and y=0 The dipole moment is 10.727467 + i 0.000000 point # 0 is located at x = 0.000000, y = 0.000000 and has potential 0.000000e+00 point 0 is associated with points and triangles (tri,P1,P2) (0,1,6), (4,5,6), point # 1 is located at x = 0.000000, y = 1.000000 and has potential 0.000000e+00 point 1 is associated with points and triangles (tri,P1,P2) (0,6,0), (1,2,7), (5,6,7), point # 2 is located at x = 0.000000, y = 2.000000 and has potential 0.000000e+00 point 2 is associated with points and triangles (tri,Pl,P2) (1,7,1), (2,3,8), (6,7,8), point # 3 is located at x = 0.000000, y = 3.000000 and has potential 0.000000e+00 point 3 is associated with points and triangles (tri,Pl,P2) (2,8,2), (3,4,8), point # 4 is located at x = 0.000000, y = 4.000000 and has potential 0.000000e+00 point 4 is associated with points and triangles (tri,P1,P2) (3,8,3), point # 5 is located at x = 1.000000, y = 0.000000 and has potential 1.159421e-01 point 5 is associated with points and triangles (tri,P1,P2) (4,6,0), (7,6,10), (10,9,10), point # 6 is located at x = 1.000000, y = 1.000000 and has potential 1.106113e-01 point 6 is associated with points and triangles (tri,P1,P2) (0,0,1), (4,0,5), (5,7,1), (7,10,5), (8,7,11), (11,10,11), point # 7 is located at x = 1.000000, y = 2.000000 and has potential 9.526785e-02 point 7 is associated with points and triangles (tri,P1,P2) (1,1,2), (5,1,6), (6,8,2), (8,11,6), (9,8,11), point # 8 is located at x = 1.000000, y = 3.000000 and has potential 6.743818e-02 point 8 is associated with points and triangles (tri,P1,P2) (2,2,3), (3,3,4), (6,2,7), (9,11,7), point # 9 is located at x = 2.000000, y = 0.000000 and has potential 2.354477e-01 point 9 is associated with points and triangles (tri,P1,P2) (10,10,5), (12,10,13), (14,12,13), point # 10 is located at x = 2.000000, y = 1.00000 and has potential 2.238533e-01 point 10 is associated with points and triangles (tri,Pl,P2) (7,5,6), (10,5,9), (11,11,6), (12,13,9), (13,11,13), point # 11 is located at x = 2.000000, y = 2.000000 and has potential 1.906126e-01 point 11 is associated with points and triangles (tri,P1,P2) (8,6,7), (9,7,8), (11,6,10), (13,13,10), point # 12 is located at x = 3.000000, y = 0.000000 and has potential 3.639483e-01 point 12 is associated with points and triangles (tri,Pl,P2) (14,13,9), (15,13,14), point # 13 is located at x = 3.000000, y = 1.000000 and has potential 3.383455e-01 point 13 is associated with points and triangles (tri,P1,P2) (12,9,10), (13,10,11), (14,9,12), (15,14,12), point # 14 is located at x = 4.000000, y = 0.000000 and has potential 5.236022e-01 point 14 is associated with points and triangles (tri,Pl,P2) (15,12,13), Fig. 3.6: Dipole Moment and Calculated Potentials for the Square Plate Defined in Fig. 3.3. Under x Excitation, T = 101, z/t = 10.

-44 - y Fig. 3.7: Perspective Plot of Potential on Top Surface of Square Plate Shown Partially Hidden by - = 0 Reference Square: x Excitation, z/t = 10, T = 2.

-45 - yI x Fig. 3.8: Perspective Plot of Potential on Top Surface of Square Plate Imposed on d = 0 Reference Square: z Excitation, z/t = 10, T = 2.

-46 - y x Fig. 3.9: Perspective Plot of Potential on Top Surface of Square Plate Shown Partially Hidden by p = 0 Reference Square: x Excitation, ~/t = 10, T = 101.

-47 - X Fig. 3.10: Perspective Plot of Potential on Top Surface of Square Plate Imposed on ~ = 0 Reference Square: z Excitation, ~/t = 10, T = 101.

CHAPTER IV. THE DIPOLE MOMENT: NUMERICAL RESULTS 4.1 Accuracy of the Results The formulations developed in Chapter III, Eqs. (3.4) and (3.6), achieve the solution for scattering from a dielectric plate by assuming that the plate is thin. The internal field is represented as = f(x,y)(a + bz) (4.1) where a or b is zero depending on whether the excitation is parallel or normal to the plane of the plate. The function f(x,y) is approximated by a set of continuous and piecewise linear triangular elements. The extent to which the function f(x,y) is indeed approximately linear determines the number of elements needed to accurately characterize f(x,y), which in turn determines the matrix size and the amount of computation time needed to solve the problem. It was felt that the approximation represented by Eq. (4.1) would be valid for thickness-to-length ratios of 0.1 or smaller. This was the region for which the thin plate formulation was developed, since other solid body formulations develop numerical difficulties in this region (e.g., Willis, 1982). As the program was developed to work in a region where existing programs do not work, there is no reliable data to which the results may be compared. Verification of the program is thus limited to searching for -48 -

-49 - inconsistencies in results for limiting cases (none were found), and then comparing specific results with existing calculations (which are of limited accuracy) for approximate agreement. Beyond indicating approximate agreement, the reference data should not be considered as a basis for determining the accuracy of the present program, since in general the results of the present program are believed to be more accurate than the reference data. Of the data available (thickness to length ratios greater than 0.1), ratios of 0.1, 0.2, and 0.5 were selected for comparison with the program. The data are only believed accurate to within a few percent, and are included as reference merely to verify general trends in the dipole moment such as variations with thickness and permittivity. The reference data for the square cylinder was obtained from Herrick (1976) and the reference data for the circular cylinder was obtained from Senior (1975). The case of tangential excitation produced reasonable agreement with the reference data. For both the circular and square cylinders (Tables 4.1 and 4.2), the largest difference at the thickness to length ratio of 0.1 was 16 percent which occurred for the square at T = 106. The next largest error was only eight percent. For thickness to length ratios of 0.5 and 0.2 the largest error was 32 percent and occurred for a thickness to length ratio of 0.5 and T = 106. The data for normal excitation was also in reasonable agreement (Table 4.3). The fact that for the z excitation the dipole moment was overestimated relative to the reference data may be indicative of insufficiently fine elements. As shown in Fig. 3.10, the z problem does produce "buckling" near the edges which exhibits the elemental nature of the plate formulation and

-50 - Table 4.1 P /V for Plates with Circular 11 thickness/l ength 0.5 0.2 -1.418 -1.234 -1.233 -1.147 Cross-Section 0 0 0.1 -1.150 -1.094 10-6 T - 1 = -1 -1.000004 2 5 10 106 0.796 0.8407 2.012 2.275 2.821 3.327 4.187 5.280 0.867 0.887 2.525 2.657 3.930 4.218 7.131 7.960 0.911 0.922 2.924 3.005 4.987 5.174 11.578 12.301 T - 1 = 1 0.999996 T - 1 = 4 3.999929 T - 1 = 9 8.999640 228,223 Top line of each row is from Senior (1975). 1 a diagonal element of the polarization tensor. *he quantity P 11 denotes

-51 - Table 4.2 T 0 P /V for Plates with Square Cross-Section 11 thickness/length 0.5 0.2 0.1 -1.470 -1.271 -1.190 -1.226 -1.139 -1.088 106 T - 1 = -1 -1.000004 2 0.805 0.850 5 2.104 2.381 0.871 0.896 2.593 2.770 4.146 4.566 8.164 9.793 0.908 0.929 2.944 3.099 5.125 5.514 13.127 15.303 T - 1 = 1 0.999996 T -1 = 4 3.999936 T - 1 = 9 8.999678 261,829 10 3.029 3.608 106 4.763 6.300 Top line of each row is from Herrick (1976). The quantity a diagonal element of the polarization tensor. P denotes 11

-52 - Table 4.3 P /V for Plates with Square Cross-Section 33 thickness/length 0.5 0.2 0.1 -2.256 -3.872 -6.43S -1.945 -3.141 -4.738 T 0 106 (T - 1)/T=+ 1,096,342 2 0.673 0.682 5 1.402 1.428 0.592 0.612 1.086 1.164 1.287 1.423 1.524 1.804 0.552 0.582 0.951 1.067 1.101 1.282 1.262 1.569 (T-1)/T = 0.5 0.554 (T-l)/T = 0.8 0.983 (T-1)/T = 0.9 1.161 (T-1)/T = 1 1.375 10 1.763 1.825 10 2.258 2.522 Top line of each row is from Herrick (1976). The quantity P denotes a diagonal element of the polarization tensor. 33

-53 - indicates that the elements may not be sufficiently fine. This hypothesis is supported by the results for the rectangle (presented in the following section), which was generated using a graded mesh which is most fine near the edges. 4.2 A Linear Predictor Model A quick inspection of the data for the circular and square cylinders presented in Tables 4.1 and 4.2 show that the dipole moment appears to be largely determined by the permittivity, T, and the thickness to length ratio. It was felt that a model based only on gross parameters such as thickness, width, length, and T would be able to generate a rough value for the dipole moment associated with a particle. Such a model would be useful for design purposes to get an approximate idea of the dipole moment associated with a particular particle. The model would be restricted to convex shapes, since an object consisting of needle like protrusions would obviously have a much different dipole moment than a convex object with similar thickness, length, width, and T parameterization. The model was constructed by restricting f(x,y) in Eq. (4.1) to be linear. For the z excitation this is equivalent to restricting f(x,y) to be constant. For either tangential or normal excitation this then permits the parameterization of the problem in terms of one unknown. The problem solved is that of a rhombus symmetric about both axes. The reduction of the problem to a linearly varying potential permits much information to be obtained from recasting the problem as an eigenvalue problem in T. After the eigenvalue is obtained, the dipole moment is then known for a given shape for all values of T. For the x excitation the problem is solved by first examining Eq. (3.4)

-54 - t/2 -x + (1 - ) V'G dsI dz' and the particle is assumed to be a rhombus as shown in Fig. 4.1. If pot is then restricted to be linear, from symmetry t consideration pt may be written as t = ax and Eq. (3.4) reduces to t/2 ax = -x + (1 - T) J f a ds' dz'. (4.2) -t/2 The choice of testing function now becomes more important than in the previous multiple element formulation of the problem. For simplicity and to accommodate the existing program the testing function was chosen as a single delta function, centered at x = Q, y = 0, and z = t/2. Then Eq. (4.2) reduces to t/2 a = -a + (1 - T) J" f a G ds' dz' -t/2 Letting t/2 I = r aG ds' dz' -t/2 S yields

-55 - z y thickn thickness width Fig. 4.1: Diagram of a Plate Used by Linear Predictor.

-56 - a: = -k + (1 - T) aI a(k + (- - 1)I) = - ((I-1 ) + ) a: I It -; -I and if the above equation is considered as an eigenvalue problem in T, X - a A A = 1 -- (T - x)a = - - - 1 a I T - T x and since t q = ax the potential on the surface of the plate is known for all values of T after x and I have been computed. For a z excitation, the simplification begins with Eq. (3.5) t/2 t 00 ot = - + (l - T) -t/2 S Gt ~ + a G )ds' dz' tO ' VS a Z Again restricting po to be linear and consistent with the symmetry of the plate the equation yields t z ~ t

-57 - so Eq. (3.5) becomes z a = u -z + (1 - T) t/2 I f -t/2 Js a DG ds' dz', t 2z' (4.3) This time the testing function is chosen as a delta function centered at z = t/2, x = y = O. Let t/2 I-= -t/2 I G ds' dz' az s then Eq. (4.3) becomes a 2 t 2 a T - T) I f and a + T I - ) L/ t 2 t2 2I a (( 1) + so = 1 - 2i a(T - x) t2 2I t2 a = -= 2 2I 1 T - X t o Z = at and the potential is known for all values of T after x and I have been calculated.

-58 - The linear predictor described above was implemented and tested on particles with length to width ratios of one to one and two to one, and with length to thickness ratios of 1, 10, and 100, for x, y and z excitation. For length to width ratios of one to one, the particles tested were the circular and square cylinders (shown in Table 4.4 and 4.5), and they were tested at T = 2, 11, and 101. For the two to one length to width ratio the particles tested were cylinders with cross sections of a triangle, semicircle, and rectangle. These shapes were tested at T = 2, 11, and 101 (Tables 4.6, 4.7, and 4.8) and T = 2 + i2, 11 + i2, and 101 + i2 (Tables 4.9, 4.10, and 4.11). One difference in notation between Tables 4.2 and 4.3 and Tables 4.4 and 4.5 must be noted. For Tables 4.2 and 4.3 the length (= width) was measured along the edge of the square while for Tables 4.4 and 4.5 the length (= width) was measured along the diagonal. This difference in notation was required to match the reference data available for the square in Table 4.2 and 4.3, and to be consistent with the linear predictor model presented in Tables 4.4 and 4.5. Although the effect is a change in the thickness by a factor of V't, this difference had a negligible effect on the resultant dipole moments. As can be seen, the results obtained from the linear predictor are in fairly good agreement with those obtained from the multi-element formulation of Chapter III. The results are particularly accurate for dipole moments near T - 1 for the tangential excitations and (T - 1)/T for the normal excitation. This is as expected and is discussed in detail in the next section. The linear predictor

-59 - Table 4.4 P /V for Several Convex Plates with Length/Width = 1 11 T = 2 11 101 T-1 = 1 T-1 = 10 T-1 = 100 S/t = 1 Linear predictor 0.858 3.77 5.72 14-pt square 0.817 3.21 4.65 12-pt circle 0.817 3.09 4.28 Sphere-Exact 0.750 2.31 2.91 z/t = 10 Linear predictor 0.941 6.15 13.8 14 pt. square 0.913 5.33 10.7 12 pt circle 0.922 5.49 10.9 Q/t = 100 Linear predictor 0.988 8.97 46.7 14 pt square 0.982 8.54 39.9 12 pt circle 0.985 8.73 42.3

-60 - Table 4.5 P /V for Several Convex Plates with Length/Width = 1 3 3 -r= 2 11 101 (T-1)/t=1/2 (T-1)/T=0.909 (T-1)/T=0.990 A/t = 1 Linear predictor 0.757 2.37 3.02 14-pt square 0.800 2.94 4.31 12-pt circle 0.771 2.59 3.59 Sphere-Exact 0.750 2.31 2.91 9/t = 10 Linear predictor 0.533 1.02 1.13 14 pt. square 0.595 1.36 1.62 12 pt circle 0.577 1.25 1.45 Z/t = 100 Linear predictor 0.503 0.919 1.002 14 pt square 0.557 1.192 1.371 12 pt circle 0.544 1.121 1.267

-61 - Table 4.6 P /V for Several Convex Plates with Length/Width = 2 11 2 T-1 = 1 11 T-1 = 10 101 T-1 = 100 k/t = 1 Linear predictor 14 pt triangle 12 pt half-circle 20 pt rectangle k/t = 10 Linear predictor 14 pt triangle 12 pt half circle 20 pt rectangle k/t = 100 Linear predictor 14 pt triangle 12 pt half circle 20 pt rectangle 0.912 0.865 0.866 0.864 0.963 0.930 0.939 0.939 0.992 0.984 0.987 0.987 5.11 4.33 4.06 4.03 9.46 7.96 6.59 6.51 7.24 6.21 6.26 6.36 9.33 8.78 8.95 9.00 20.8 16.3 15.0 15.8 58.8 48.6 49.5 53.4

-62 - Table 4.7 P /V for Several Convex Plates with Length/Width = 2 22 = 2 11 101 T-1 = -1 T = 10 T-1 = 100 z/t = 1 Linear predictor 0.794 2.78 3.72 14 pt triangle 0.772 2.73 3.90 12 pt half-circle 0.768 2.55 3.38 20 pt rectangle 0.768 2.57 3.43 k/t = 10 Linear predictor 0.886 4.38 7.23 14 pt triangle 0.858 4.09 7.29 12 pt half circle 0.870 4.13 6.80 20 pt rectangle 0.879 4.34 7.40 Z/t = 100 Linear predictor 0.974 7.94 27.8 14 pt triangle 0.963 7.44 26.5 12 pt half circle 0.970 7.73 26.8 20 pt rectangle 0.973 8.03 31.8

-63 - Table 4.8 P /V for Several Convex Plates with Length/Width = 2 33 T = 2 11 101 (T-1)/T=1/2 (T-1)/T=0.909(T-1)/T=0.990 /t = 1 Linear predictor 0.809 2.98 4.07 14 pt triangle 0.851 3.74 6.15 12 pt half circle 0.826 3.29 5.01 20 pt rectangle 0.802 2.94 4.20 Z/t = 10 Linear predictor 0.550 1.09 1.21 14 pt triangle 0.641 1.63 2.04 12 pt half circle 0.621 1.50 1.81 20 pt rectangle 0.582 1.26 1.46 k/t = 100 Linear predictor 0.504 0.925 1.01 14 pt triangle 0.594 1.38 1.66 12 pt half circle 0.583 1.32 1.54 20 pt rectangle 0.526 1.02 1.13

-64 - Table 4.9 P /V for Several Convex Plates with Length/Width = 2 11 T - T-1 = and Complex Permittivity 2+i2 ll+i2 1+i2 10+i2 k/t = 1 Linear predictor 14,14 pt triangle 12,12 pt half circle 20 pt rectangle k/t = 10 Linear predictor 14 pt triangle 12 pt half circle 20 pt rectangle Q/t = 100 Linear predictor 14 pt triangle 12 pt half circle 20 pt rectangle 1.19+il.61 1.1 9+il.41 1.22+il.40 1.22+il.39 1.09+i 1.84 1.13+il.69 1.13+il.73 1.13+il.73 1.02+il.97 1.04+il.93 1.03+i 1.94 1.03+il 1.94 5.16+i.517 4.37+i.413 4.10+i.337 4.07+i.333 7.30+il.04 6.26+i.820 6.31+i.798 6.42+i.833 9.36+il.71 8.81+il.56 8.98+il.61 9.03+il.64 101+i2 100+i2 9.46+i.017 7.96+i.015 6.59+i.009 6.51 +i.009 20.8+i.086 16.3+i.063 1 5.0+i.047 15.8+i.053 58.5+i.686 48.6+i.521 49.5+i.510 53.4+i.601

-65 - Table 4.10 P /V for Several Convex Plates with Length/Width = 2 22 and Complex Permittivity T = T-1 = 2+i2 l+i2 11+i2 10+i2 101+i2 100+i2 Q/t = 1 Linear predictor 14 pt triangle 12 pt half circle 20 pt rectangle z/t = 10 Linear predictor 14 pt triangle 12 pt half circle 20 pt rectangle z/t = 100 Linear predictor 14 pt triangle 12 pt half circle 20 pt rectangle 1.23+il.07 1.18+il.00 1.20+i.977 1.20+i.982 1.22+il.49 1.21+il.37 1.22+i 1.42 1.21+il.46 1.07+il.89 1.09+i 1.84 1.07+il.87 1.06+i 1.88 2.81+i.152 2.75+i167 2.57+i.133 2.59+i.137 4.42+i.379 4.13+i.368 4.17+i.348 4.29+i.386 7.99+il.26 7.49+il.13 7.78+il.20 8.08+il.31 3.72+i.002 3.90+i.004 3.38+i.002 3.43+i.002 7.23+i.010 7.29+i.014 6.80+i.010 7.40+i.011 27.8+i.155 26.5+i.170 26.8+i.153 31.8+i.216

-66 - Table 4.11 P /V for Several Convex Plates with Length/Width = 2 33 and Complex Permittivity (T-1)/T = k/t = 1 Linear predictor 14 pt triangle 12 pt half circle 20 pt rectangle Z/t = 10 Linear predictor length/width 2+i2 0.75+i.25 1.24+il.14 1.23+i.133 1.24+il.21 1.23+il.11 0.852+i.335 1.009+i.571 0.978+i.517 0.911+i.411 0.759+i.257 0.926+i.468 0.907+i.442 0.801+i.303 ll+i2 0.912+i.016 3.00+i.174 3.77+i.291 3.32+i.224 2.96+i.178 1.09+i.023 1.64+i.063 1.51+i.052 1.26+i.034 0.928+i.016 1.39+i.045 1.32+i.039 1.02+i.021 101+i2 0.990+i.0002 4.07+i.003 6.15+i.011 5.01+i.007 4.20+i.005 1.21+i.000 2.04+i.001 1.81+i.001 1.46+i.000 1.01+i.000 1.65+i.001 1.54+i.000 1.13+i.000 14 12 pt pt triangle half circle rectangle 100 20 pt k/t = Linear predictor 14 pt triangle 12 pt half circle 20 pt rectangle

-67 - uniformly underestimates the dipole moment calculated for z excitation (for Real T > 1). This is a result of choosing the testing function as a delta function centered on the middle of the plate thus eliminating any effects of the "buckling" near the edges of the plate. This difference between the linear predictor and the multi-element algorithm is in the opposite direction of the difference between the multi-element algorithm and the reference data, thus the result of the linear predictor for the z excitation might be better than Tables 4.5, 4.8, and 4.11 would indicate. Also, the dipole moments calculated for the rectangle for z excitation are lower than those calculated for the triangle and semicircle with the same thickness and T. This supports the idea that the multi-element formulation for the z excitation might be using elements which are too coarse to accurately represent the buckling near the edges. Finally, it is interesting to examine the eigenvalues extracted by the linear predictor algorithm, and these are shown in Table 4.12. These seem to parallel the path of the eigenvalues associated with the spheroid (e.g., Senior and Weil, 1982) and would seem to have physical significance. Since there is only a single eigenvalue the interpretation is simply as an estimate for the location of the resonance region. The most striking result is that the eigenvalue for the z excitation, though negative as expected, is extremely close to zero, while the x and y components are much farther from zero. From the perspective of designing particles with high absorption (operating in the resonance region), the important part of the dipole moment may be the z component. This is in contradistinction to the region

-68 - Table 4.12 Eigenvalues Extracted by Linear Predictor length/width = 1 t/z = 1.0 t/k = 0.1 t/k = 0.01 t/k = 1.0 t/k = 0.1 t/k = 0.01 x excitation -5.07 -15.02 -86.62 z excitation -2.12 -0.14 -0.01 length/width = 2 x excitation y excitation -9.45 -2.86 -25.32 -6.80 -140.44 -37.66 z excitation -3.25 -0.22 -0.01

-69 - of Real:> O, where the normal component is usually negligible compared to the tangential components. 4.3 Shape Effects As may be seen from Tables 4.4 through 4.11, the shape of a plate does have an effect on its associated dipole moments. Shape is important when the product of thickness/length (t/z) and T is large. For (Tt/z) small the tangential components of the dipole moment may be approximated by T-l, and the normal component may be approximated by (T-1)/T. The physical significance of these approximations for the tangential excitation is a total electric field Et which equals the incident electric field E. For the normal excitation, the total electric field Dt is assumed to equal the incident electric field D1. These results may be obtained by examining the formulae for the dipole moment. From Senior (1976) X = (1 - n * xi t ds' (4.4) B starting with x. = x for tangential excitation. 1 i+ s = _x + s = -x + (on the plate) X = (T - 1) fn. x ds' + (1 - T) x B B = (r - 1) J v x dv' + (1 - T) (xS) dv' (cont.) V V 1 B V

-70 - =(T - l)vB + (1 - T1) BV 4S p1 DXdv' So, for S ~X <<1 xl z (T1-I)V For normal excitation,. x. = Z 1 ~. t -I O + 0p~ 1 T 1 1 Tr 3 on the plate from (4.4) = — i z ds ' + (1 -T) nJ>Iz2 B B ds' T -1 V * 7 dv' + B (V T) f V, ) dv' B T- 1r - T VB + (V T)Jf VB dv' So, for IT << 1 T-1 V 33~ T vB

-71 - Thus both approximations assume the scattered field on the plate is negligible compared to the incident field. The motivation for the linear predictor model was to account for small but non-zero scattered fields. Since the incident field would still be the dominant component, the total potential on the plate should still be effectively linear. The linear predictor provides a next order approximation to the scattered field. Instead of assuming the scattered field is negligible, it is linear. This permits an increase in accuracy over the small scattered field approximation and incorporates gross shape effects. However, as the scattered field increases the importance of the non-linear components of the scattered field also increases. An indication of the size of the scattered field is the amount by which the computed dipole moment (using the linear predictor model) differs from the small scattered field model (i.e., T - 1 or (T-1)/T). The primary benefits of the linear predictor model then are to give a rough estimate of the dipole moment, and when such estimates differ from the small field estimates of the dipole moment to indicate that use of the multi-element formulation is appropriate. The region where the use of the multi-element formulation is appropriate constitutes a very important region, particularly in regard to the analysis of aerosols. Particles which remain suspended in the air generally have a large surface area relative to their weight. Thin particles consitute a common realization of this characteristic. In the case of aerosols which efficiently absorb or reflect radio waves, the constituent particles must either have a relatively large permittivity or a permittivity in the resonance region of the particle. The dipole moments

-72 -assocated with both of these cases are strongly influenced by the shape of the particle, and this may be further investigated through the use of the multi-element program developed in Chapter III.

CHAPTER V. SCATTERING BY DISTRIBUTIONS OF PARTICLES 5.1 The Clausius-Mosotti-Lorentz-Lorenz Equation The problem of scattering by sparse distributions of particles is now considered. Although the formulae used are classical (Mosotti, 1850), various modifications to the standard form have surfaced to describe experimental results for distributions with high densities. In this chapter the Clausius-Mosotti-Lorentz-Lorenz equation will be derived, and it will be shown that although modifications to this formula may approximately describe experimental results, such extensions are not rigorous and hence should not be expected to work in cases other than those for which they were developed. The final section discusses scattering by plates as formulated in the two preceding chapters, and the validity of a limit case is shown. Assume a uniform applied field, E, in a distribution of particles. For each particle, a dipole moment V is generated, where p = aE where a is the polarizability of the particle, and EP is the local electric field acting on the particle. E differs from E as a result of near-field disturbances of nearby particles. As the distribution of particles becomes sparse, E will converge to E. -73 -

-74 - Now, if the density of particles is N per unit volume, then P = Np = NaEE However, from the definition of P. P = D- OE E ( ) ECo = r 1) To find a relation between the polarizability per unit volume, No, and the relative permittivity, r, assume a capacitor as shown in Fig. 5.1. The distribution of particles is now located between a pair of infinite parallel plates. The region between the plates is uniformly filled by the distribution of particles. The voltage between the plates is V = d|E|, E = E(-z) And, E (A) = E (A)+ E (A)+ E (A), 1 2 3 where, E (A) is the local electric field acting upon the particle located at "A". E (A) is the electric field resulting from the charges located on the metal plates. E (A) is the electric field resulting from the polarized 2 particle distribution exterior to the sphere "S". The sphere "S"

-75 - V=dIEl A. 9.4 d x ZI v=o Fig. 5.1: Parallel Plate Capacitor with Spherical Exclusion Volume.

is taken large enough so that the near field effect of the individual polarized particles is not apparent when viewed from "A". E (A) is the electric field resulting from the polarized 3 particle distribution interior to the sphere "S". In general, it is not possible, or at least not easy, to calculate E. For a 3 cubical array of particles (i.e., particles lying on the vertices of adjoining cubes) E = 0, as shown in the following section. The Mosotti (1850) approximation is to simply set E = 0. 3 The next step is to calculate E and E. 1 2 E = P (-z), p - charge density on top plate 1 o E has two components. The first component is due to a bound charge 2 distribution in the vicinity of the plates. The second is the bound charge distribution surrounding the sphere. Since the sphere is large enough to hide individual particle effects, it is permissible to consider the distribution of particles exterior to the sphere as a dielectric continuum and use terms such as "bound charge". So E = b( ( - zPb2 ds 2 C ) 4b c R2 S 0 Then, by linear superposition,

-77 - E = E+ E ~ 1 2 = (-) + (-z)- cp ds 0 o S 47FE R2 b ds 0 0 CO 4TWc R2 = + Pbl (__) + i' ds Since v * P = -Pb1, and further, since the effects of the polarized particles within the sphere are accounted for with E, the 3 dielectric representation of the polarized particles external to the sphere may be assumed to terminate on the surface of the sphere. Hence, P 0 within the sphere. Now, despite the segmentation of the particle distribution with the sphere, the distribution is, in fact, uniform. Therefore, P = IPI (-z) and P n = -IPI cos e so = P flPicos e ds (-z) ds o 0 4scoR2 But, by the symmetries of the problem, EQ has only a z component. Therefore

-78 - E = bi (-z) 0 bi(-Z) 0 0 P + %ib-z 0 P + Pb 0 - I1 cos2 e J) ds z 4Tr c R2 0 0 Cos 2 a R2 sin e de d~ 2c0 (2TV) [- ~-cos2 ~ 0 2 F IL 0 and + P~bl 0 so p - E + 3 D - cE D - E +360E = E + - 0 30 3 0 D+2sE ~E- C~ 3 0 3 0 ~ + 2 = p r3 Thus P = NctE z = N 3 but P = D - E = c - c ) o =E[Cr(c - 1)

-79 - So E (c - 1) Co(~r - 1) Cr + 2 = NaE 3 r + 2 + Na 3 - and rearranging C - 1 Nu = 3c whh is the rel on des red which is the relation desired. Nc 3 co Cr - 1 cr +2 is the standard form of writing the Clausius-Mosotti-Lorentz-Lorenz equation. Solving for c, r ( r + 2) 3NC = r r - 1 CrN( ~r 3 c 'r ) 1) = -1 - 2 Na3 3so 1+2 Na 1 N3o - No 3c = 1 +- Na 3co

-80 - = 1+ Nu I Nu 5.2 Cubical Distributions of Particles As noted in the previous section the exclusion volume is a virtual cavity which has no physical significance, and is merely a computational artifice. The exclusion volume must satisfy two constraints. First, all particles outside the exclusion volume must be far enough away to be accurately represented by an equivalent, homogeneous, dielectric. Second, the combined effect of the particles within the exclusion volume should not alter the electric field at the point of computation. The first constraint may be satisfied immediately by letting the exclusion volume be arbitrarily large. The second constraint is not usually possible to satisfy exactly for an arbitrary particle distribution, so that an approximate solution is used which works reasonably well. Since the task of determining the field contributions of an arbitrary particle distribution is in general very difficult, a highly symmetric particle distribution is used. The standard reasoning is that if the field contribution of a symmetric particle distribution is zero, then the field contribution of a random distribution of particles should also be zero. The derivation for a cubical distribution is given below and follows the approach of Jackson (1975).

-81 - The field resulting from a dipole is E 3n(p n) - Ix - xo3 where p = the dipole moment, x = the location of the dipole, x = the point of field measurement and n = points from x to x. Then the total E field due to a cubical array of dipoles, each with identical p, is E =2 ijk 3(xi,,z k)[ (xi,yzk)] - P 2 ij k |jk To show E = E = E = 0, write x y z EX ijk 3x.[p * (xi yj Zk)] - P I 12 I~15 The summation is taken over all vertices of the array within the exclusion volume, and (x.,yj,zk) are the coordinates of each dipole. E2 3xp 3 p + 3x + ip - Pxl12 i - ix 1 5 ijk lXI

-82 - Since this equation must hold for all choices of p xpY pz, then x iyj xiz E E —= 0 -1 | |5 Z- y|5 ijk l1 ijk l and from similar formulae for E = Ez= 0,. we get yzk E YiZk ijk |x|5 This represents a rather strong symmetry constraint, and is satisfied by the sphere, the cube, and an infinity of other symmetrical shapes. However, since the only purpose of the exclusion volume is to yield the Clausius-Mosotti equation, which is shape independent, there is no point in finding all the acceptable shapes. In the following section the computations are performed using a cubical exclusion volume. 5.3 A Cubical Exclusion Volume The standard derivation of the Clausius-Mosotti relation used a spherical virtual cavity (see Section 5.1). However, the only reason for using the sphere is a general symmetry argument, and a cubical cavity should work equally well. If the spherical cavity is replaced by the cubical cavity in the Clausius-Mosotti derivation, two of the field components

-83 - may be affected. The first is the field due to the particles within the cavity. However, due to symmetry considerations (see the preceding section) this field component is zero forboth the spherical cavity and the cubical cavity. The second field component is due to the bound charge on the virtual cavity surface. For a sphere, the resulting field was shown in Section 5.1 to be - I 1 E = E where E is the uniform macroscopic field. It will now be shown that the cube yields the same result. Assume a virtual cube in a uniform electric field E = Eo(-z) similar to Fig. 5.1. Then the bound charge, Pb = -v * P, will accumulate on the top and bottom surfaces. We wish to find the E field due to this charge distribution at the center of the cube. -P E = - Rds, E t 47wcR2 where R is measured from the center of the cube. E v' PR ds s 4wrR2 A ^ = + R ds + m Rds top 47rR2 bottom 4 where P denotes a "macroscopic" P. Since P - 0 inside the cube,

-84 - top - m R ds - 27TrR2 Further, by symmetry: EQz Z z top - R*z ds = 2rER2 Il l ^ - 2F z 271 f ds top c -1 = + r 2Tr top z R3 ds q For a sphere, E = - E So, if 2Tr 3 top z R3 ds then the field in a cubical virtual cavity is the same as that in a spherical virtual cavity. Let the dimensions of the cube be 2. x 2z x 2z. Then rz ds top R3 fZ f dx dy - R3 z= *I Q -4: 2 + y2 + Z2)3/2dx dy _ (x2 + y2 + z2)3/2 2=k

-85 - Then, from Gradshteyn and Ryzhik (1980, Sec. 2.271.5) -9PI z x y2+ z2 V1X2+ —2 + Z2 dy Z=, 9, = 1j 2 1 dy 9, 0 49 21, ~~ 1 -dy 2+ Z2 Vy=+-2 97' Let = 2 - Y2 + 2z2. 2 =2 y2 = uz -22 U = - Z 2 1 V-y +2 k dy = k(u2 - 2) 1/2 u du Then top Z ds u 2z2 - 9,2 1 -U v'u27 7, du Vu2 ~- 2 4 1 du Expanding in partial fractions yields - xf + l 7 (x - ])A

-86 - l et u = x + 1 and v = x - 1 = 2.~ ~~ -1 - du + f 2vl '2 vvV —+2 civ Then, again from Gradshteyn and Ryzhik (1980, Sec. 2.266) = 2~sin1(-2 -2x) + sin- (-2+ 2x) Si 2.f 1i ( -2 - 2(43 + i (4~-+ 1)4V8 + sin -1 (-2 - 2'~ 1)) - sinl(1 2 + 2 /2- - ( Y/2 -1I)4v8 + sin 1 (-2 + 2(4 3-1) (4~ —1I)4A =2{ s in-1 F 5 7r = < 2~ (-42 - )3 + sini + sin -'(-1 )-sin -1 _) - 2(j) 2gLL Thus the spherical and cubical exclusion volumes yield the same resul t.

-87 - 5.4 Restriction of the Clausius-Mosotti Equation to Distribution of Plates The Clausius-Mosotti relation is valid for any sparse distribution of particles. Thus = 1 + Na Na co1 - 3 may be used with a obtained from the algorithms presented in Chapters III and IV. The validity of the formula can be verified for very thin plates in a dense cubical distribution. a = ( - 1)v = (T - 1) tk2 for excitation parallel to the plate and T-1 T-l 2 = = I t v2 T T for excitation normal to the plate, where the particle is assumed square with dimensions z x z x t. The above formulae were obtained in Chapter IV with a small scattered field approximation which is satisfied as long as Tt/k <<1. If the cubical distribution of plates is now made dense so that the plates form a multiple layer dielectric, then r as seen under tangential excitation would be t r

-88 - and for normal excitation (11 For parae excitation, r For parallel excitation, Na = N(T - 1)t2 = c = 1 +3(T - 1)t/ r 3-(T -1)t/z t (T - 1) t 1 + (T - 1 )t/ = (1 - t/z) + T(t/Z) as expected. For normal excitation C = 1 + r 3 - 1 t T - 3T -1 t T Z 3- 1 t+ 3 - t T k T z 3 - T - 1 IF t 9 3 —1 t 1 -3 T Z j - 3+2 T 1- t T j

-89 - [(l - t) (3 + 2 - r ) + Q (3+2T-l t) 3T1 -1 3 + 2T -1 t T Z J 3 + 2 - 1 t__ = + T 3 3+2T - t T t k (_ -t + t 1 3 + 2(T - 1)t/ 1 T 3+2 T 1 +t 1 as expected. Thus the effective permittivity of a dense cubical distribution of thin plates as predicted by the Clausius-Mosotti theory is in agreement with the exact average permittivity of a layered dielectric when excited either normally or tangentially to the plane of the layers. Further, the Clausius-Mosotti theory is entirely adequate for describing scattering by distributions of particles as long as the electrical interaction between the particles may be accurately described with only the dipole term. The low frequency expansion permits the Clausius-Mosotti formula to be used with dynamic

-90 - electromagnetic radiation as well as for the electrostatic and magnetostatic problems. For dynamic radiation, the dipole moment may contain an imaginary component which then results in an imaginary component in the effective permittivity. On a macroscopic level this corresponds to a lossy medium. Since these results can also be applied to magnetic scattering, the entire distribution can then be characterized by complex electric and magnetic polarization tensors, Although the analysis of electromagnetic propagation through a complex tensor medium is not particularly simple, it is trivial compared to the problem of analyzing propagation through a distribution of particles by computing the scattering of the electromagnetic wave from each individual particle in the distribution.

CHAPTER VI. CONCLUSIONS The problem of constructing the low frequency expansion which is valid for scattering by an open surface has been solved with the simultaneous solution of a problem from classical physics, find F given v x F = f. A highly efficient and accurate numerical program has been developed to solve the problem of static scattering from a thin plate. The potential on flat plates with thickness to length ratios less than 0.1 and arbitrary shape may be solved and the resulting dipole moments computed, with the plate being described via an arbitrary number of triangular elements with arbitrary shape and size. This flexibility permits the description of an arbitrary plate with a minimum number of elements. Varying the size of the elements is important near the edges of a plate, primarily when computing the dipole moments associated with normal excitation of the plate. The range of operation of the program is very broad. The permittivity can assume a wide range of values, real and complex, small and large. The result is that the scattering from a particle composed of any homogeneous material can be analyzed, including perfect conductors. The thickness of the material can also assume a wide range of values, and the program has been successfully tested for length to thickness ratios ranging from 1:1 to in excess of 106:1. Finally, distributions of particles have been discussed, and the completeness of the Clausius-Mosotti-Lorentz-Lorenz formulation has been -91 -

-92 - shown. The validity of this formulation for a cubical dense distribution of thin plates has also been shown. Although the content of Chapter V was restricted to sparse (particle volume as a percentage of total volume) distributions of particles, an area of considerable interest is dense distributions of particles. The primary difference between dense and non-dense scattering theories is that sparse scattering theories assume the particle interaction is limited to far-field interaction. In Section 5.4 the thin plate approximation for the dipole moment guaranteed that the scattered field was negligible and that only far field interactions were being considered. There have been several attempts to modify the ClausiusMosotti-Lorentz-Lorenz formula to account for non-sparse distributions. The theories are designed to account for experimental results on dense distributions of particles and in general they may be characterized as attempting to describe higher order interactions (quadrupole, etc.) without specifically calculating these components. Discussions of the relative merits of these methods are contained in Granqvist and Hunderi (1977) and Bohren and Battan (1980). Future work may use the program developed for this dissertation to investigate the resonance region associated with thin plates. To this end, the matrix problem solved might be decomposed into an eigenfunction expansion. This would permit the simultaneous calculation of the dipole moments associated with a particular shape for all values of T. The linear predictor described in Chapter IV, although quite satisfactory for the purpose of determining when to use the full multi-element program, could perhaps be further refined by modifying the testing function. Specifically, it might be desirable

-93 - to use a Galerkin's method approach (e.g., Harrington, 1982). Although the results presented in Chapter IV are believed accurate and show qualitative agreement with the results of Senior (1975) and Herrick (1976), it would be useful to verify some of the thin plate calculations for which the program was developed (length to thickness ratios greater than 10), either experimentally or numerically. Efforts are currently underway to obtain verification of the results using a version of the program developed by Senior and Willis (1982) specifically modified to permit the accurate analysis of thin rotationally symmetric, plates (Weil, 1984). It is hoped that the developed program (a listing of which is provided in the Appendix will be incorporated into emerging theories on low frequency scattering by distributions of particles. Since the shape of a particle affects the dipole moments which are associated with a particle which in turn determines parameters such as the effective dielectric constant of the distribution, the developed program should provide a valuable tool in determining the electrical properties of particle distributions. As shape becomes extremely important in the resonance region, the shape of the particle may aid in the identification of different particles in cases where the particles are illuminated with relatively low frequency electromagnetic radiation. This may be particularly useful in regions where identification through high frequency electromagnetic radiation, which relies primarily on shape effects which become apparent only when the particle size is comparable to a wavelength, is either not possible or inconvenient.

APPENDIX PROGRAM LISTING 1 /* function definitions */ 2 #include <math.h> 3 double ppcon(), ptcon(),u(),v(),ur(),vr(),intdvz(),intdz(),msqrt(),mlog(). 4 area() sumang(); 6 /* macro definitions */ 6 #define Darray(Ara,I1,I2.L1.L2) Ara[(Il)+(I2)*(L1)] /* Macro to generate other\ 7 macros to mimic fortran\ 8 arrays */ 9 #define Mnumpoi 100 /* maximum number of points */ 10 #define Mnumtri 100 /* maximum number of triangular patches */ 11 #define Mnumatri 10 /* maximum number of triangles which may share a \ 12 common point, 6 is average for internal points, 8 accounts for\ 13 most schemes, thus 10 is a reasonable upper bound */ 14 #define True 1 /* value of logical true */ 15 #define False 0 /* value of logical false */ 16 #define Tri(J1,J2) Darray(tri,J1,J2,Mnumatri,Mnumpoi) /* set up tri as two\ 17 dimensional array with limits (Mnumatri, Mnumpoi) */ 18 #define Poil(Jl,J2) Darray(poil,Jl,J2,Mnumatri,Mnumpoi) /* set up poll as two\ Is dimensional array with limits (Mnumatri, Mnumpoi) */ 20 #define Poi2(Jl,J2) Darray(poi2,Jl,J2,Mnumatri,Mnumpoi) /* set up poi2 as two\ 21 dimensional array with limits (Mnumatri, Mnumpoi) */ 22 #define Poi(J1,J2) Darray(poi,Jl,J2,Mnumtri.3) /* set up poi as two dimensional\ 23 array with limits (Mnumtri,3) */ 24 #define Matrix(Jl,J2) Darray(matrix,J1,J2.Mnumpoi,Mnumpoi) /* set up matrix as\ 25 two dimensional array with limits (Mnumpoi,Mnumpoi) */ 26 #define Cmatrix(Jl,J2) Darray(cmatrix,J1,J2,Mnumpoi,Mnumpoi) /* complex 27 version of Matrix, used with the complex 28 structure */ 29 struct complex { 30 double real; 31 double imag; 32 33 /* */ 34 /* */ 35 /* End of Definitions / Beginning of external variable declarations 36 / */

37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 /* int flagsym; / int flagsyms; Int plotsym; int eof; /* en int potflag; / */ * flagsym may assume values of 0,1,2, or 3. Flagsym is used by low level routines to incorporate symmetry in a manner invisible to the rest of the program. 0 - indicates no symmetry 1 - indicates object possesses mirror symmetry about x=O 2 - indicates object possesses mirror symmetry about y=O 3 - indicates object possesses mirror symmetry about x=O and y=O Flagsym is set at the beginning of the main program and after that is never changed */ /* flagsyms is used in conjunction with flagsym to generate mirror images of the source patches while calculating contributions to the field point. The initial object description is assumed to lie in the first quadrant, however if flagsym is 0,. this is not necessary. Flagsyms is always less than or equal to flagsym, and may assume the values of 0,1,2, or 3. These values are given the following meanings: 0 - source patch is original patch (presumably first quadrant) 1 - source patch is in second quadrant 2 - source patch is in fourth quadrant 3 - source patch is in third quadrant */ /* plotsym is used to denote what type of mirroring is desired in the perspective plot. 0 - no mirroring 1 - mirror about x=0 2 - mirror about y=0 3 - mirror about x=0 and y=0 Any mirroring that is done must be consistent with the mirroring specified in generating the plate. plotsym is specified in rgrph. */ Id of file variable (=0 means end of file) */ '* flag marks whether or not potentials have been computed 0 - potential has not been computed 1 - potential has been computed assuming x excitation 2 - potential has been computed assuming y excitation 3 - potential has been computed assuming z excitation */ x[Mnumpoi]; y[Mnumpoi]; poten[Mnumpoi]; complex cpoten[Mnumpoi]; /* complex version of poten */ /* point is a structure which contains the locations of all of the points used in defining the triangular patches. A point number 0 is permitted as C defines arrays beginning with 0. Q0 struct { double double double struct }point

81 X and y are the coordinates of the points, and 82 poten is the potential of each point, as determined from the 83 matrix problem. */ 84 struct { 85 int numatri[Mnumpoi]; 86 int Tri(Mnumatri,Mnumpoi-1); 87 int Poil(Mnumatri,Mnumpoi-1); 88 int Poi2(Mnumatri,Mnumpoi-1); 89 } apoint; /* apoint is a structure which contains lists of triangles 90 associated with each point. numatri contains the number of 91 triangles associated with each point; tri contains the list of 92 triangle numbers associated with each point; and poil and poi2 93 contain the other two vertices used in defining the triangle. 94 poil and poi2 are indices to point. tri is an index to 95 triang. 96 struct { 97 struct { 98 double x[Mnumtri]; 99 double y Mnumtri]; 100 } centroid; 101 int Poi(Mnumtri,3-1); 102 double poten[Mnumtri]; 103 struct complex cpoten[Mnumtri]; /* complex version of poten */ 104 } triang; /* triang is used to solve the scattering problem with z 105 excitation. centroid.x and.y contain the coordinates 106 of the centroid of each triangle. poi contains the list of 107 points (vertices) associated with each triangle, and is an 108 index to point. poten is the potential of each triangle as 109 determined from solution of the matrix problem. */ 110 double Matrix(Mnumpoi,Mnumpoi-l); /* This is OTHE' matrix. Trivial points 111 (ie, those which have zero potential from 112 symmetry considerations) are skipped when 113 reading or filling the matrix, which is 114 always done sequentially. */ 115 struct complex Cmatrix(Mnumpoi,Mnumpoi-1); /* complex version of Matrix */ 116 double fvect[Mnumpoi]; /* fvect is the forcing vector for the matrix problem 117 and after solution of the matrix problem contains 118 the solution vector. */ 119 struct complex cfvect[Mnumpoi]; /* complex version of fvect */ 120 int mnumpoi; /* This is the total number of points. Points are assumed to be 121 numbered sequentially from 0. mnumpoi <= Mnumpoi. */ 122 int mnumtri; /* This is the total number of triangles. Triangles are assumed 123 to be numbered sequentially from 0. mnumtri <= Mnumtri. */ 124 union {

125 char str[2]; 126 char let; 127 } corn; 128 /* This is the first character of each input line, and is 129 interpreted as a command. Valid commands are: 130 d - display linkup of points and triangles 131 and display potentials if defined 132 e - specify direction of exciting field, 133 and solve the resultant matrix problem. 134 g - graph potentials of points/triangles 135 h - enter heading, ie, a descriptive one line title. 136 m - enter material parameters of plate: t and tau. 137 p - enter coordinates of next point 138 r - regenerate matrix problem using finer mesh 139 s - define symmetry to be assumed in solving problem 140 t - enter definition of next triangle * 141 char hedstr[82]; /* contains one line heading, description of data */ 142 double t; /* Thickness of the plate */ 143 double tau; /* permittivity of the plate */ 144 double taui; /* imaginary part of the permittivity. This variable is also 145 used as a flag: if taui is zero, computations are performed 146 assuming "tau" is purely real; if taui is non-zero, the routines 147 appropriate for a complex tau are invoked */ 148 double dipmom; /* contains the real part of the computed dipole moment */ 149 double dipmomi; /* contains the imaginary part of the computed dipole moment */ 150 double Epsilon = le-10; /* A very small number, used for approximate equality */ 151 double Pi = 3.1415926535; /* Pi */ 152 /******************************************************************************/ 153 /* */ 154 *M A I N PROGRAM */ 155 /* */ 156 /* This is a program which calculates scattering by an arbitrarily shaped, */ 157 /* thin, flat, dielectric plate. Excitation may be specified in the x, y, */ 158 /* or z directions. The plate is made up of an arbitrary number (nominally */ 159 /* less than 50) of triangular patches, upon which a method of moments */ 160 /* solution is obtained. Contributions from each patch is calculated via */ 161 /* surface integrals which are evaluated analytically, thus contributions */ 162 /* from each patch is obtained exactly (at least to 10 plus digits). The */ 163 /* only approximations which are made are that the potential varies linearly */ 164 /* with z inside the plate, and the division of the plate into triangular */ 165 /* patches, inside of which the field varies linearly with x and y. The */ 166 /* shape and size of the triangular patches are completely arbitrary, and it */ 167 /* is noted that the patches need not be contiguous. To speed computations, */ 168 /* and to enhance that accuracy obtainable with a given number of patches, */

169 /* the user may specify that the plate possesses mirror symmetry about x=O, */ 170 /* y=O. or both x=O and y=0. */ 171 /* */ 172 /****************************************************************************/ 173 main() { 174 hedstr[0] = '\0'; 175 for(eof=scanf("%ls",com.str);eof == 1;) 176 switch (com.let) { 177 case P': 178 case 'p': 179 case 'T': 180 case 't': 181 rdata(); /* read in the data */ 182 continue; 183 case 'D': 184 case 'd': 185 ddata(); /* display the data */ 186 break; 187 case 'G': 188 case 'g': 189 rgrph(; /* read in mirroring specification */ cc 190 break; 191 case 'H': 192 case 'h': 193 gethed (); 194 break; 195 case 'M': 196 case 'm': 197 rmp(); /* read in material parameters */ 198 break; 199 case 'S': 200 case 's': 201 rsym(); /* read in the symmetry of the plate */ 202 break; 203 case 'E': 204 eigen(); /* solve eigenvalue problem */ 205 break; 206 case 'e': 207 rexcit(); /* read in direction of excitation 208 and solve the resulting matrix problem */ 209 switch (potflag) { 210 case 1: 211 xsolv(); 212 pexc();

213 break; 214 case 2: 215 ysolv(); 216 pexc() 217 break; 218 case 3: 219 zsolv(); 220 pexc(); 221 break; 222 default: 223 printf("invalid excitations); 224 break; 225 226 break; 227 default: 228 printf("%c is an invalid command\n',com.let); 229 while(scanf("%c",com.str),com.let!= '\n'); 230 /* flush out bad command */ 231 break; 232 } 233 eof=scanf("ls",corn. sr); 234 } 235 236 ***************************************************************************** 237 /* 238 /* ROUTINE: ddata */ 239 /* */ 240 /* ddata is called by the main program to display data. ddata display */ 241 /* locations of points, connections of points, and potentials of points if */ 242 /* they have been computed. ddata also displays status of various flags */ 243 /* which indicate symmetries of plate and direction of excitation. * 244 /* */ 245 ****************************************************************************** 246 ddata() { 247 int i.,; /* loop variables */ 248 if (t!= 0) printf("t = %lf, tau = %lf + i %lf \n'.t,tau,taui); 249 switch (potflag) { 250 case 0: 251 printf("potential has not been computed\n'); 252 break; 253 case 1: 254 printf("potential has been computed assuming x excitation\n'); 255 break; 256 case 2:

257 printf (potential has been computed assuming y excitation\n'); 258 break; 259 case 3: 260 printf (potential has been computed assuming z excitation\n'); 261 break; 262 263 switch (flagsym) { 264 case 0: 265 printf('plate is not symmetric\n'); 266 break; 267 case 1: 268 printf('plate has mirror symmetry about x=0\n'); 269 break; 270 case 2: 271 printf("plate has mirror symmetry about y=O\n"); 272 break; 273 case 3: 274 printf('plate has mirror symmetry about x=0 and y=O\n"); 275 break; 276 > 277 for (i=0; i<mnumpoi; i++) - 278 printf('point # %d is located at x = %lf, y = lf\n'"i,point.x[i] ~ 279 point.y[i]); 280 if (potflag) ( 281 if(!taul) printf ( and has potential %e\n" point.poten[i]); 282 else printf(' and has potential %lf + i %lf\n', 283 point.cpoten[i].real,poinpoint.cpoteni].imag); 284 } 285 printf ('point %d is associated with points and triangles (tri,P1,P2)\n' 286 i); 287 for (j=0; J<apoint.numatri[i];j++) { 288 printf("(%d,%d,%d), ',apoint.Tri(j,i), apoint.Poil(J,i), 289 apoint.Poi2(,i)); 290 if (j%6 = 5 == apoint.numatri[i]-l) printf('\n'); 291 > 292 } 293 294 /*************************************************************************** 295 /* */ 296 /* ROUTINE: rsym */ 297 / */ 298 /* rsym is called by the main program to read in the symmetry of the plate. */ 299 /* rsym assumes no symmetry if response is inappropriate. */ 300 /* */

301 /*******************************************************************/ 302 rsym() { 303 scanf (%d",flagsym); 304 switch(flagsym) { 305 default: flagsym=O; 306 case 1: 307 case 2: 308 case 3: 309 310 } 311 } 312 /******************************************************************************/ 313 /* */ 314 /* ROUTINE: rgrph */ 315 /* */ 316 /* rgraph is called by the main program to read in the mirroring desired. */ 317 /* rgrph then calls graph to produce the graph. plotsym is enforced to */ 318 /* be consistent with flagsym. */ 319 /* */ 320 /******************************************************************************/ 321 rgrph() { 322 scanf ('%d",.plotsym); 323 switch(plotsym) { 324 default: plotsym=0; 325 case 1: 326 case 2: 327 case 3: 328 329 } 330 plotsym = plotsym t flagsym; 331 graph(); 332 > 333 /******************************************************************************/ 334 /* */ 335 /* ROUTINE: rmp */ 336 /* */ 337 /* rmp is called by the main program to read in the material parameters of */ 338 /* the plate, specifically thickness and permittivity. */ 339 /* */ 340 /******************************************************************************/ 341 rmp() { 342 scanf("%lf %lf %lf.,&t,ttau,.taui); 343 dipmomi = 0; 344 }

345 /*****************************************************************************/ 346 /* */ 347 /* ROUTINE: rexcit */ 348 /* */ 349 /* rexcit is called by the main program to read in the direction of */ 350 /* excitation. rexcit assumes x excitation in the event of an invalid */ 351 /* response. */ 352 /* */ 353 /******************************************************************************/ 354 rexcit() 355 char charscr[2]; /* scratch variable, contains direction of excitation */ 356 scanf("%ls",charscr); 357 switch (*charscr) { 358 default: 359 potflag = 0; 360 break; 361 case 'X': 362 case 'x': 363 potflag = 1; 364 break; 365 case 'Y': 366 case 'y': ' 367 potflag = 2; 368 break; 369 case 'Z': 370 case 'z': 371 potflag = 3; 372 break; 373 } 374 } 375 376 /****************************************************************************** 377 /* */ 378 /* ROUTINE: xsolv */ 379 /* */ 380 /* xsolv is called by main to fill THE matrix. xsolv assumes excitation is */ 381 /* in the x direction. Points which lie on x=0 are assumed to have zero */ 382 /* potential and are skipped in both the row and column loops. The outside */ 383 /* loop is the column loop and is associated with the field point. The */ 384 /* inside loop is associated with the source point. Images of the source */ 385 /* point are not apparent at this level and are handled in ppcon. */ 386 /* After the matrix and forcing vector are filled, solvem is called to */ 387 /* the matrix problem. The solution vector is then read in the same manner */ 388 /* it was generated, ie, trivial points are zeroed and skipped in reading */

389 /* the solution vector. 390 /* */ 391 /******************************************************************************/ 392 xsolv() 393 int flagtriv, flagtrif; /* flagtriv and flagtrif are used by the source loop 394 and the field loop respectively to flag trivial 395 points. Values are either True (defined above as 396 1) or False: 0. 397 int ipois, ipoif; /* ipois and ipoif denote the source and field points 398 respectively. This is in contradistinction to irow 399 and icol. */ 400 int irow, icol; /* icol and irow denote positions within THE matrix. 401 If x symmetry is not being used, these variables 402 will contain the same values as ipois and ipoif. 403 If x symmetry is being used, then icol and irow will 404 gradually fall behind ipois and ipoif as trivial 405 points are encountered. */ 406 for ( flagtrif = False, ipoif = 0, irow = 0; ipoif < mnumpoi; ipoif++, irow += 407!flagtrif, flagtrif = False ) ( 408 if (flagsym e 1 &&!point.x[ipoif]) {flagtrif = True; continue;} 409 for ( flagtriv=False,ipois=O,icol=O;ipois < mnumpoi; ipois++, icol += 410!flagtriv, flagtriv=False) { 411 if (flagsym & 1 t!point.x[ipois]) {flagtriv=True; continue;} 412 if (!taui) Matrix(irow,icol)=(tau-l)*ppcon(ipoif,ipois); 413 else { 414 Cmatrix(irow,icol).real=(tau-l)*ppcon(ipoif,ipois); 415 Cmatrix(irow,icol).imag= taui*ppcon(ipoif,ipois); 416 } 417 if (irow == icol) { 418 if (!taui) Matrix(irow,icol) += 1; 419 else Cmatrix(irow,icol).real += 1; 420 } 421 422 if (!taui) fvect[irow]=point.x[ipoif]; 423 else { 424 cfvect[irow].real=point.x[ipoif]; 425 cfvect[irow].imag=0; 426 } 427 } 428 if (!taui) ( 429 pmat(); 430 pvec(; 431 e 432 solvem(irow);

433 if (!taui) { 434 pmat(); 435 pvecO; 436 ( 437 for ( flagtrif=False, ipoif=0, irow=O; ipoif < mnumpop- ipoif++, irow += 438!flagtrif, flagtrif = False ) { 439 if (flagsym e 1 U!point.x[ipoif]) flagtrif = True; 440 if (!taui) point.poten[ipoif] = (flagtrif)? 0.: fvect[irow]; 441 else { 442 point.cpoten[ipoif].real = (flagtrif)?O.:cfvect[irow].real; 443 point.cpoten[ipoif].imag = (flagtrif)?O.:cfvect[irow].imag; 444 f 445 } 446 > 447 /***************************************************************************** 448 /* 449 /* ROUTINE: ysolv */ 450 /* */ 451 /* ysolv is called by main to fill THE matrix. ysolv assumes excitation is */ 452 /* in the y direction. Points which lie on y=O are assumed to have zero 453 /* potential and are skipped in both the row and column loops. The outside */ 454 /* loop is the column loop and is associated with the field point. The */ 455 /* inside loop is associated with the source point. Images of the source */ 456 /* point are not apparent at this level and are handled in ppcon. */ 457 /* After the matrix and forcing vector are filled, solvem is called to */ 458 /* the matrix problem. The solution vector is then read in the same manner */ 459 /* it was generated, ie, trivial points are zeroed and skipped in reading */ 460 /* the solution vector. 461 /* 462 /********************************************************************* 463 ysolv() { 464 int flagtriv, flagtrif; /* flagtriv and flagtrif are used by the source loop 465 and the field loop respectively to flag trivial 466 points. Values are either True (defined above as 467 1) or False: 0. */ 468 int ipos, ipoif; /* ipois and ipoif denote the source and field points 469 respectively. This is in contradistinction to irow 470 and icol. */ 471 int irow, icol; /* icol and irow denote positions within THE matrix. 472 If y symmetry is not being used, these variables 473 will contain the same values as ipois and ipoif. 474 If y symmetry is being used, then icol and irow will 475 gradually fall behind ipois and ipoif as trivial 476 points are encountered. */

477 for ( flagtrif=False. ipoif=0. irow=0; ipoif < mnumpoi; ipoif++, irow += 478!flagtrif, flagtrif = False ) { 479 if (flagsym & 2 &!point.y[ipoif]) {flagtrif = True; continue;} 480 for ( flagtriv=False,ipois=O,icol=0;ipois < mnumpoi; ipois++, icol += 481!flagtriv, flagtriv=False) ( 482 if (flagsym 2 &&!point.y[ipois]) (flagtriv=True; continue;} 483 if (!taui) Matrix(irow,icol)=(tau-l)*ppcon(ipoifipois); 484 else ( 485 Cmatrix(irow.icol).real=(tau-l)*ppcon(ipoif.ipois); 486 Cmatrix(irow,icol).imag=taui*ppcon(ipoif,ipois); 487 488 if (irow == icol) { 489 if (!taui) Matrix(irow,icol) += 1; 490 else Cmatrix(irow,icol).real += 1; 491 } 492 } 493 if (!taui) fvect[irow]=point.y[ipoif]; 494 else { 495 cfvect[irow].real=point.y[ipoif]; 496 cfvect[irow].imag=0; 497 } 498 } 499 solvem(irow); 500 for ( flagtrif=False, ipoif=0, irow=0; ipoif < mnumpoi; ipoif++, irow += 501!flagtrif, flagtrif = False ) { 502 if (flagsym k 2 &a!point.y[ipoif]) flagtrif = True; 503 if (!taui) point.poten[ipoif] = (flagtrif)? 0.: fvect[irow]; 604 else { 505 point.cpoten[ipoif]. real=(flagtrif)?0.:cfvect[irow]. real; 506 point.cpotenipoif].imag=(flagtrif)?0.:cfvect[irow].imag; 507 } 508 509 } 510 /*****************************************************************************/ 511 /* 512 /* ROUTINE: ppcon(ipoif,ipois) */ 513 /* 514 /* ppcon is called by xsolv and ysolv to calculate the point to point 515 /* contribution. ppcon processes all triangles associated with the source */ 516 /* point. Imaging of the source patches due to symmetry is not explicitly */ 517 /* performed by pppcon but is handled by lower level routines. All */ 518 /* processing is done by calling ptcon */ 519 /* */ 520 ***************** *********** ******** i

521 double ppcon(ipoif,ipois) 522 int ipoif, ipois; /* field point index, source point index */ 523 { 524 int i; /* loop variable */ 525 double sum; /* scratch variable, holds sum of contributions */ 526 for (i=O,sum=O; i<apoint.numatri[ipois]; i++) 527 sum += ptcon(ipoif, ipois, apoint.Poii(i,ipois), 528 apoint.Poi2(i.ipois)); 529 return(sum/(4*Pi)); 530 } 531 /*******************************************************************************/ 532 /* */ 533 /* ROUTINE: ptcon(ipoif,ipois,ipoil,ipoi2) */ 534 /* */ 535 /* ptcon calculates the contribution from a source triangle to a field */ 536 /* point. The source triangle is delimited by the vertices associated with */ 537 /* ipois, ipoil, and ipoi2. Depending upon the shape of the triangle, the */ 538 /* field contribution is calculated using one of three formulas. */ 539 /* The coordinates of all points are transformed to the u-v coordinate */ 540 /* system using the functions u, v, ur, and vr. ur and vr are used to */ 541 /* calculate coordinates of patches in reflected quadrants. traninit is */ 542 /* used to initialize the translation. */ 543 /* */ 544 /******************************************************************************/ 545 double ptcon(ipoif,ipois,ipoil,ipoi2) 546 int ipoif, ipois, ipoil, ipoi2; /* indices for: field point, source point, and 547 two delimiting points of source triangle */ 548 { 549 int i; /* scratch variable */ 550 double bsum; /* stores sum of contribution from first quadrant and all reflected 551 quadrants */ 552 double sum; /* stores line integral contributions to surface integral */ 553 for (flagsyms=O,bsum=0; flagsyms<=flagsym; flagsyms += (flagsym&i)?1:2) { 554 traninit(ipoilipoi2); 555 if (ur(ipois) < 0) { 556 i=ipoil; 557 ipoil=ipoi2; 558 ipoi2=i; 559 traninit(ipoil,ipoi2); 560 } 561 sum = (fabs(ur(ipois)) * Epsilon > fabs(vr(ipois)))?0.: 562 intdvz(ur(ipois)/vr(ipois),O.,00.,vr(lipois),ipoif); 563 sum += (fabs(ur(ipois)) * Epsilon > fabs(vr(ipoi2)-vr(ipois)))?0.: 564 intdvz(-ur(ipois)/(vr(ipoi2)-vr(ipois)),ur(ipois)*vr(ipoi2)

565 /(vr(ipoi2)-vr(ipois)),vr(ipois),vr(ipoi2),ipoif); 566 sum -= intdvz(0..,0.,O.vr(ipoi2),ipoif); 567 sum /= fabs(ur(ipois)); 568 bsum += (potflag!= 3 &e potflag t flagsyms)? -sum: sum; 569 } 570 return(bsum); 571 } 572 /******************************************************************************/ 573 /* */ 574 /* Additional External Variables */ 575 /* */ 576 /* Used in routines: traninitu,v,ur,vr */ 577 /1 */ 578 /********************************************************************* 579 double xO,myO,calp,salp; /* variables used in mapping to u-v coordinate system. 580 xO, and myO are the coordinates of the origin, and 681 calp and salp are the cosine and sine of the angle 582 alpha used to rotate the x-y coordinate system. */ 583 /* myo is used instead of yO so as to not be 584 confused with the system supplied bessel 585 function which is accessed by yO. */ 586 687 /**************************************************************************** 688 /* */ 689 /* ROUTINE: traninit(ipoil.ipoi2) */ 590 /* */ 591 /* traninit initializes the translation routines u,v.ur, and vr. Initial */ 592 /* calculations are made to obtain the appropriate displacements and */ 593 /* rotations which will be needed in the translation routines. */ 594 /* */ 595 ******************************************************************************/ 596 traninit(ipoilipoi2) 697 int ipoil,ipoi2; /* indices to two auxiliary points used in defining source 598 source triangle. u-v coordinate system will be centered 699 on the point associated with ipoil, with the point 600 associated with ipoi2 lying on line u=O. / 601 { 602 double alp; /* this is the angle of rotation */ 603 switch (flagsyms) { 604 case 0: xO=point.x[ipoil]; 605 myO=point.y[ipoil]; 606 alp=atan2 (point.y [ipoi2] -myO,point.x[ipoi2]-xO); 607 calp=cos(alp); 608 salp=sin(alp);

609 break; 610 case 1: xO= -point.x[ipoil]; 611 myO=point.y[ipoil]; 612 alp=atan2(point.y[ipoi2]-myO,-point.x[ipoi2]-xO); 613 calp=cos(alp); 614 salp=sin(alp); 615 break; 616 case 2: x0=point.x[ipoil]; 617 myO= -point.y[ipoil]; 618 alp=atan2(-point.y[ipoi2 -myO.point.x[ipoi2]-xO); 619 calp=cos(alp); 620 salp=sin(alp); 621 break; 622 case 3: xO= -point.x[ipoil]; 623 myO= -point.y[ipoil]; 624 alp=atan2(-point.y[ipoi2]-myO,-point.x[ipoi2]-xO); 625 calp=cos(alp); 626 salp=sin(alp); 627 break; 628 629 return; 630 } 631 /******************************************************************************/ 632 /* */ 633 /* ROUTINES: u(i),v(i),ur(),vr(i) */ 634 */ 635 /* These routines are used to calculate coordinates of points in the u-v */ 636 /* coordinate system. i is an index to point. u and v return the */ 637 /* coordinates of the point, and ur and vr return the coordinates of an */ 638 /* image of the point, the location of the image depending on the value of */ 639 /* flagsyms. */ 640 */ 641 /****************************************************************************** 642 double u(i) 643 int i; /* index to point */ 644 { 645 return((pont.x -xO)*(-salp)+(point.y[]-myO)*c-alp)(po i; 646 } 647 double v(i) 648 int i; /* index to point */ 649 { 650 return((point.x[i]-xO)*calp+(point.y[i]-myO)*salp); 651 } 652 double ur(i)

653 int i; /* index to point */ 654 { 655 switch (flagsyms) ( 656 case 0: return((point.x[i]-xO)*(-salp)+(point.y[i]-myO)*calp); 657 case 1: return((-point.x[i]-xO)*(-salp)+(point.y[i]-myO)*calp); 658 case 2: return((point.x[i]-xO)*(-(salp)+(-point.y[i]-my)*calp); 659 case 3: return((-point.x[i]-xO)*(-salp)+(-point.y[i]-myO)*calp); 660 } 661 } 662 double vr(i) 663 int i; /* index to point */ 664 { 665 switch (flagsyms) { 666 case 0: return((point.x[i]-xO)*calp+(point.y[ i]-myO)*salp); 667 case 1: return((-point.x[i]-xO)*calp+(point.y[i]-myO)*salp); 668 case 2: return((point.x[i]-xO)*calp+(-point.y[i]-myO)*salp); 669 case 3: return((-point.x[i]-xO)*calp+(-point.y[i]-myO)*salp); 670 671 } 672 ************************************************************** **** 673 /* 674 /* ROUTINE: rdata 67 /* */ 676 /* rdata is called by the main program to read in point and triangle 677 /* definitions. Upon encountering a non- "p" or 't' command, rdata calls */ 678 /* ldata to link the data, and then returns to the main program. 679 /* */ 680 /****************************************** ** ***** ***** 681 rdata() { 682 int i; /* loop variable */ 683 int inumb; /* number of next point or triangle */ 684 do { 685 switch (com.let) { 686 case P': 687 case 'p': 688 scanf("8%d",inumb); 689 scanf("%le %le',point.x+inumbpoint.y+inumb); 690 mnumpoi=inumb+l; 691 continue; 692 case 'T': 693 case t': 694 scanf("%d",,inumb); 695 for (i=0;i<3;i++) scanf('%d',&triang.Poi(inumbi)); 696 mnumtri=inumb+l;

697 continue; 698 default: 699 ldata(); 700 return; 701 702 703 while (scanf("ls",com.str) == 1); 704 705 /**************************************************************************** 706 /* */ 707 /* ROUTINE: ldata 708 /* */ 709 /* Idata links all information concerning points and triangles. This */ 710 /* information is needed to precisely define each patch on the plate. */ 711 /** 712 /************************* ****** *** * ******* *********** 713 ldata() { 714 int i,j; /* loop variables */ 715 int isl; /* scratch variable */ 716 for (i=O;i<mnumpoi;i++) apoint.numatri[i]=0; 717 for (i=O;i<mnumtri;i++) { 718 for (triang.centroid.x[i]=triang.centroid.y[i]=,0J=0;J<3;J++) { 719 isl=triang.Poi(i. ); 720 triang.centroid.x i] += point.x[isl]/3; 721 triang.centroid.y[i] += point.y[isl]/3; 722 apoint.Tri(apoint.numatri[isl] isl) = i; 723 apoint.Poil(apoint.numatri[isl],isl) = triang.Poi(i,(J+l)%3); 724 apoint.Poi2(apoint.numatri[isl]++,isl) = triang.Poi(i,(j+2)%3); 725 } 726 727 return; 728 } 729 ****************************************************************************** 730 /* */ 731 /* ROUTINE: intdvz(a.b.vl,v2,ipoif) */ 732 / */ 733 /* intdvz performs integration in the v and z directions. The limits of */ 734 /* integration in the v direction are vl and v2. vl is nominally less than */ 735 /* v2 however this may be reversed depending on the geometry of the triangle */ 736 /* The integration in the z direction is from -t/2 to t/2, the observation */ 737 /* point lying on the z=t/2 plane. The u,v coordinates of the observation */ 738 /* point are given by the index ipoif to the structure point. a and b denote*/ 739 /* the variation of the u' variable with v. Specifically, u'=av'+b. 740 /* */

741 /****************************************************************************** 742 double intdvz(a,b,vl,v2,ipoif) 743 double a,b,vl,v2; /* parameters of integration, see above */ 744 int ipoif; /* index to point, denotes observation point */ 745 { 746 double sum; /* contain sum of integrations in the z direction, corresponding 747 to the integral evaluated for v'=vl and v'=v2 */ 748 sum=intdz(msqrt(l+a*a),(l+a*a)*v2*v2+(-2*v(ipoif)+2*a*(b-u(ipoif)))* 749 v2+v(ipoif)*v(ipoif)+(b-u(ipoif))*(b-u(ipoif)),2*(1+a*a)*v2+(-2*v(ipoif)+ 750 2*a* (b-u(ipoif)))); 751 sum -= intdz(msqrt(l+a*a),(l+a*a)*vl*vl+(-2*v(ipoif)+2*a*(b-u(ipoif)))* 752 vl+v(ipoif)*v(ipoif)+(b-u(ipoif))*(b-u(ipoif) )2*(l+a*a)*vl+(-2*v(ipoif)+ 753 2*a*(b-u(ipoif)))); 754 return(sum); 755 } 756 *****************************************************************************/ 757 /* 758 /* ROUTINE: intdz(A,B,C) 759 /* */ 760 /* intdz performs the z integration. The integration is done analytically, */ 761 /* following Gradshteyn and Ryzhik, formulae #2.267-1,#2.261,#2.266. */ 762 /* The integral which is evaluated is */ 763 /* 1/A*ln(2A(b+z**2)**.5+C) from -t to 0. */ 764 /* * 765 /*************************************************************************** 766 double intdz(A,B,C) 767 double A,B,C; /* parameters of integration */ 768 { 769 double a,b.zl,z2; /* new parameters of integration */ 770 int ianO; /* False (=0) if a is effectively zero */ 771 double sum; /* scratch variable, holds sum of contributions */ 772 double troot; /* scratch variable, holds root of R=a+b*x+c*z**2 */ 773 if (fabs(C) < Epsilon * fabs(A)) { 774 /* short cut for C = 0 */ 775 sum = (fabs(B) < Epsilon * fabs(A))? (2*t*mlog(t)-2*t): 776 (t*mlog(B+t*t)-2*t+2*msqrt (B) atan(t/msqrt(B))); 777 sum /= 2*A; 778 sum += t*mlog(2*A)/A; 779 return(sum); 780 } 781 if (fabs(B) < Epsilon * fabs(A)) { 782 /* short cut for B = 0 */ 783 z2 = 2*A*t+C; 784 zl = C;

785 sum = z2*mlog(z2)-z2; 786 sum -= zl*mlog(zl)-zl; 787 sum /= 2*A*A; 788 return(sum); 789 } 790 /* no short cuts */ 791 a=C*C-4*B*A*A; 792 ianO = fabs(a)>20*Epsilon*(fabs(C*C)+fabs(4*B*A*A)); 793 b= -2*C; 794 zl=2*A*msqrt(B)+C; 795 z2=2*A*msqrt(t*t+B) +C; 796 sum=troot=msqrt(a+b*z2+z2*z2); 797 if (ianO) sum += (a>O)? 798 -msqrt(a)*mlog((2*a+b*z2+2*msqrt(a)*troot)/z2): 799 -msqrt(-a)*((fabs(troot)<fabs(a)*Epsillon)? 800 (((2*a+b*z2)>0)?Pi/2: -Pi/2): 801 atan((2*a+b*z2) / (2*msqrt(-a)*troot))); 802 sum += b/2*mlog(2*troot+2*z2+b); 803 sum -= troot=msqrt(a+b*zl+zl*zl); 804 if (ianO) sum -= (a>O)? 805 -msqrt(a)*mlog((2*a+b*zl+2*msqrt(a)*troot)/zl): ~ 806 -msqrt(-a)*((fabs(troot)<fabs (a)*Epsilon)? 807 (((2*a+b*zl)>0)?Pi/2: -Pi/2): 808 atan ((2*a+b*zl) / (2*msqrt(-a)*troot))); 809 sum -= b/2*mlog(2*troot+2*zl+b); 810 sum /= -2*A*A; 811 sum += t/A*mlog(2*A*msqrt(B+t*t)+C); 812 return(sum); 813 } 814 /************************** ** ***** ******** **/ 815 /* */ 816 /* ROUTINE: solvem(N) */ 817 /*,/ 818 /* solvem is used to solve the matrix problem. The matrix and forcing vector*/ 819 /* are created in xsolv or ysolv. N is the dimension of the vector. solvem */ 820 /* solves the matrix problem by calling the two fortran routines dgeco and */ 821 /* dgesl which preform a decomposition and back substitution. solvem is the */ 822 /* only c routine which calls a fortran routine. */ 823 /* 824 /*****************************************************************************/ 825 solvem(N) 826 int N; /* N is the order of the matrix stored in the array matrix. N is less 827 than or equal to isize */ 828 ( i i

829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 int isize,job; /* arguments for fortran routines, isize is the size of the array matrix, job (=0) indicates type of matrix problem */ int ipvt[Mnumpoi]; /* an array used by the fortran routines to store the pivoting vector */ double rcond; /* the condition number of the matrix. */ double zCMnumpoi*2]; /* scratch vector, lengthened for complex case */ isize=Mnumpoi; job=0; if (taui) dgeco_(matrix,kisize,&N,ipvt,krcond,z); else cgeco (cmatrix,kisize,tN,ipvt,trcond,z); printf(Tcondition number is %e\n",1/rcond); if (!taui) dgesl (matrix,kisize,kNipvt,fvect,kjob); else cgesl (cmatrix,kisize,.N,ipvtcfvect.&job); /* */ /* ROUTINE: dipole /* */ /* dipole is called by the main program after the potentials over the plate */ /* have been computed. dipole examines each point to determine which points */ /* are exterior points and then uses these points to define the 'edge" of */ /* of the plate. This is needed to perform the surface integral, see Senior */ /* 1976, Low-frequency scattering by a dielectric body. The surface integral*/ /* gives the dipole moment. This is normalized by the volume of the body */ /* which is obtained by summing the area of the individual triangles and /* multiplying by the thickness t. In the process of isolating the exterior */ /* points, any point which is incorrectly linked will be identified. dipole() { double dip,vol; /* dip contains the dipole moment and vol contains the volume of the plate */ double scrat,sum; /* scrat and sum are scratch variables */ double sdip; /* sdip is a scratch variable */ double dipi.sdipi; /* imaginary part of dip and sdip */ double ycom.xcom,ycomp,xcomp,ycomn,xcomn; /* more scratch variables, used to determine sign of dipole contributions */ int plistl[Mnumpoi],plist2[Mnumpoi],plist3[Mnumpoi]; /* plisti and plist2 are used to hold the linked list of exterior points, plist3 hold the list of associated triangles. The contents of these arrays are pointers to point and, for plist3, to triang. */ int elist[Mnumpoil]; /* elist contains the original. unlinked, list of -- I

873 exterior points. */ 874 int i.J,k,l,m; /* loop variables */ 875 /* 876 fill elist 877 */ 878 for (i=O,k=O;i<mnumpoi;i++) { 879 scrat=sumang(i); 880 if (scrat>360.01) ( 881 printf(">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>); 882 printf("E R R 0 R <<<<<<<<<<<<<<<<<<<<<<<<<<<<< \n'); 883 printf(">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>*); 884 printf("point # %d is defined incorrectly:\n",i); 886 printf (">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>); 886 printf("total associated angle is %lf \n",scrat); 887 continue; 888 } 889 if (scrat<359.99) elist[k++] = i; 890 } 891 /* 892 if z excitation, skip the following stuff and goto zcase 893 */ 894 if (potflag == 3) goto zcase; 895 /* 896 create linked list 897 */ 898 for (i=O,m=O;i<k-1;i++) { 899 for (j=i+l; j<k; j++) { 900 for (l=0;l<apoint.numatri[elist(i]];l++) { 901 if (apoint.Poil(l,elist[i]) == elist[j]) { 902 plistl[m] = elist[i]; 903 plist2[m] = elist[j]; 904 plist3[m++] = apoint.Tri(l,elist[i]); 905 } 906 if (apoint.Poi2(l,elist[il) == elist[j]) { 907 plistl[m] = elist[i]; 908 plist2[m] = elist[j]; 909 plist3[m++] = apoint.Tri(l,elist[i]); 911 } 912 } 913 } 914 /* 915 calculate dipole moment, assuming potflag=l or 2. 916 */

917 for Ci0.~dip=0,di pi=O i<m; iI-+){ 918 if (!taui) sdip = ((potfiag == 1)? 919 fabs(point. yiplistil[ii] -point.y plist2[i]]): 920 fabs(point.x~plistl~i]]-point.x plist2[i]])) 921*ponpoeplslilpitptnpit[]/2 922 else*pinponpit[i)ontptnpitiJ/2 923 sdip = ((potfiag = ) 924 fabs~point.y~plistl[i]]-point.y~plist2(i)]): 925 fabs~point.x~plistl(i]]-point.xtplist2(i]])) 926 *(point.cpoten~p listi (1]] real 927 +point.cpoten Lp list2 Ci]].real)/2; 928 sdipi = ((Potfiag == 1)? 929 fabs~point.y~plistl~i]J-point.y plist2[il])): 930 fabsC(point.xlplistl~i]]-point.xlplist2Ei]J)) 931 *Cpoint.cpoten p listiI~i]].mag 932 +point.cpoten Lp list2 Ci]].ia)2 9iag/2 934 / 935 now determine correct sign 936 * 937 ycom = point.y plistl(i]]-point.y~plist2(iJ]; 938 xcom = point.x[plistl~i]]-point.x~plist2(i]J; 939 ycomp = point.Y yplistl CiJJ-triang.centroid.y~plist3[iJJ; 940 xcomp = point.xtplistl~i]J-triang.centroid.xlplist3CiJ]; 941 ycomn = -xcom; 942 xcomn = ycom; 943 scrat = Cycomp*ycomn+xcomp*xcomn) *(C(potf lag ==1)?xcomn:ycomn); 944 dip += sdip*C(scrat>0)?l:-l); 945 if Ctaui) dipi += sdipi*C(scrat>0)?1:-1); 946} 947 948 / 949 calculate volume of plate 950 * 951 for Ci=O~vol=O;i<mnumtri;i++) 952 vol '= area~point.x~triang.Poi~i,0)J,point.y~triang.Poi~i,0)], 953 point.x~triang.Poi~i 1)],point.y~triang.Poi~i.1)). 954 point.x~triang.Poi (1.2)],point.y~triang.Poi~i,2)]); 955 dipmom (tau-1)*dip/vol-taui*dipi/vol; 956 dipmomi ( Ctaui)?0.:Ctau-l)*dipi/vol+taui*dip/vol; 957 return; 958 /********** 959 end of x or y dipole computation 960 beginning of z dipole computation I -j -i _n I

961 label zcase is only branched to from 962 one spot, immediately following error 963 checking routine. 964 *********************/ 965 zcase: 966 for (i=0,vol=O,dip=O,dipi=O;i<mnumtri;i++) { 967 scrat = area(point.x[triang. Poi(i,0) ],point.y[triang.Poi(i,0)], 968 point.x[triang. Po(i. 1)i ]point.y triang.Poi(i.1) ] 969 point.x[triang.Poi(i.2) ].point.y triang.Poi(i2) ]); 970 vol += scrat; 971 if (!taui) { 972 sdip = point.poten[triang.Poi(i,0) ]; 973 sdip += point.poten[triang.Pol(i,1)]; 974 sdip += point.poten[triang. Poi (i.2) ]; 975 sdip *= scrat/3; 976 dip += sdip; 977 } 978 else ( 979 sdip = point.cpoten[triang.Po(i1,0)].real; 980 sdipi = point.cpoten(triang.Poi(i,0) ].imag; 981 sdip += point.cpoten[triang.Poi(i.)].real; 982 sdipi += point.cpoten[triang.Poi(i,1)]. imag; G 983 sdip += point.cpoten[triang.Poi (1,2)].real; 984 sdipi += point.cpoten[triang.Pol(i,2)]. imag; 985 sdip *= scrat/3; 986 sdipi *= scrat/3; 987 dip += sdip; 988 dipi += sdipi; 989 } 990 } 991 /* note that vol actually contains area */ 992 dipmom = -2*(tau-l)*dip/(vol*t)+2*taui*dipi/(vol*t); 993 dipmomi = (!taui)?0:-2*(tau-1)*dipi/(vol*t)-2*taui*dip/(vol*t); 994 return; 995 > 996 /******************************************************************************/ 997 /* */ 998 /* ROUTINE: sumang(ipoint) */ 999 /* */ 1000 /* sumang is called by dipole to produce the sum of the angles associated */ 1001 /* with the point ipoint. If the sum is 360, then the point is an interior */ 1002 /* point, if the sum is less than 360. then the point is an exterior point */ 1003 /* and defines the edge of the plate */ 1004 /* */ i i

1005 ****************************************************************************** 1006 #define Pxl (polnt.x[apoint.Poll(iipoint)]-point.x[ipoint]) /* shorthand */ 1007 #define Px2 (point.x[apoint.Poi2(i,ipoint)]point.xpoint]) /* shorthand */ 1008 #define Pyl (point.y[apoint.Poil(i,ipoint)]-point.y[ipoint]) /* shorthand */ 1009 #define Py2 (point.y[apoint.Poi2(i,ipointoint)]pont[ipoint]) /* shorthand */ 1010 double sumang(ipoint) 1011 int ipoint; /* pointer to the structure point */ 1012 { 1013 int i; /* loop variable */ 1014 double scale,prod,sum; /* scale contains the product of the lengths of the 1015 two edges, prod contains the dot product of the 1016 two edges, and sum contains a running total of the 1017 angles. */ 1018 for (i=O,sum=O;i<apoint.numatri[ipoint];i++) { 1019 scale = (Pxl*Pxl+Pyl*Pyl); 1020 scale *= (Px2*Px2+Py2*Py2); 1021 scale = sqrt(scale); 1022 prod = Pxl*Px2+Pyl*Py2; 1023 sum += acos(prod/scale); 1024 } 1025 if (!point.x[ipoint] U flagsym k 1) sum *= 2; 1026 if (!point.y[ipoint] e flagsym & 2) sum *= 2; 1027 sum *= 180/Pi; 1028 return(sum); 1029 } 1030 *************************************************************************** 1031 /* */ 1032 /* ROUTINE: gethed */ 1033 /* */ 1034 /* gethed is called by the main program to get the heading of the data. */ 1035 /* gethed fills hedstr with the remainder of the line (the entire line */ 1036 /* except for the letter H which must be in column 1. */ 1037 /* */ 1038 *****************************************************************************/ 1039 gethed() { 1040 int i; /* loop variable */ 1041 for (i=0; (hedstr[i]=getchar())!= '\n'; i++); 1042 hedstr[i] = '\0'; 1043 return; 1044 } 1045 /***************************************************************************** 1046 * */ 1047 /* ROUTINE: pexc */ 1048 /**

1049 /* pexc is called by the main program immediately after calling xsolv or */ 1050 /* ysolv. pexc calls dipole to calculate the dipole moment, and prints the */ 1051 /* result, along with information such as the heading, the material */ 1052 /* parameters, the direction of excitation and the implicit symmetry of the */ 1053 /* plate. 1054 / */ 1055 /**************************************************************************** 1056 pexc() { 1057 printf("%s \n',hedstr); 1058 printf("t = %lf, tau = %lf + i %lf\n',t,tau.taui); 1059 switch (potflag) ( 1060 case 0: 1061 printf('potential has not been computed\n'); 1062 break; 1063 case 1: 1064 printf('potential has been computed assuming x excitation\n"); 1065 break; 1066 case 2: 1067 printf("potential has been computed assuming y excitation\n"); 1068 break; 1069 case 3: 1070 printf("potential has been computed assuming z excitation\n'); c 1071 break; 1072 } 1073 switch (flagsym) { 1074 case 0: 1075 printf("plate is not symmetric\n'); 1076 break; 1077 case 1: 1078 printf('plate has mirror symmetry about x=0\n"); 1079 break; 1080 case 2: 1081 printf('plate has mirror symmetry about y=O\n"); 1082 break; 1083 case 3: 1084 printf("plate has mirror symmetry about x=O and y=O\n'); 1085 break; 1086 ) 1087 dipole(); 1088 printf("The dipole moment is %lf + i lf\n,dipmom,dipmo pmomi); 1089 } 1090 1091 pmat()O 1092 int ij; __ )0

1093 for(i=;i<lO;printf("\n"),i++) 1094 for(j=0;j<10;j++) printf("%lf."Matrix(i,j)); 1095 } 1096 pvec()( 1097 int i; 1098 for(i=0;i<10;i++) printf("%lf *,fvect[i]); 1099 printf("\n'); 1100 } 1101 double msqrt(x) 1102 double x; 1103 { 1104 if (x> -le-15) return(sqrt(x)); 1105 return(sqrt(x)); 1106 } 1107 double mlog(x) 1108 double x; 1109 { 1110 if (x> -le-16) return(log(x)); 1111 return(log(x)); 1112 } 1113 /* function definitions / - 1114 #include <math.h> 1115 double zppcon().zptcono()zplcon(),zintl(),zint2(),p().1(),alpo(); 1116 /* macro definitions */ 1117 #define Darray(Ara Il,I2,L1,L2) Ara[(Il)+(I2)*(L1)] /* Macro to generate other\ 1118 macros to mimic fortran\ 1119 arrays */ 1120 #define Mnumpoi 100 /* maximum number of points */ 1121 #define Mnumtri 100 /* maximum number of triangular patches */ 1122 #define Mnumatri 10 /* maximum number of triangles which may share a \ 1123 common point. 6 is average for internal points, 8 accounts for\ 1124 most schemes, thus 10 is a reasonable upper bound */ 1125 #define True 1 /* value of logical true */ 1126 #define False 0 /* value of logical false */ 1127 #define Tri(Jl,J2) Darray(tri,J1,J2,Mnumatri.Mnumpoi) /* set up tri as two\ 1128 dimensional array with limits (Mnumatri. Mnumpoi) */ 1129 #define Poil(Jl.J2) Darray(poil,JlJ2,Mnumatri.Mnumpoi) /* set up poll as two\ 1130 dimensional array with limits (Mnumatri, Mnumpoi) */ 1131 #define Poi2(JlJ2) Darray(poi2.Jl.J2.Mnumatri,Mnumpoi) /* set up poi2 as two\ 1132 dimensional array with limits (Mnumatri, Mnumpoi) */ 1133 #define Poi(J1,J2) Darray(poi.Jl,J2.Mnumtri,3) /* set up poi as two dimensional\ 1134 array with limits (Mnumtri,3) */ 1135 #define Matrix(Jl,J2) Darray(matrix,J1,J2,Mnumpoi.Mnumpoi) /* set up matrix as\ 1136 two dimensional array with limits (Mnumpoi.Mnumpoi) */ i i

1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 #define Cmatrix(Jl,J2) Darray(cmatrix,Jl,J2,Mnumpoi,Mnumpoi) /* complex version of Matrix, used with the complex structure */ struct complex { double real; double imag; }; /* /* /* End of Definitions / Beginning of external variable declarations /* /* extern *1 *1 *1 *1 int flagsym; /* extern int flagsyms; /, flagsym may assume values of 0,1,2, or 3. Flagsym is used by low level routines to incorporate symmetry in a manner invisible to the rest of the program. 0 - indicates no symmetry 1 - indicates object possesses mirror symmetry about x=0 2 - indicates object possesses mirror symmetry about y-0 3 - indicates object possesses mirror symmetry about x=0 and y=0 Flagsym is set at the beginning of the main program and after that is never changed */ '* flagsyms is used in conjunction with flagsym to generate mirror images of the source patches while calculating contributions to the field point. The initial object description is assumed to lie in the first quadrant, however if flagsym is 0, this is not necessary. Flagsyms is always less than or equal to flagsym, and may assume the values of 0,1,2, or 3. These values are given the following meanings: 0 - source patch is original patch (presumably first quadrant) 1 - source patch is in second quadrant 2 - source patch is in fourth quadrant 3 - source patch is in third quadrant */ x[Mnumpoi]; y [Mnumpoi]; poten [Mnumpoil]; complex cpoten[Mnumpoi]; /* complex version of poten */; /* point is a structure which contains the locations of all of the points used in defining the triangular patches. A point number 0 is permitted as C defines arrays beginning with 0. X and y are the coordinates of the points, and r) I extern struct { double double double struct }point

1181 poten is the potential of each point, as determined from the 1182 matrix problem. */ 1183 extern 1184 struct { 1185 int numatri[Mnumpoi]; 1186 int Tri(Mnumatri,Mnumpoi-1); 1187 int Poil(Mnumatri,Mnumpoi-1); 1188 int Poi2(Mnumatri,Mnumpoi-l); 1189 } apoint; /* apoint is a structure which contains lists of triangles 1190 associated with each point. numatri contains the number of 1191 triangles associated with each point; tri contains the list of 1192 triangle numbers associated with each point; and poll and poi2 1193 contain the other two vertices used in defining the triangle. 1194 poll and poi2 are indices to point. tri is an index to 1195 triang. */ 1196 extern 1197 struct { 1198 struct { 1199 double x(Mnumtri]; 1200 double y[Mnumtri]; 1201 } centroid; 1202 int Poi(Mnumtri,3-1); r 1203 double poten[Mnumtri]; 1204 struct complex cpoten[Mnumtri]; /* complex version of poten */ 1205 } triang; /* triang is used to solve the scattering problem with z 1206 excitation. centroid.x and.y contain the coordinates 1207 of the centroid of each triangle. poi contains the list of 1208 points (vertices) associated with each triangle, and is an 1209 index to point. poten is the potential of each triangle as 1210 determined from solution of the matrix problem. */ 1211 extern 1212 double Matrix(Mnumpoi,Mnumpoi-l); /* This is OTHE' matrix. Trivial points 1213 (ie, those which have zero potential from 1214 symmetry considerations) are skipped when 1215 reading or filling the matrix, which is 1216. always done sequentially. */ 1217 extern 1218 struct complex Cmatrix(Mnumpoi,Mnumpoi-1); /* complex version of Matrix */ 1219 extern 1220 double fvect[Mnumpoi]; /* fvect is the forcing vector for the matrix problem 1221 and after solution of the matrix problem contains 1222 the solution vector. */ 1223 extern 1224 struct complex cfvect[Mnumpoi]; /* complex version of fvect */ -I

1225 extern 1226 int mnumpoi; /* This is the total number of points. Points are assumed to be 1227 numbered sequentially from 0. mnumpoi <= Mnumpoi. */ 1228 extern 1229 int mnumtri; /* This is the total number of triangles. Triangles are assumed 1230 to be numbered sequentially from 0. mnumtri <= Mnumtri. */ 1231 extern 1232 double t; /* Thickness of the plate */ 1233 extern 1234 double tau; /* permittivity of the plate */ 1235 extern 1236 double taui; /* imaginary part of the permittivity. This variable is also 1237 used as a flag: if taui is zero, computations are performed 1238 assuming "tau" is purely real; if taui is non-zero, the routines 1239 appropriate for a complex tau are invoked */ 1240 extern 1241 double dipmom; /* contains the real part of the computed dipole moment */ 1242 extern 1243 double dipmomi; /* contains the imaginary part of the computed dipole moment */ 1244 extern 1245 double Epsilon; /* A very small number, used for approximate equality */ r 1246 extern 1247 double Pi; /* Pi */ 1248 /* 1249 local global definitions 1250 */ 1251 #define Pointx(I) (point.x[I]*(flagsyms t 1?-1:1)) /* reflect x coordinate */ 1252 #define Pointy(I) (point.y[I]*(flagsyms a 2?-1:1)) /* reflect y coordinate */ 1253 /* 1254 variable declarations 1255 */ 1256 double d; /* set to either t or 0, depending on which side of the 1257 surface is being integrated */ 1258 double zxO,zyO,alp,calp,salp; 1259 /* these variables are set by zrotat when the p-1 coordinate 1260 system is defined. They are used by p() and 10. */ 1261 /******************************************************************************/ 1262 /* */ 1263 /* ROUTINE: zsolv() */ 1264 /* zsolv is called by the main program to calculate the potential on a */ 1265 /* resulting from z excitation. zsolv fills Matrix using zppcon (point- */ 1266 /* point) contributions. The matrix problem is solved by solvem. */ 1267 /* */ 1268 /******************************************************************************/ 0.m

1269 zsolv() ( 1270 int ipoif,ipois; /* pointers to field point and source point, respectively */ 1271 for (ipoif=O; ipoif < mnumpoi; ipoif++) ( 1272 for(ipois=0; ipois<mnumpoi; ipois++) { 1273 if (!taui) Matrix(ipoif,ipois) = (tau-1)/t*zppcon(ipoif,ipois); 1274 else ( 1275 Cmatrix(ipoif,ipois).real=(tau-l)/t*zppcon(ipoif,ipois); 1276 Cmatrix(ipoif, ipois).imag=taui/t*zppcon(ipoif ipois); 1277 } 1278 if (ipoif == ipois) { 1279 if (!taui) Matrix(lpoif,ipois) +=.5; 1280 else Cmatrix(ipoif,ipois).real +=.5; 1281 } 1282 } 1283 if (!taui) fvect[ipoif] = -t/2; 1284 else { 1285 cfvect[ipoif].real = -t/2; 1286 cfvect[ipoif].imag = 0; 1287 } 1288 } 1289 if (!taui) { 1290 pmat(); 1291 vec(); 1292 e 1293 solvem(ipoif); 1294 if (!taui) { 1295 pmat(); 1296 pvec); 1297 > 1298 for (ipoif=O; ipoif < mnumpoi; ipoif++) 1299 if (itaui) point.poten[ipoif] = fvect[ipoif]/2; 1300 else { 1301 point.cpoten[ipoif].real = cfvect[ipoif].real/2; 1302 point.cpoten[ipoif].imag = cfvect[ipoif].imag/2; 1303 o 1304 } 1305 ****************************************************************************** 1306 /* * 1307 /* ROUTINE: zppcon(ipoif,ipois) */ 1308 /* zppcon is called by zsolv to calculate the point to point contribution */ 1309 /* to the integral. zppcon evaluates the surface integral on the top and */ 1310 /* the bottom by setting d=0,t and calling zptcon. */ 1311 /* */ 1312 /************************** ******************** **************

1313 double zppcon(ipoif,ipois) 1314 int ipoif,ipois; /* pointers to field and source points, respectively */ 1315 { 1316 int i; /* loop variable */ 1317 double sum; /* holds sum of contributions from both sides of all 1318 relevant triangles */ 1319 for (i=O,sum=O; i<apoint.numatri[ipoisl;i++) { 1320 d=O; 1321 sum += zptcon(ipolf,ipois,apoint.Poil(i,ipois),apoint.Poi2(i,ipois)); 1322 d=t; 1323 sum -= zptcon(ipoif,ipois,apoint.Poil(i,ipois),apoint.Poi2(i,ipois)); 1324 } 1325 return(sum/(4*Pi)); 1326 } 1327 /*************************************************************************** 1328 /* */ 1329 /* ROUTINE: zptcon(ipoif,ipois,ip,ipoeip ) */ 1330 /* zptcon is called by zppcon to calculate the point-triangle contribution */ 1331 /* to the surface integral. zptcon accomplishes this by calling zplcon */ 1332 /* to evaluate the line integral about the triangle (the surface integral */ 1333 /* is converted to a line integral via Wilton's method). */ 1334 /* */ 1335 /*******************************************************************************/ 1336 double zptcon(ipoif,ipois,ipoel,ipoe2) 1337 int ipoif,ipois,ipoellipoe2; /* pointers to the field point, source point, 1338 first point of edge opposing ipois in desired triangle, and 1339 second point of edge opposing ipois in desired triangle. */ 1340 { 1341 double sum; /* holds sum of integration on three edges of triangle */ 1342 for (sum=0,flagsyms=0; flagsyms <= flagsym; flagsyms += (flagsym & 1)?1:2) { 1343 sum += zplcon(ipoif,ipois,ipoelipoe 2,ipois.ipoel) 1344 sum += zplcon(ipoifipois,ipoel,ipoe2,ipoelipoe2) 1345 sum += zplcon(ipoif,ipois,ipoeli,ipoe2,ipois) 1346 } 1347 return(sum); 1348 } 1349 /****************************************************************************** 1350 /* */ 1351 /* ROUTINE: zplcon(ipoifipoisipoes,lpoel,2poe2,ipo ipo2) */ 1352 /* zplcon is called by zptcon to calculate the contribution of one line of */ 1353 /* a triangle to a point. zplcon performs this by evaluating two line */ 1354 /* integrals via the routines zinti and zint2. */ 1355 /* */ 1356 ************************************************************************* I

1357 #define Dot(Xl,Y1,X2,Y2) (X1*X2+Y1*Y2) /* define dot product */ 1358 #define Vx(I2,Il) (Pointx(I2)-Pointx(I1)) /* x component of vector */ 1359 #define Vy(I2,Il) (Pointy(I2)-Pointy(I1)) /* y component of vector */ 1360 #define Vfx(I1) (Pointx(Il)-point.xLipoif]) /* x component relative to ipoif */ 1361 #define Vfy(Ii) (Pointy(I1)-point.y[ipoif]) /* y component relative to ipoif */ 1362 #define Unit(X,Y) temp=sqrt(X*X+Y*Y); \ 1363 X / temp; \ 1364 Y /= temp; /* normalize the vector pair X.Y */ 1365 double zplcon(ipoif,ipois,ipoel,ipoe2,ipoil.ipoi2) 1366 int ipoif,ipois,ipoel,ipoe2,ipoil,ipoi2; 1367 /* ipoif is the pointer to the field point, 1368 ipois is the pointer to the source point, 1369 ipoel is the pointer to the first point associated with the 1370 edge which opposes the source point in the current triangle. 1371 ipoe2 is the pointer to the second point associated with the 1372 edge which opposes the source point in the current triangle. 1373 ipoil is the pointer to the first point associated with the 1374 line over which the integration is to be performed in the 1375 current triangle, 1376 ipoi2 is the pointer to the second point associated with the 1377 line over which the integration is to be performed in the 1378 current triangle, u 1379 */ 1380 { 1381 int ipoi3; /* ipoi3 is the third point in the triangle ipoil.ipoi2.ipoi3 */ 1382 double xO,yO; /* will contain coordinates of the point on line ipoil1383 ipoi2 which is closest to ipoif */ 1384 double a,b,c; /* contains the parameters which describe line ipoil-ipoi2 */ 1385 int iscr; /* scratch variable */ 1386 double sum; /* holds result of integrations */ 1387 double tx.ty; /* temporary variables, used to hold x and y components of 1388 a vector */ 1389 double temp; /* used only in the macro Unit, hold intermediate computation */ 1390 double uxuy; /* holds x and y components of the unit vector u */ 1391 double qx,qy; /* holds x and y components of the unit vector q */ 1392 double scalel,scale2; /* these are the factors by which the two integrations, 1393 zintl and zint2. are scaled */ 1394 double mq; /* holds the magnitude of the q vector */ 1395 double tl; /* scratch variable */ 1396 /* calculate scaling factors */ 1397 tx = Vx(ipoi2,ipoil); 1398 ty = Vy(ipoi2,ipoil); 1399 ux = ty; 1400 uy = -tx; n) rI

1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 Unit(ux,uy); ipoi3 = (!((ipoil-ipoel)*(ipoi2-ipoel))? ((ipoil-ipoe2)*(ipoi2-ipoe2)?ipoe2 ipos):ipoel); if (Dot(Vx(ipoil,ipoi3),Vy(ipoil,ipoi3),ux,uy)<0) ( ux = -ux; uy = -uy; } /* now find q, first find qhat */ tx = Vx(ipoe2,ipoel); ty = Vy(ipoe2,ipoel); qx = ty; qy = -tx; Unit(qx,qy); if (Dot(Vx(ipoel,ipois),Vy(ipoel,ipois),qx,qy)>0) ( qx = -qx; qy = -qy; /* next find mag of q, mq */ tl = -Dot(Vx(ipoel,ipois),Vy(ipoel,ipois),qx,qy); mq = 1/tl; /* now have enough to get u dot q, next get offset */ tl = Dot(Vfx(ipoel),Vfy(ipoel),qx,qy)*mq; scalel = -tl; scale2 = mq*Dot(ux,uy,qx,qy); /* find xO,yO: needed to set up 1-p coordinate system integration, cy=ax+b, describes line ipoil - ipoi2 if (fabs(Pointx(ipoill)-Pointx(ipoi2)) < Epsilon) { c=O; a=l; b = -Pointx(ipoil); else { I oa I to facilitate */ c=l; a=(Pointy(ipoil)-Pointy(ipoi2))/(Pointx(ipoil)-Pointx(ipoi2)); b=Pointy(ipoil)-a*Pointx(ipoil); if (c==O) { xO=Pointx(ipoil); yO=point.y[ipoif]; else { /* find closest point on line to ipoif */ xO = (point.x[ipoif]+a*otpoint.y[ipof]-a*b)/l+a*a); y = a*xO+b;

1445 /* generate (p.1) coordinate system with ipoif at center and xO,yO on p axis */ 1446 zrotat(point.x ipoif],point.y(ipoif],xO,yO,ipoil.ipoi2); 1447 if (l(ipoil) > l(ipoi2)) ( 1448 iscr = ipoil; 1449 ipoil=ipoi2; 1450 ipoi2=iscr; 1451 } 1452 if (p(ipoil)>Epsilon) 1453 sum = scalel*zintl(ipoif.ipois,ipoel,ipoe2,ipoilipoi2) 1454 * ((p(ipoil)>p(ipoi3))?1:-1); 1455 else sum = 0; 1456 sum += scale2*zint2(ipoifi,poisipoel,ipoe2,ipoilipoi2); 1457 if (point.x[ipoif] == Pointx(ipoil) & point.y[ipoif == Pointy(ipoil)) 1458 sum -= scalel*d*alpo(ipoi ipois,ipoel,ipoe2)/2; 1459 if (point.xCipoif] == Pointx(ipoi2) U& point.y[ipoif] == Pointy(ipoi2)) 1460 sum -= scalel*d*alpo(ipoi2,ipois,ipoel,ipoe2)/2; 1461 return(sum); 1462 } 1463 /******************** *******************************************************'*/ 1464 /* */ 1465 /* ROUTINE: zintl(ipoif.ipoisipoel,ipoe2,ipoil.ipoi2) */ 1466 /* zintl is called by zplcon to calculate the line integral resulting from */ 1467 /* the surface integral 1/R. The line integral is evaluated over the line */ 1468 /* segment ipoil - ipoi2. */ 1469 /* */ 1470 /******************************************************************************/ 1471 #define R(I) (sqrt(l(I)*l(I) p (I) *p(I)+d*d)) 1472 double zintl(ipoif,ipois,ip,ipoe2l,poie2,pol ipo2) 1473 int ipoif,ipoisipoelipoe2.ipoil,ipoi2; 1474 /* ipoif is the pointer to the field point, 1475 ipois is the pointer to the source point, 1476 ipoel is the pointer to the first point associated with the 1477 edge which opposes the source point in the current triangle, 1478 ipoe2 is the pointer to the second point associated with the 1479 edge which opposes the source point in the current triangle, 1480 ipoil is the pointer to the first point associated with the 1481 line over which the integration is to be performed in the 1482 current triangle. 1483 ipoi2 is the pointer to the second point associated with the 1484 line over which the integration is to be performed in the 1485 current triangle, 1486 */ 1487 ( 1488 double pO; /* distance from ipoif to line ipoilipoi2 */ ) 4

1489 double sum; /* holds contributions to integral */ 1490 sum=O; 1491 pO=p(ipoil); 1492 if (pO < Epsilon) return(sum); 1493 sum += log(l(ipoi2)+sqrt(pO*pO+d*d+l(ipoi2)*l(ipoi2))); 1494 sum -= log(l(ipoil)+sqrt(pO*pO+d*d+l(ipoil)*l(ipoil))); 1495 sum *= pO; 1496 if (d > Epsilon) { 1497 sum += d*atan(d*l(ipoi2)/(pO*R(ipoi2))); 1498 sum -= d*atan(d*l(ipoil)/(pO*R(ipoil))); 1499 } 1500 return(sum); 1501 } 1602 /******************************************************************************/ 1603 /* */ 1604 /* ROUTINE: zint2(ipoif,ipois.ipoel,ipoe2.ipoil,ipoi2) */ 1605 /* zint2 is called by zplcon to calculate the line integral resulting from */ 1506 /* the surface integral x/R. The line integral is evaluated over the line */ 1607 /* segment ipoil - ipoi2. */ 1608 /* */ - 1609 /******************************************************************************/ 1610 double zint2(1poif,ipols.ipoelipoe2.1poll,ipoi2) 1511 int ipoif,ipois.ipoel.ipoe2,ipoil.ipoi2; 1612 /* ipoif is the pointer to the field point. 1513 ipois is the pointer to the source point, 1514 ipoel is the pointer to the first point associated with the 1615 edge which opposes the source point in the current triangle, 1516 ipoe2 is the pointer to the second point associated with the 1617 edge which opposes the source point in the current triangle, 1518 ipoil is the pointer to the first point associated with the 1519 line over which the integration is to be performed in the 1520 current triangle, 1521 ipoi2 is the pointer to the second point associated with the 1622 line over which the integration is to be performed in the 1523 current triangle, 1524 */ 1625 < 1526 double pO; /* distance from ipoif to line ipoil-ipoi2 */ 1527 double sum; /* holds contributions to integral */ 1528 double r02; /* may equal pO**2 or pO**2+t**2, depending on d */ 1629 sum=0; 1530 p0=p(ipoil); 1531 rO2=pO*pO+d*d; 1532 if(r02<Epsilon){ i

1533 sum += l(ipoi2)*l(ipoi2); 1534 sum -= ipol(ipol )*l(iipol)(l(ipol)l(po2)>0?1:-); 1535 sum = fabs(sum); 1536 sum /= 2; 1537 } 1538 else { 1539 sum += r02*log(R(ipoi2)+1(ipoi2))+l(ipoi2)*R(ipoi2); 1540 sum -= r02*log(R(ipoil)+l(ipoil))+l(ipoil)*R(ipoil); 1541 sum /= 2; 1542 } 1543 return(sum); 1544 } 1545 ************************************************************************** 1546 /* */ 1547 /* ROUTINE: zrotat(xl,yl,x2,y2,ipoil,poi2),/ 1548 /* zrotat is called by zplcon to initialize the (p,l) coordinate system. / 1549 /* zrotat initializes several external variables which are used by p() and */ 1550 /*. */ 1551 /* */ 1552 ************************************************************************** 1553 zrotat(xl,yl,x2,y2ipoil,ipoi2) 1554 double xl,yl,x2,y2; /* xl,yl are the coordinates of the origin of the (p,1) 1555 coordinate system. x2,y2 are the coordinate of a point 1656 which should lie on the p axis */ 1557 int ipoil,ipoi2; /* pointers, used if xl=x2 and yl=y2 */ 1658 { 1559 zxO=xl; 1560 zyO=yl; 1561 if (fabs(x2-xl)+fabs(y2-yl)<Epsilon) 1562 /* line up with perpendicular to ipoil-ipoi2 */ 1563 alp=atan2(Pointx(ipoi2)-Pointx(ipoil),-(Pointy(ipoi2)-Pointy(ipoil))); 1564 else /* use standard orientation */ 1565 alp=atan2(y2-yl,x2-xl); 1566 calp=cos(alp); 1567 salp=sin(alp); 1568 return; 1569 } 1570 ***************************************************************************** 1571 /* 1572 /* ROUTINES: l(il) and p(il) */ 1573 /* l(il) and p(il) return, respectively, the reflected coordinates of the */ 1574 /* point ii. Since 1() and p() are never used with ipoif, it is safe to */ 1575 /* reflect all points. */ 1576 /* */

1577 /***************************************************************************** 1578 double l(il) 1579 int ii; 1580 { 1581 return((Pointx(il)-zx0)*(-salp)+(Pointy(il)-zyO)*calp); 1582 } 1583 double p(il) 1584 int 11; 1585 { 1586 return((Pointx(il)-zxO)*calp+(Pointy(il)-zyO)*salp); 1587 } 1588 /******************************************************************************/ 1589 /* / 1590 /* ROUTINE: alpo(i,ipoisiel.ie2) */ 1591 /* alpo returns the angle associated with point i. The angle is delimited */ 1592 /* by ipoisiel,ie2. of which one equals i. */ 1593 /* */ 1594 /******************************************************************************/ 1595 #define Spx(I) (point.x[I]-point.x[i]) /* shorthand */ 1596 #define Spy(I) (point.y[I]-pointy[i]) /* shorthand */ 1597 double alpo(i,ipois,iel,ie2) 1598 int i.ipoisiel.ie2; /* see above, ipois, iel and ie2 delimit the triangle, 1599 of which i is one of the vertices */ 1600 ( 1601 double scale,prod,ans; /* scale contains the product of the lengths, 1602 prod contains the dot product. 1603 and ans contains the angle, ie the answer */ 1604 if (iel == i) iel=ipois; 1605 if (ie2 == i) ie2=ipois; 1606 scale = Spx(iel)*Spx(iel)+Spy(iel)*Spy(iel); 1607 scale *= Spx(ie2)*Spx(ie2)+Spy(ie2)*Spy(le2); 1608 scale = sqrt(scale); 1609 prod = Spx(iel)*Spx(ie2)+Spy(iel)*Spy(ie2); 1610 ans = acos(prod/scale); 1611 return(ans); 1612 } 1613 /* function definitions */ 1614 #include <math.h> 1615 double ppcon(). ptcon().u(),v(),ur() vr(),intdvz(),intdzo()msqrt().mlog(). 1616 area(). sumang(),zppcon (); 1617 /* macro definitions */ 1618 #define Darray(Ara.I1,I2,Ll,L2) Ara[(Il)+(I2)*(L1)] /* Macro to generate other\ 1619 macros to mimic fortran\ 1620 arrays */ I __j

1621 1622 1623 1624 1625 1626 1627 1628 1629 1630 1631 1632 1633 1634 1635 1636 1637 1638 1639 1640 1641 1642 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 1656 1657 1658 1659 1660 1661 1662 1663 1664 #define #define #define #define #define #define #define #define #define #define #define Mnumpoi 100 /* maximum number of points */ Mnumtri 100 /* maximum number of triangular patches */ Mnumatri 10 /* maximum number of triangles which may share a \ common point. 6 is average for internal points, 8 accounts for\ most schemes, thus 10 is a reasonable upper bound */ True 1 /* value of logical true */ False 0 /* value of logical false */ Tri(Jl,J2) Darray(tri,J1,J2,Mnumatri,Mnumpoi) /* set up tri as two\ dimensional array with limits (Mnumatri. Mnumpoi) */ Poil(Jl,J2) Darray(poil,J1,J2,Mnumatri,Mnumpoi) /* set up poll as two\ dimensional array with limits (Mnumatri, Mnumpoi) */ Poi2(JIJ2) Darray(poi2,Jl,J2,Mnumatri,Mnumpoi) /* set up poi2 as two\ dimensional array with limits (Mnumatri, Mnumpoi) */ Poi(J1,J2) Darray(poi,J1,J2,Mnumtri,3) /* set up poi as two dimensional\ array with limits (Mnumtri,3) */ Matrix(Jl,J2) Darray(matrix,J1,J2,Mnumpoi,Mnumpoi) /* set up matrix as\ two dimensional array with limits (Mnumpoi,Mnumpoi) */ Cmatrix(Jl, J2) Darray(cmatrix,Jl,J2,Mnumpoi,Mnumpoi) /* complex version of Matrix, used with the complex structure */ struct complex { double real; double imag; /* /, /* End of Definitions / Beginning of external variable declarations /* /* extern int potflag; /* flag marks whether or not potentials have been compul 00._. *1 *1 *1 ted 0 - potential has not been computed I - potential has been computed assuming x excitation 2 - potential has been computed assuming y excitation 3 - potential has been computed assuming z excitation,,/ extern struct { double double double struct }point x[Mnumpoil]; y[Mnumpoi]; poten[Mnumpoi]; complex cpoten[Mnumpoi]; /* complex version of poten */; /* point is a structure which contains the locations of all of the points used in defining the triangular patches. A point number 0 is permitted as C defines arrays beginning with 0.

1665 X and y are the coordinates of the points, and 1666 poten is the potential of each point, as determined from the 1667 matrix problem. */ 1668 extern 1669 double t; /* Thickness of the plate */ 1670 extern 1671 double tau; /* permittivity of the plate */ 1672 extern 1673 double taui; /* imaginary part of the permittivity. This variable is also 1674 used as a flag: if taui is zero, computations are performed 1676 assuming "tau" is purely real; if taui is non-zero, the routines 1676 appropriate for a complex tau are invoked */ 1677 extern 1678 double dipmom; /* contains the real part of the computed dipole moment */ 1679 extern 1680 double dipmomi; /* contains the imaginary part of the computed dipole moment */ 1681 extern 1682 double Epsilon; /* A very small number, used for approximate equality */ 1683 extern 1684 double Pi; /* Pi */ 1685 /**************************************************************************** 1686 /* */ 1687 /* ROUTINE: eigen 1688 /* */ 1689 /* eigen is called by the main program to calculate eigenvalues. eigen is */ 1690 /* is designed to function on a simplified version of the general plate */ 1691 /* problem, so as to affect a linear predictor model of the dipole moment. */ 1692 /* The plate must have full x-y symmetry (ie, s3), and should be described */ 1693 /* via a single triangle, with point 0 falling at the origin, point 1 lying */ 1694 /* on the y axis, and point 2 lying on the x axis. The material parameters */ 1695 /* should be specified before invoking eigen() via the E command. eigen */ 1696 /* uses the thickness specified in the m command, but ignores the */ 1697 /* permittivity. eigen examines x, y, and z excitation; and will print out */ 1698 /* the eigenvalue for each case along with a dipole moment which would result*/ 1699 /1 if tau=2 and the field variation within the plate was unity, or rather, */ 1700 /* the parameter a (see write up) equals 1. */ 1701 /* */ 1702 /******************************************************************************/ 1703 eigen() { 1704 double integra; /* holds result of integration */ 1705 double leng; /* holds x or y position of test point */ 1706 double lambda; /* holds eigenvalue */ 1707 /* first do x excitation *7 1708 potflag=1;

1709 integra = ppcon(2,2)*point.x[2]; /* normalize for unit slope */ 1710 leng = point.x[2]; 1711 lambda = 1-leng/integra; 1712 point.poten[0]=point.poten[1]=0; 1713 point.poten[2]=point.x[2]; 1714 tau=2; 1715 taui=0; 1716 dipole(); 1717 printf('thickness = %lf\n',t); 1718 printf(' x excitation, eigenvalue = Xlf, I = Xlf, punit = %lf\n',lambda,integradipmom); 1719 /* now for y excitation */ 1720 potflag=2; 1721 integra = ppcon(1,1)*point.y[1]; 1722 leng = point.y[1]; 1723 lambda = 1-leng/integra; 1724 point.poten [0]=point.poten 2]=0; 1725 point.poten[l]=point.y[1]; 1726 dipole); 1727 printf(' y excitation, eigenvalue = lf, I = Xlf, punit = %lf\n,lambda,integra.dipmom); 1728 /* and now for z excitation */ 1729 potflag=3; cc 1730 integra = zppcon(0,0)+zppcon(0,1)+zppcon(0,2); 1731 lambda = 1 - t/(2*integra); 1732 point.poten[0 =point.poten [l]=point.poten[2]=t/2; 1733 dipole ); 1734 printf(" z excitation, eigenvalue = %lf, I = Xlf, punit = Xlf\n",lambda,integradipmom); 1735 return; 1736 } 1737 /* function definitions */ 1738 #include <math.h> 1739 #include <stdio.h> 1740 double area(); 1741 /* macro definitions */ 1742 #define Darray(Ara,I.,I2,LIL2) Ara[(Il)+(I2)*(L1)] /* Macro to generate other\ 1743 macros to mimic fortran\ 1744 arrays */ 1745 #define Mnumpoi 100 /* maximum number of points */ 1746 #define Mnumtri 100 /* maximum number of triangular patches */ 1747 #define Mnumatri 10 /* maximum number of triangles which may share a \ 1748 common point, 6 is average for internal points. 8 accounts for\ 1749 most schemes, thus 10 is a reasonable upper bound */ 1750 #define Mnumplo 51 /* number of lines to plot in graph, must be odd. < 60 */ 1751 #define True 1 /* value of logical true */ 1752 #define False 0 /* value of logical false */ J D

1753 1754 1755 1756 1757 1758 1759 1760 1761 1762 1763 1764 1765 1766 1767 1768 1769 1770 1771 1772 1773 1774 1775 1776 1777 1778 1779 1780 1781 1782 1783 1784 1785 1786 1787 1788 1789 1790 1791 1792 1793 1794 1795 1796 #define Para(Jl,J2) Darray(para,J1, J2,Mnumplo,Mnumplo) /* set up para as two dimensional array with limits (Mnumplo, Mnumplo) */ #define Tri(J1.J2) Darray(tri,J1,J2,Mnumatri,Mnumpoi) /* set up tri as two\ dimensional array with limits (Mnumatri, Mnumpoi) */ #define Poi(Jl,J2) Darray(poi,J1,J2,Mnumtri,3) /* set up poi as two dimensional\ array with limits (Mnumtri,3) */ #define XO point.x[triang.Poi(i,0)] /* set up shorthand notations */ #define X1 point.x[triang.Poi(i,1)] #define X2 point.x triang.Poi(i,2)) #define YO point.y triang.Poi(i,0)] #define Y1 point.y[triang.Poi(i,l)] #define Y2 point.y[triang.Poi(i,2)] /* End of Definitions / Beginning of e /* End of Definitions / Beginning of external variable declarations */ */ */ */ extern int flagsym;, extern int flagsyms; extern int plotsym; /* flagsym may assume values of 0,1,2, or 3. Flagsym is used by low level routines to incorporate symmetry in a manner invisible to the rest of the program. 0 - indicates no symmetry 1 - indicates object possesses mirror symmetry about x=O 2 - indicates object possesses mirror symmetry about y=O 3 - indicates object possesses mirror symmetry about x=0 and y=O Flagsym is set at the beginning of the main program and after that is never changed */ /* flagsyms is used in conjunction with flagsym to generate mirror images of the source patches while calculating contributions to the field point. The initial object description is assumed to lie in the first quadrant, however if flagsym is 0. this is not necessary. Flagsyms is always less than or equal to flagsym, and may assume the values of 0,1,2, or 3. These values are given the following meanings: 0 - source patch is original patch (presumably first quadrant) 1 - source patch is in second quadrant 2 - source patch is in fourth quadrant 3 - source patch is in third quadrant */ /* plotsym indicates what mirroring is desired for plotting 0 - no mirroring 1 - mirrored about x=0 2 - mirrored about y=0 -

1797 3 - mirrored about x=O and y=O 1798 Any mirroring specified must also have been previously 1799 specified in the generation of the object (see flagsym). */ 1800 extern 1801 int potflag; /* flag marks whether or not potentials have been computed 1802 0 - potential has not been computed 1803 1 - potential has been computed assuming x excitation 1804 2 - potential has been computed assuming y excitation 1805 3 - potential has been computed assuming z excitation */ 1806 float Para(MnumploMnumplo-1); /* Para contains the potentials which will 1807 be plotted */ 1808 extern 1809 struct { 1810 double x[Mnumpoi]; 1811 double y[Mnumpoi]; 1812 double poten[Mnumpoi]; 1813 }point; /* point is a structure which contains the locations of all 1814 of the points used in defining the triangular patches. A point 1815 number 0 is permitted as C defines arrays beginning with 0. 1816 X and y are the coordinates of the points, and 1817 poten is the potential of each point, as determined from the c 1818 matrix problem. / c 1819 extern 1820 struct { 1821 struct { 1822 double x[Mnumtri]; 1823 double y[Mnumtri]; 1824 } centroid; 1825 int Poi(Mnumtri,3-1); 1826 double poten[Mnumtri]; 1827 } triang; /* triang is used to solve the scattering problem with z 1828 excitation. centroid.x and.y contain the coordinates 1829 of the centroid of each triangle. poi contains the list of 1830 points (vertices) associated with each triangle, and is an 1831 index to point. poten is the potential of each triangle as 1832 determined from solution of the matrix problem. */ 1833 extern 1834 int mnumpoi; /* This is the total number of points. Points are assumed to be 1835 numbered sequentially from 0. mnumpoi <= Mnumpoi. */ 1836 extern 1837 int mnumtri; /* This is the total number of triangles. Triangles are assumed 1838 to be numbered sequentially from 0. mnumtri <= Mnumtri. */ 1839 extern 1840 double Epsilon: /* A very small number, used for approximate equality */ J -I

1841 extern 1842 double Pi; /* Pi */ 1843 /**************** **************** **** ***** *** 1844 /* 1845 /* ROUTINE: graph */ 1846 /* */ 1847 /* graph is called by the main program to graph the calculated potentials. */ 1848 /* graph accomplishes this by rerouting the standard output to the input */ 1849 /* of a call to imprint. imprint in turn produces a graph and queues it for */ 1850 /* printing. All variables which must be transferred to graph and its */ 1851 /* subordinate routines are passed from the main program via external */ 1852 /* variables. */ 1853 /* */ 1854 *********************************************************** 1855 graph() ( 1856 FILE *p_imprint; 1857 FILE *popen(); 1858 int retl, oldstdout; 1859 printf('This shouldn't go to the file\n'); 1860 fflush(stdout); 1861 oldstdout=dup(fileno(stdout));; 1862 if (oldstdout<O) perror("dup failed'); 1863 /* save stdout pointer */ 1864 p imprint=popen("imprint -n -I -Ltektronix'); 1865 iT (p imprint == NULL) perror('popen died'); 1866 retl=aup2(fileno(p imprint),fileno(stdout)); /*make stdout the same as imprint*/ 1867 if (retl<O) perrorT"dup2 failed'); 1868 graphp(); 1869 fflush(stdout); 1870 retl=close(fileno(stdout)); 1871 if (retl<O) perror("close of stdout file"); 1872 retl=dup2(oldstdout,fileno(stdout)); 1873 if (retl<O) perror("restore of stdio failed'); 1874 retl=pclose(p imprint); 1875 if (retl<O) perror("pclose failed"); 1876 retl=close(oldstdout); 1877 if (retl<O) perror("close failed'); 1878 printf("This shouldn't go to the file either.\n'); 1879 } 1880 graphy() 1881 int i,j; 1882 float nara[2500],angle,zmax,zscale; 1883 double Pi=3.14159265; 1884 Pi *= 2./5.; A) yo

1885 1886 1887 1888 1889 1890 1891 1892 1893 1894 1895 1896 1897 1898 1899 1900 1901 1902 1903 1904 1905 1906 1907 1908 1909 1910 1911 1912 1913 1914 1915 1916 1917 1918 1919 1920 1921 1922 1923 1924 1925 1926 1927 1928 zscale = 1.; for (i=O;i<50;i++) { for (j=0;j<50;j++) { nara[i+50* ]=zscale*sin(i*Pi/10)*sin(j*Pi/10)+zscale; } angle=30.; zmax=2.; dgsurf(nara,50,angle,zmax,1,0); /* */ /* ROUTINE: graphp */ /* */ /* graphp is called by graph. graphp performs the necessary calculations to */ /* fill para in a manner acceptable to dgsurf. */ /* */ graphp() { double m12,m20,mO1; /* these variables contain the slope of the lines double myl2,my20,myO1; double bl2,b20,bOl; /* double mbyl2,mby20,mby0l connecting points 12, 20. and 01, respectively. */ /* These variables are used if the line is vertical, implying that m** should be infinity. Since this is not possible, that formulation of the line is changed to my*y=m*x+b, with my=O. Except for this case of a vertical line, my=l. */ these variables contain the y offset of the lines connecting points 12, 20, and 01, respectively. */ L; /* these variables contain the evaluation of the expression mx+b-y, where x and y are the coordinates of the centroid of the triangle, and m and b are defined above by point pairs using the formulation for the line as y=mx+b. The only importance of the mby** variables is the sign of the number they contain, since if an arbitrary point (x*,y*) substituted into the expression mx+b-y produces a similar sign, then both the centroid and the point (x*,y*) are on the same side of the line. By making this comparison for all three lines, it is determined whether or not an arbitrary point lies within a given triangle. */ /* zofset contains the offset which is added to all potentials as they are stored in para. This is!-j float zofset;

1929 1930 1931 1932 1933 1934 1935 1936 1937 1938 1939 1940 1941 1942 1943 1944 1945 1946 1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972 necessary because dgsurf requires all points to be positive. zofset is defined as float so that it may be compared with elements of para to determine if they have been filled. An exact match indicates that they have not been filled */ double zmax,zmin,zscrat; float zpmax,angle; double xmin,xmax,ymin,ym int mirx,miry; /* zmax and zmin are the maximum and minimum potentials and are, of course, needed to determine what is the proper offset to be added to the potentials. zpmax is the maximum value to be plotted (usually zmax+zofset). zscrat is a scratch variable used in performing the above described calculations. */ /* parameters which are passed to dgsurf. zpmax is the maximum value to be graphed. angle is the angle at which the perspective plot is to be made. */ max; /* These variables delimit the region occupied by the plate */ /* mirx and miry are flags used to note whether a mirror image of calculated potentials is to be plotted. This may be used in conjunctions with the symmetry specification (see main program). Thus, using symmetry, only, say, 1/4 of the potentials need be computed, and using mirror plotting, the additional points may be generated. mirx and miry may assumed values of 1,0, or -1. 0 indicates no mirroring is desired. 1 indicates exact mirroring -1 indicates that the mirrored image should be inverted The mirroring may be specified independently for the x and y directions. The type of mirroring is determined by the program by considering direction of mirroring and direction of excitation.,/ (A) 0o CO I int ixmin,iymin,ixmax,iymax;/* These variables are used to delimit the region of para which is currently being filled with values from a triangular patch. */ double txmin,txmax,tymin,tymax; /* These variables are used to delimit the triangular patch which is currently being mapped to para */ double xi,yi; /* These variables specify coordinates within the triangular patch which are being mapped to para */

1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 double 11,12,10,delta; /* These are the area coordinates: L1,L2, and L3, see, for example Zienkiewicz, The Finite Element Method for an explanation. delta is the area of the triangle and 11, 12, and 10 are the area coordinates associated with points 1, 2, and 0. respectively. */ int i,j,ix,iy; /* loop variables. ix and iy are used specifically to index para */, if (potflag == 3) for (i=O;i<mnumpoi;i++) point.poten[i: mirx = (!(plotsym e 1))?0:((potflag == 1)? -1:); miry = (!(plotsym & 2))?0:((potflag == 2)? -1:1); zmax = -le20; zmin = le20; xmax = -1e20; xmin = le20; ymax = -le20; ymin = le20; for (i=0; i < mnumpoi; i++) { if (zmax < point.poten[i]) zmax = point.poten[i]; if (zmin > point.poten[i]) zmin = point.poten[i]; if (xmax < point.x[i]) xmax = point.x[i]; if (xmin > point.x[i]) xmin = point.x[i]; if (ymax < point.y[i]) ymax = point.y[i]; if (ymin > point.y[i]) ymin = point.y[i]; zscrat=(fabs(zmax) > fabs(zmin))?fabs(zmax):fabs(zmin); if (mirx == -1 II miry == -1) { zofset = zscrat; zpmax = 2*zofset; else { zofset = (zmin>0)?0:-zmin; zpmax = zofset+zmax; for (i=0; i < Mnumplo; i++) { for (j=0; j < Mnumplo; j++) { Para(i.J) = zofset; = -point.potenLIJ; (0 I clc for (i=O; i<mnumtri; i++) { for (txmin=tymin=le20,txmax=tymax= -le20j,=0; j<3; if (txmin > point.x[triang.Poi(i,)]) txmin = point.x[triang.Poi(i J)]; if (tymin > point.y[triang.Pol(i,jJ)] j++) {

2017 tymin = point.y~triang.Poi(i jDI; 2018 if (txmax < point.xltriang.Poi(ij)]) 2019 txmax = point.x[triang.Po i(iJ)); 2020 if (tynax < point.y~triang.Poi~i,j))) 2021 tymax = point~y[triang.Poi(i,j)]; 2022 2023 zscrat= X1-XO; 2024 if (fabs~zscrat) > Epsilon) { 2025 iol = CY1-YO) /zscrat; 2026 myOl = 1; 2027 2028 else ( 2029 mol = 1; 2030 myOl = 0; 2031 2032 zscrat = X2-X1; 2033 if (fabs~zscrat) > Epsilon) { 2034 m12 = CY2-Y) /zscrat; 2035 myl2 = 1; 2036 2037 else { 2038 m12 = 1; 2039 myl2 = 0; 2040 2041 zscrat = XO-X2; 2042 if (fabs(zscrat) > Epsilon) { 2043 m20 = CYO-Y2) /zscrat; 2044 my20 = 1; 2045 2046 else { 2047 m2O = 1; 2048 my2O = 0; 2049 2050 bOl = my0l*YO-mOl*XO; 2051 b12 = my12*YI-m12*X1; 2052 b20 = my2O*Y2-m2O*X2; 2053 mbyOl = TOl*triang.centroid.x[i]+bOl-myyl*triang.centroid.yi]; 2054 mbyl2 = ml2*triang.centroid.x[i)+b12-myl2*triang.centroid.y[i]; 2055 mby2o = m20*triang.cntroid.x [i]+b20-my20m*triang.centroid.y[i]; 2056 delta = area(XO,YO,X1DY1,X2,Y2); 2057 txmin * (txmin < 0)? (1+Epsilon): (I-Epsilon); 2058 tymin *= (tynin < 0)? (1+Epsilon): (1-Epsilon); 2059 txmax * (txmax > 0)? (1+Epsilon): (l-Epsilon); 2060 tymax * (tymax > 0)? CI+Epsilon): (I-Epsilon); I

2061 2062 2063 2064 2065 2066 2067 2068 2069 2070 2071 2072 2073 2074 2075 2076 2077 2078 2079 2080 2081 2082 2083 2084 2085 2086 2087 2088 2089 2090 2091 2092 2093 2094 2095 2096 2097 2098 2099 2100 2101 2102 2103 2104 ixmin = (mirx)? ceil((txmin-xmin)*(Mnumplo-1)/2/(xmax-xmin)): ceil ((txmin-xmin)*(Mnumplo-1) / (xmax-xmin)); ixmax = (mirx)? floor((txmax-xmin)*(Mnumplo-1)/2/(xmax-xmin)): floor((txmax-xmin)*(Mnumplo-1)/(xmax-xmin)); iymin = (miry)? ceil((tymin-ymin)*(Mnumplo-1)/2/(yax-ymin)): ceil((tymin-ymin)*(Mnumplo-1)/(ymax-ymin)); iymax = (miry)? floor((tymax-ymin)*(Mnumplo-1)/2/(ymax-ymin)): floor((tymax-ymin)*(Mnumplo-l)/(ymax-ymin)); for (ix=ixmin; ix <= ixmax; ix++) ( xi = (mirx)? (ix*(xmax-xmin)*2/(Mnumplo-1)+xmin): (ix*(xmax-xmin) / (Mnumplo-1)+xmin); for (iy=iymin; iy <= iymax; iy++) { yi = (miry)? (iy*(ymax-ymin)*2/(Mnumplo-l)+ymin): (iy*(ymax-ymin)/(Mnumplo-1)+ymin); if (Para(ixiy)!= zofset II mby01*(m01*xi+b01-my0O*yi) < -Epsilon II mbyl2*(ml2*xi+bl2-myl2*yi) < -Epsilon 11 mby20*(m20*xi+b20-my20*yi)< -Epsilon ) continue; 10 = area(xi,yi,Xl,Yl,X2,Y2)/delta; 11 = area(XO,YO,xi,yi,X2,Y2)/delta; 12 = area(XO.YO,X1,Yl,xi,yi)/delta; Para(ix,iy) = 10*point.poten[triang.Poi(iO)] + 1l*point.poten[triang.Poi(i.1)] + 12*point.poten[triang.Poi(i,2) ] + zofset; } }:) { for (ix=(Mnumplo-1)/2; ix >= 0; ix —) for (iy=0;iy<Mnumplo;iy++) Para(ix+(Mnumplo-1)/2,iy)=Para(ix,iy); for (ix=O;ix<(Mnumplo-1)/2;ix++) for (iy=O;iy<Mnumplo;iy++) Para(ix.iy) = (mirx == 1)? Para(Mnumplo-l-ix,iy): -Para(Mnumplo-l-ix.iy)+2*zofset; if (mirx if (mirJ r) < for (iy=(Mnumplo-1)/2; iy >= 0; iy —) for (ix=O;ix<Mnumplo;ix++) Para(ix,iy+(Mnumplo-1)/2)=Para(ix,iy); for (iy=0;iy<(Mnumplo-1)/2;iy++) for (ix=O;ix<Mnumplo;ix++) Para(ix,iy) = (miry == 1)? Para(ix.Mnumplo-l-ly):

2105 -Para(ix,Mnumplo-l-iy)+2*zofset; 2106 ) 2107 angle=(potflag == 3)?30:75; 2108 dgsurf(para,Mnumplo,angle,zpmax. 10); 2109 } 2110 ******************************************************************************/ 2111 /* */ 2112 /* ROUTINE: area */ 2113 /* */ 2114 /* area is called by graphp to calculate the area of a triangle. This is */ 2115 /* needed in the computation of the area coordinates. */ 2116 /* */ 2117 /**************************************************************************** 2118 double area(xl,yl,x2,y2,x3,y3) 2119 double xl,yl,x2,y2,x3,y3; /* These are the coordinates of the three 2120 points which delimit the triangle */ 2121 { 2122 return(fabs(x2*y3-y2*x3-xl*y3+yl*x3+xl*y2-yl*x2)/2); 2123 ) 2124 **************************************************************************** 2125 /**/ 2126 /* ROUTINES: DSGSURF and DGZPLT */ 2127 /* */ 2128 /* These routines are proprietary and source listings are not available. */ 2129 /* */ 2130 **************************************************************** 2131 C NAASA 2.1.028 DGECO FTN-A 05-02-78 THE UNIV OF MICH COMP CTR 2132 SUBROUTINE DGECO(A,LDA,N,IPVT,RCOND,Z) 2133 INTEGER LDA,N,IPVT(1) 2134 DOUBLE PRECISION A(LDA,1),Z(1) 2135 DOUBLE PRECISION RCOND 2136 C 2137 C DGECO FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION 2138 C AND ESTIMATES THE CONDITION OF THE MATRIX. 2139 C 2140 C IF RCOND IS NOT NEEDED, DGEFA IS SLIGHTLY FASTER. 2141 C TO SOLVE A*X = B, FOLLOW DGECO BY DGESL. 2142 C TO COMPUTE INVERSE(A)*C, FOLLOW DGECO BY DGESL. 2143 C TO COMPUTE DETERMINANT(A), FOLLOW DGECO BY DGEDI. 2144 C TO COMPUTE INVERSE(A), FOLLOW DGECO BY DGEDI. 2145 C 2146 C ON ENTRY 2147 C 2148 C A DOUBLE PRECISION(LDA, N)

2149 2150 2151 2152 2153 2154 2155 2156 2157 2158 2159 2160 2161 2162 2163 2164 2165 2166 2167 2168 2169 2170 2171 2172 2173 2174 2175 2176 2177 2178 2179 2180 2181 2182 2183 2184 2185 2186 2187 2188 2189 2190 2191 2192 C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C THE MATRIX TO BE FACTORED. LDA INTEGER THE LEADING DIMENSION OF THE ARRAY A. N INTEGER THE ORDER OF THE MATRIX A. ON RETURN A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT. THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE L IS A PRODUCT OF PERMUTATION AND UNIT LOWER TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. IPVT INTEGER(N) AN INTEGER VECTOR OF PIVOT INDICES. RCOND DOUBLE PRECISION AN ESTIMATE OF THE RECIPROCAL CONDITION OF A. FOR THE SYSTEM A*X = B, RELATIVE PERTURBATIONS IN A AND B OF SIZE EPSILON MAY CAUSE RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND. IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION 1.0 + RCOND.EQ. 1.0 IS TRUE, THEN A MAY BE SINGULAR TO WORKING PRECISION. IN PARTICULAR, RCOND IS ZERO IF EXACT SINGULARITY IS DETECTED OR THE ESTIMATE UNDERFLOWS. Z DOUBLE PRECISION(N) A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS AN APPROXIMATE NULL VECTOR IN THE SENSE THAT NORM(A*Z) = RCOND*NORM(A)*NORM(Z). LINPACK. THIS VERSION DATED 07/14/77. CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. SUBROUTINES AND FUNCTIONS LINPACK DGEFA BLAS DAXPY,DDOT,DSCAL,DASUM 4 -

2193 C FORTRAN DABS,DMAX1.DSIGN 2194 C 2195 C INTERNAL VARIABLES 2196 C 2197 DOUBLE PRECISION DDOT,EK,T,WK,WKM 2198 DOUBLE PRECISION ANORM, S,DASUM,SM,YNORM 2199 INTEGER INFO,J,K,KB,KP1,L 2200 C 2201 DOUBLE PRECISION DSIGN 2202 C 2203 C COMPUTE 1-NORM OF A 2204 C 2205 ANORM = O.ODO 2206 DO 10 J = 1, N 2207 ANORM = DMAX1(ANORM, DASUM(N, A (1,J),1)) 2208 10 CONTINUE 2209 C 2210 C FACTOR 2211 C 2212 CALL DGEFA(A,LDA,N,IPVT,INFO) 2213 C 2214 C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) 2215 C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND TRANS(A)*Y = E 2216 C TRANS(A) IS THE TRANSPOSE OF A. THE COMPONENTS OF E ARE 2217 C CHOSEN TO CAUSE MAXIMUM LOCAL GROWTH IN THE ELEMENTS OF W WHERE 2218 C TRANS(U)*W = E. THE VECTORS ARE FREQUENTLY RESCALED TO AVOID 2219 C OVERFLOW. 2220 C 2221 C SOLVE TRANS(U)*W = E 2222 C 2223 EK = 1.ODO 2224 DO 20 J = 1, N 2225 Z(J) = O.ODO 2226 20 CONTINUE 2227 DO 100 K = 1, N 2228 IF (Z(K).NE. 0.ODO) EK = DSIGN(EK,-Z(K)) 2229 IF (DABS(EK-Z(K)).LE. DABS(A(K,K))) GO TO 30 2230 S = DABS(A(K,K) )/DABS (EK-Z(K)) 2231 CALL DSCAL(N,S.Z,1) 2232 EK = S*EK 2233 30 CONTINUE 2234 WK = EK - Z(K) 2235 WKM = -EK - Z(K) 2236 S = DABS(WK)

2237 SM = DABS(WKM) 2238 IF (A(K,K).EQ. O.ODO) GO TO 40 2239 WK = WK/A(K,K) 2240 WKM = WKM/A(K,K) 2241 GO TO 50 2242 40 CONTINUE 2243 VK = 1.ODO 2244 WKM = 1.ODO 2245 50 CONTINUE 2246 KP1 = K + 1 2247 IF (KP1.GT. N) GO TO 90 2248 DO 60 J = KP1, N 2249 SM = SM + DABS(Z(J)+WKM*A(K,J)) 2250 Z(J) = Z(J) + WK*A(K,J) 2251 S = S + DABS(Z(J)) 2252 60 CONTINUE 2253 IF (S.GE. SM) GO TO 80 2254 T = KM - WK 2255 WK = WKM 2256 DO 70 J = KP1, N 2257 Z(J) = Z(J) + T*A(K,J) 2258 70 CONTINUE 2259 80 CONTINUE 2260 90 CONTINUE 2261 Z(K) = WK 2262 100 CONTINUE 2263 S = 1.ODO/DASUM(N,Z,1) 2264 CALL DSCAL(N,S,Z,1) 2265 C 2266 C SOLVE TRANS(L)*Y = V 2267 C 2268 DO 120 KB = 1, N 2269 K = N + 1 - KB 2270 IF (K.LT. N) Z(K) = Z(K) + DDOT(N-K,A(K+1,K),1,Z(K+1),1) 2271 IF (DABS(Z(K)).LE. 1.ODO) GO TO 110 2272 S = 1.ODO/DABS(Z(K)) 2273 CALL DSCAL(N.SZ,1) 2274 110 CONTINUE 2275 L = IPVT(K) 2276 T = Z(L) 2277 Z(L) = Z(K) 2278 Z(K) = T 2279 120 CONTINUE 2280 S = 1.ODO/DASUM(N,Z,1) 41

2281 2282 2283 2284 2285 2286 2287 2288 2289 2290 2291 2292 2293 2294 2295 2296 2297 2298 2299 2300 2301 2302 2303 2304 2305 2306 2307 2308 2309 2310 2311 2312 2313 2314 2315 2316 2317 2318 2319 2320 2321 2322 2323 2324 C C C C C C C C C CALL DSCAL(N,S,Z,1) YNORM = 1.ODO SOLVE L*V = Y DO 140 K = 1, N L = IPVT(K) T = Z(L) Z (L) = Z(K) Z(K) = T IF (K.LT. N) CALL DAXPY(N-K,T,A(K+1,K),1,Z(K+i),1) IF (DABS(Z(K)).LE. 1.ODO) GO TO 130 S = 1.ODO/DABS(Z(K)) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM 130 CONTINUE 140 CONTINUE S = 1.ODO/DASUM(N,Z,1) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM SOLVE U*Z = V DO 160 KB = 1, N K = N + 1 - KB IF (DABS(Z(K)).LE. DABS(A(K,K))) GO TO 150 S = DABS (A(K,K))/DABS(Z(K)) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM 150 CONTINUE IF (A(K,K).NE. O.ODO) Z(K) = Z(K)/A(K,K) IF (A(K,K).EQ. O.ODO) Z(K) = 1.ODO T = -Z(K) CALL DAXPY(K-1,T,A(1,K),1,Z(1),1) 160 CONTINUE MAKE ZNORM = 1.0 S = 1.ODO/DASUM(N,Z,1) CALL DSCAL(N,S,Z,1) YNORM = S*YNORM IF (ANORM.NE. O.ODO) RCOND = YNORM/ANORM IF (ANORM.EQ. O.ODO) RCOND = O.ODO RETURN I -

2325 END 2326 C NAASA 2.1.029 DGEFA FTN-A 05-02-78 THE UNIV OF MICH COMP CTR 2327 SUBROUTINE DGEFA(A,LDA,N,IPVT,INFO) 2328 INTEGER LDA,N,IPVT(l),INFO 2329 DOUBLE PRECISION A(LDA,1) 2330 C 2331 C DGEFA FACTORS A DOUBLE PRECISION MATRIX BY GAUSSIAN ELIMINATION. 2332 C 2333 C DGEFA IS USUALLY CALLED BY DGECO, BUT IT CAN BE CALLED 2334 C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. 2335 C (TIME FOR DGECO) = (1 + 9/N)*(TIME FOR DGEFA) 2336 C 2337 C ON ENTRY 2338 C 2339 C A DOUBLE PRECISION(LDA, N) 2340 C THE MATRIX TO BE FACTORED. 2341 C 2342 C LDA INTEGER 2343 C THE LEADING DIMENSION OF THE ARRAY A 2344 C 2345 C N INTEGER 2346 C THE ORDER OF THE MATRIX A 2347 C 2348 C ON RETURN 2349 C 2350 C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS 2351 C WHICH WERE USED TO OBTAIN IT. 2352 C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE 2353 C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER 2354 C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. 2355 C 2356 C IPVT INTEGER(N) 2357 C AN INTEGER VECTOR OF PIVOT INDICES. 2358 C 2359 C INFO INTEGER 2360 C = 0 NORMAL VALUE. 2361 C = K IF U(KK).EQ. 0.0. THIS IS NOT AN ERROR 2362 C CONDITION FOR THIS SUBROUTINE, BUT IT DOES 2363 C INDICATE THAT DGESL OR DGEDI WILL DIVIDE BY ZERO 2364 C IF CALLED. USE RCOND IN DGECO FOR A RELIABLE 2365 C INDICATION OF SINGULARITY. 2366 C 2367 C LINPACK. THIS VERSION DATED 07/14/77 2368 C CLEVE MOLER: UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. I -.. I

2369 C 2370 C SUBROUTINES AND FUNCTIONS 2371 C 2372 C BLAS DAXPY,DSCAL,IDAMAX 2373 C 2374 C INTERNAL VARIABLES 2375 C 2376 DOUBLE PRECISION T 2377 INTEGER IDAMAX,J,K,KPI.L,NM1 2378 C 2379 C 2380 C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING 2381 C 2382 INFO = 0 2383 NM1 = N - 1 2384 IF (NM1.LT. 1) GO TO 70 2385 DO 60 K = 1, NM1 2386 KP1 = K + 1 2387 C 2388 C FIND L = PIVOT INDEX 2389 C L 2390 L = IDAMAX(N-K+1,A(K,K),l) + K - 1 2391 IPVT(K) = L 2392 C 2393 C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED 2394 C 2395 IF (A(L,K).EQ. O.ODO) GO TO 40 2396 C 2397 C INTERCHANGE IF NECESSARY 2398 C 2399 IF (L.EQ. K) GO TO 10 2400 T = A(L,K) 2401 A(L,K) = A(K,K) 2402 A(K,K) = T 2403 10 CONTINUE 2404 C 2405 C COMPUTE MULTIPLIERS 2406 C 2407 T = -I.ODO/A(K,K) 2408 CALL DSCAL(N-K,T,A(K+,1K),1) 2409 C 2410 C ROW ELIMINATION WITH COLUMN INDEXING 2411 C 2412 DO 30 J = KP1, N

2413 2414 2415 2416 2417 2418 2419 2420 2421 2422 2423 2424 2425 2426 2427 2428 2429 2430 2431 2432 2433 2434 2435 2436 2437 2438 T = A(L,J) IF (L.EQ. K) GO TO 20 A(L,J) = A(K,J) A(K,J) = T 20 CONTINUE CALL DAXPY(N-K,T,A(K*1,K),1.A(K+1, J),1) 30 CONTINUE GO TO 50 40 CONTINUE INFO = K 50 CONTINUE 60 CONTINUE 70 CONTINUE IPVT(N) = N IF (A(N,N).EQ. O.ODO) INFO = N RETURN END C NAASA 2.1.030 DGESL FTN-A 05-02-78 THE UNIV SUBROUTINE DGESL(A,LDA,N,IPVT,B,JOB) INTEGER LDA,N,IPVT(1),JOB DOUBLE PRECISION A(LDA,1),B(1) C C DGESL SOLVES THE DOUBLE PRECISION SYSTEM C A * X = B OR TRANS(A) * X = B C USING THE FACTORS COMPUTED BY DGECO OR DGEFA. C OF MICH COMP CTR 2439 C ON ENTRY 2440 2441 2442 2443 2444 2445 2446 2447 2448 2449 2450 2451 2452 2453 2454 2455 2456 C C C C C C C C C C C C C C C C C A DOUBLE PRECISION(LDA. N) THE OUTPUT FROM DGECO OR DGEFA. LDA INTEGER THE LEADING DIMENSION OF THE ARRAY A N INTEGER THE ORDER OF THE MATRIX A IPVT INTEGER (N) THE PIVOT VECTOR FROM DGECO OR DGEFA. B DOUBLE PRECISION(N) THE RIGHT HAND SIDE VECTOR. JOB INTEGER k.

2457 C = 0 TO SOLVE A*X = B, 2458 C = NONZERO TO SOLVE TRANS(A)*X = B WHERE 2459 C TRANS(A) IS THE TRANSPOSE. 2460 C 2461 C ON RETURN 2462 C 2463 C B THE SOLUTION VECTOR X 2464 C 2465 C ERROR CONDITION 2466 C 2467 C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A 2468 C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY 2469 C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER 2470 C SETTING OF LDA. IT WILL NOT OCCUR IF THE SUBROUTINES ARE 2471 C CALLED CORRECTLY AND IF DGECO HAS SET RCOND.GT. 0.0 2472 C OR DGEFA HAS SET INFO.EQ. 0. 2473 C 2474 C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX 2475 C WITH P COLUMNS 2476 C CALL DGECO(A,LDA,N,IPVT,RCOND,Z) 2477 C IF (RCOND IS TOO SMALL) GO TO... 2478 C DO 10 J = 1, P 2479 C CALL DGESL(A,LDA,N,IPVT,C(1,J),O) 2480 C 10 CONTINUE 2481 C 2482 C LINPACK. THIS VERSION DATED 07/14/77 2483 C CLEVE MOLER. UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. 2484 C 2485 C SUBROUTINES AND FUNCTIONS 2486 C 2487 C BLAS DAXPY,DDOT 2488 C 2489 C INTERNAL VARIABLES 2490 C 2491 DOUBLE PRECISION DDOT,T 2492 INTEGER K,KB,L,NM1 2493 C 2494 NM1 = N- 1 2495 IF (JOB.NE. 0) GO TO 50 2496 C 2497 C JOB = 0, SOLVE A * X = B 2498 C FIRST SOLVE L*Y = B 2499 C 2500 IF (NM1.LT. 1) GO TO 30 -I

2501 DO 20 K = 1, NM1 2502 L = IPVT(K) 2503 T = B(L) 2504 IF (L.EQ. K) GO TO 10 2505 B(L) = B(K) 2506 B(K) = T 2507 10 CONTINUE 2508 CALL DAXPY(N-K,T,A(K+1,K),,B(K+I),1) 2509 20 CONTINUE 2510 30 CONTINUE 2511 C 2512 C NOW SOLVE U*X = Y 2513 C 2514 DO 40 KB = 1, N 2515 K = N + 1 - KB 2516 B(K) = B(K)/A(K,K) 2517 T = -B(K) 2518 CALL DAXPY(K-1,T,A(1,K),1,B(1),1) 2519 40 CONTINUE 2520 GO TO 100 2521 50 CONTINUE 2522 C 2523 C JOB = NONZERO, SOLVE TRANS(A) * X = B 2524 C FIRST SOLVE TRANS(U)*Y = B 2525 C 2526 DO 60 K = 1, N 2527 T = DDOT(K-1,A(1,K),,B(l),1) 2528 B(K) = (B(K) - T)/A(K,K) 2529 60 CONTINUE 2530 C 2531 C NOW SOLVE TRANS(L)*X = Y 2532 C 2533 IF (NM1.LT. 1) GO TO 90 2534 DO 80 KB = 1, NM1 2535 K = N - KB 2536 B(K) = B(K) + DDOT(N-K,A(K+1,K),1,B(K+1),1) 2537 L = IPVT(K) 2538 IF (L.EQ. K) GO TO 70 2539 T = B(L) 2540 B(L) = B(K) 2541 B(K) = T 2542 70 CONTINUE 2543 80 CONTINUE 2544 90 CONTINUE -j n7

2545 100 CONTINUE 2546 RETURN 2547 END 2548 C NAASA 1.1.001 DASUM FTN-A 05-02-78 THE UNIV OF MICH COMP CTR 2549 DOUBLE PRECISION FUNCTION DASUM(N,DX,INCX) 2550 C 2551 C TAKES THE SUM OF THE ABSOLUTE VALUES. 2552 C JACK DONGARRA, LINPACK, 6/17/77. 2553 C 2554 DOUBLE PRECISION DX(1),DTEMP 2555 INTEGER I,INCX,M,MP1,N,NINCX 2556 C 2557 DASUM = O.ODO 2558 DTEMP = O.ODO 2559 IF(N.LE. O)RETURN 2560 IF(INCX.EQ.1)GOTO 20 2561 C 2562 C CODE FOR INCREMENT NOT EQUAL TO 1 2563 C 2564 NINCX = N*INCX 2565 DO 10 I = 1,NINCXINCX r 2566 DTEMP = DTEMP + DABS(DX(I)) 2667 10 CONTINUE 2568 DASUM = DTEMP 2569 RETURN 2570 C 2571 C CODE FOR INCREMENT EQUAL TO 1 2572 C 2573 C 2574 C CLEAN-UP LOOP 2575 C 2576 20 M = MOD(N,6) 2577 IF( M.EQ. 0 ) GO TO 40 2578 DO 30 I = 1,M 2579 DTEMP = DTEMP + DABS(DX(I)) 2580 30 CONTINUE 2581 IF( N.LT. 6 ) GO TO 60 2582 40 MP1 = M + 1 2583 DO 50 I = MP1,N,6 2584 DTEMP = DTEMP + DABS(DX(I)) + DABS(DX(I + 1)) + DABS(DX(I + 2)) 2585 * + DABS(DX(I + 3)) + DABS(DX(I + 4)) + DABS(DX(I + 5)) 2586 50 CONTINUE 2587 60 DASUM = DTEMP 2588 RETURN I -n \3

2589 2590 2591 2592 2593 2594 2595 2596 2597 2598 2599 2600 2601 2602 2603 2604 2605 2606 2607 2608 2609 2610 2611 2612 2613 2614 2615 2616 2617 2618 2619 2620 2621 2622 2623 2624 2625 2626 2627 2628 2629 2630 2631 2632 C C C C C C C C C C C C C C C C C END NAASA 1.1.004 DAXPY FTN-A 05-02-78 SUBROUTINE DAXPY(N, DA, DX, INCX, DY, INCY) THE UNIV OF MICH COMP CTR CONSTANT TIMES A VECTOR PLUS A VECTOR. USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. JACK DONGARRA, LINPACK, 6/17/77. DOUBLE PRECISION DX (1),DY (1),DA INTEGER I,INCX,INCY,IXIY.M,MP1,N IF(N.LE.O)RETURN IF (DA.EQ. O.ODO) RETURN IF (INCX.EQ.. AND. INCY.EQ.1)GOTO 20 CODE FOR UNEQUAL INCREMENTS OR NOT EQUAL TO 1 EQUAL INCREMENTS IX = 1 IY = I IF(INCX.LT.O)IX = IF(INCY.LT.0)IY = DO 10 I = 1,N DY(IY) = DY(IY) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN (-N+I)*INCX + 1 (-N+1)*INCY + 1 + DA*DX (IX) C0 -GOA. CODE FOR BOTH INCREMENTS EQUAL TO 1 CLEAN-UP LOOP 20 M = MOD(N,4) IF( M.EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DY(I) + DA*DX(I) 30 CONTINUE IF( N.LT. 4 ) RETURN 40 MPI = M + 1 DO 50 I = MP1,N,4 DY(I) = DY(I) + DA*DX(I) DY(I + 1) = DY(I + 1) + DA*DX(I + 1)

2633 2634 2635 2636 2637 2638 2639 2640 2641 2642 2643 2644 2645 2646 2647 2648 2649 2650 2651 2652 2653 2654 2655 2656 2657 2658 2659 2660 2661 2662 2663 2664 2665 2666 2667 2668 2669 2670 2671 2672 2673 2674 2675 2676 DY(I + 2) = DY(I + 2) + DA*DX(I + 2) DY(I + 3) = DY(I + 3) + DA*DX(I + 3) 50 CONTINUE RETURN END C NAASA 1.1.003 DDOT FTN-A 05-02-78 THE UNIV OF MICH COMP CTR DOUBLE PRECISION FUNCTION DDOT(N,DX,INCX,DY,INCY) C C FORMS THE DOT PRODUCT OF A VECTOR. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 6/17/77. C DOUBLE PRECISION DX(1),DY(1),DTEMP INTEGER I,INCX,INCY,IX,IY,M,MP1,N C DDOT = O.ODO DTEMP = O.ODO IF(N.LE.O)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GOTO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C C C C C C C IX = 1 IY = I IF(INCX.LT.O)IX = IF(INCY.LT.O)IY = DO 10 I = 1,N DTEMP = DTEMP + IX = IX + INCX IY = IY + INCY 10 CONTINUE DDOT = DTEMP RETURN (-N+1)*INCX + 1 (-N+1)*INCY + 1 DX(IX)*DY(IY) CODE FOR BOTH INCREMENTS EQUAL TO 1 CLEAN-UP LOOP 20 M = MOD(N,5) IF( M.EQ. 0 ) GO TO 40 DO 30 I = 1,M DTEMP = DTEMP + DX(I)*DY(I)

2677 30 CONTINUE 2678 IF( N.LT. 5 ) GO TO 60 2679 40 MP1 = M + 1 2680 DO 50 I = MP1,N,5 2681 DTEMP = DTEMP + DX(I)*DY(I) + DX(I + 1)*DY(I + 1) + 2682 * DX(I + 2)*DY(I + 2) + DX(I + 3)*DY(I + 3) + DX(I + 4)*DY(I + 4) 2683 60 CONTINUE 2684 60 DDOT = DTEMP 2685 RETURN 2686 END 2687 C NAASA i.1.009 DSCAL FTN-A 05-02-78 THE UNIV OF MICH COMP CTR 2688 SUBROUTINE DSCAL (N,DA,DX, INCX) 2689 C 2690 C SCALES A VECTOR BY A CONSTANT. 2691 C USES UNROLLED LOOPS FOR INCREMENT EQUAL TO ONE. 2692 C JACK DONGARRA, LINPACK, 6/17/77. 2693 C 2694 DOUBLE PRECISION DA,DX(1) 2695 INTEGER I,INCX,M,MP1,N,NINCX 2696 C 2697 IF (N.LE.O)RETURN 2698 IF(INCX.EQ.l)GOTO 20 2699 C 2700 C CODE FOR INCREMENT NOT EQUAL TO 1 2701 C 2702 NINCX = N*INCX 2703 DO 10 I = 1,NINCX,INCX 2704 DX(I) = DA*DX(I) 2705 10 CONTINUE 2706 RETURN 2707 C 2708 C CODE FOR INCREMENT EQUAL TO 1 2709 C 2710 C 2711 C CLEAN-UP LOOP 2712 C 2713 20 M = MOD(N,5) 2714 IF( M.EQ. 0 ) GO TO 40 2715 DO 30 I = 1,M 2716 DX(I) = DA*DX(I) 2717 30 CONTINUE 2718 IF( N.LT. 5 ) RETURN 2719 40 MP1 = M + 1 2720 DO 50 I = MP1,N,5

2721 2722 2723 2724 2725 2726 2727 2728 2729 2730 2731 2732 2733 2734 2735 2736 2737 2738 2739 2740 2741 2742 2743 2744 2745 2746 2747 2748 2749 2750 2751 2752 2753 2754 2755 2756 2757 2758 2759 2760 2761 2762 2763 2764 C C C C C C C C C DX(I) = DA*DX(I) DX(I + 1) = DA*DX(I + 1) DX(I + 2) = DA*DX(I + 2) DX(I + 3) = DA*DX(I + 3) DX(I + 4) = DA*DX(I + 4) 50 CONTINUE RETURN END NAASA 1.1.020 IDAMAX FTN-A 05-02-78 THE UNIV OF MICH COMP CTR INTEGER FUNCTION IDAMAX(N,DX,INCX) FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. JACK DONGARRA, LINPACK, 6/17/77. DOUBLE PRECISION DX(1),DMAX INTEGER I,INCX,IX,N IDAMAX = 1 IF (N.LE. 1)RETURN IF(INCX.EQ.1)GOTO 20 CODE FOR INCREMENT NOT EQUAL TO 1 cn IX = 1 DMAX = DABS(DX(1)) IX = IX + INCX DO 10 I = 2,N IF(DABS(DX(IX)).LE.DMAX) GO TO 5 IDAMAX = I DMAX = DABS(DX(IX)) IX = IX + INCX CONTINUE RETURN 5 10 C C C CODE FOR INCREMENT EQUAL TO 1 20 DMAX = DABS(DX(1)) DO 30 I = 2,N IF(DABS(DX(I)).LE.DMAX) GO TO 30 IDAMAX = I DMAX = DABS(DX(I)) 30 CONTINUE RETURN END

2765 2766 2767 2768 2769 2770 2771 2772 2773 2774 2775 2776 2777 2778 2779 2780 2781 2782 2783 2784 2785 2786 2787 2788 2789 2790 2791 2792 2793 2794 2795 2796 2797 2798 2799 2800 2801 2802 2803 2804 2805 2806 2807 2808 C NAASA 2.1.042 CGECO FTN-A 05-02-78 THE UNIV OF MICH COMP CTR SUBROUTINE CGECO(A,LDA,N,IPVTRCONDZ) IMPLICIT REAL*8(A-H,O-Z) INTEGER LDA,N,IPVT(1) COMPLEX*16 A(LDA,1),Z(1) REAL*8 RCOND C C CGECO FACTORS A COMPLEX MATRIX BY GAUSSIAN ELIMINATION C AND ESTIMATES THE CONDITION OF THE MATRIX. C C C C c c c c c c c c c c c c c c c c c c c c c c c c c c c c C C C IF TO TO TO TO ON RCOND IS NOT NEEDED, CGEFA IS SLIGHTLY FASTER. SOLVE A*X = B, FOLLOW CGECO BY CGESL. COMPUTE INVERSE(A)*C, FOLLOW CGECO BY CGESL. COMPUTE DETERMINANT(A), FOLLOW CGECO BY CGEDI. COMPUTE INVERSE(A), FOLLOW CGECO BY CGEDI. ENTRY A COMPLEX(LDA, N) THE MATRIX TO BE FACTORED. LDA INTEGER THE LEADING DIMENSION OF THE ARRAY A N INTEGER THE ORDER OF THE MATRIX A en I ON RETURN A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT. THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE L IS A PRODUCT OF PERMUTATION AND UNIT LOWER TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. IPVT INTEGER(N) AN INTEGER VECTOR OF PIVOT INDICES. RCOND REAL AN ESTIMATE OF THE RECIPROCAL CONDITION OF A FOR THE SYSTEM A*X = B, RELATIVE PERTURBATIONS IN A AND B OF SIZE EPSILON MAY CAUSE RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION

2809 C 1.0 + RCOND.EQ. 1.0 2810 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING 2811 C PRECISION. IN PARTICULAR, RCOND IS ZERO IF 2812 C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE 2813 C UNDERFLOWS. 2814 C 2815 C Z COMPLEX(N) 2816 C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. 2817 C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS 2818 C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT 2819 C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) 2820 C 2821 C LINPACK. THIS VERSION DATED 07/14/77 2822 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. 2823 C 2824 C SUBROUTINES AND FUNCTIONS 2825 C 2826 C LINPACK CGEFA 2827 C BLAS CAXPY,CDOTC,CSSCAL,SCASUM 2828 C FORTRAN DABS,DIMAG,DMAX1,DCMPLX,DCONJG,REAL 2829 C 2830 C INTERNAL VARIABLES c 2831 C 2832 IMPLICIT REAL*8(A-H,O-Z) 2833 COMPLEX*16 CDOTC,EK,T,WK,WKM 2834 REAL*8 ANORM,S,SCASUM,SM,YNORM 2835 INTEGER INFO,J,K,KB,KP1,L 2836 C 2837 COMPLEX*16 ZDUM,ZDUM1,ZDUM2,CSIGN1 2838 REAL*8 CABS1 2839 CABS1(ZDUM) = DABS(DREAL(ZDUM)) + DABS(DIMAG(ZDUM)) 2840 CSIGN1(ZDUMI,ZDUM2) = CABS1(ZDUM1)*(ZDUM2/CABS1(ZDUM2)) 2841 C 2842 C COMPUTE 1-NORM OF A 2843 C 2844 ANORM = O.ODO 2845 DO 10 J = 1, N 2846 ANORM = DMAX1(ANORM,SCASUM(N,A(1,J),1)) 2847 10 CONTINUE 2848 C 2849 C FACTOR 2850 C 2851 CALL CGEFA(A,LDA,N,IPVT,INFO) 2852 C 'I 3n )o

2853 C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))). 2854 C ESTIMATE = NORM (Z)/NORM (Y) WHERE A*Z = Y AND CTRANS(A)*Y = E. 2855 C CTRANS(A) IS THE CONJUGATE TRANSPOSE OF A. 2856 C THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL 2857 C GROWTH IN THE ELEMENTS OF W WHERE CTRANS(U)*W = E. 2858 C THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW. 2859 C 2860 C SOLVE CTRANS(U)*W = E 2861 C 2862 EK = DCMPLX (1.ODO O. ODO) 2863 DO 20 J = 1, N 2864 Z(J) = DCMPLX(O. ODO O.ODO) 2865 20 CONTINUE 2866 DO 100 K = 1, N 2867 IF (CABS1 (Z (K)).NE. O.ODO) EK = CSIGN1(EK,-Z(K)) 2868 IF (CABS1(EK-Z(K)).LE. CABS1 (A(K,K))) GO TO 30 2869 S = CABS1(A(K,K))/CABS1(EK-Z(K)) 2870 CALL CSSCAL(N.S.Z.1) 2871 EK = DCMPLX(S.O.ODO)*EK 2872 30 CONTINUE 2873 WK = EK - Z(K) 2874 WKM = -EK - Z(K) L 2875 S = CABS1(WK) 2876 SM = CABS1 (WKM) 2877 IF (CABS1(A(K,K)).EQ. O.ODO) GO TO 40 2878 WK = WK/DCONJG(A(K,K)) 2879 WKM = WKM/DCONJG(A(K K)) 2880 GO TO 50 2881 40 CONTINUE 2882 WK = DCMPLX(1. ODO, O. ODO) 2883 WKM = DCMPLX(1.ODO, O.ODO) 2884 50 CONTINUE 2885 KP1 = K + 1 2886 IF (KP1.GT. N) GO TO 90 2887 DO 60 J = KP1, N 2888 SM = SM + CABS (Z(J)+WKM*DCONJG(A(KJ))) 2889 Z(J) = Z(J) + WK*DCONJG(A(K,J)) 2890 S = S + CABS1(Z(J)) 2891 60 CONTINUE 2892 IF (S.GE. SM) GO TO 80 2893 T = WKM- WK 2894 WK = WKM 2895 DO 70 J = KP1, N 2896 Z(J) = Z(J) + T*DCONJG(A(K.J) ) 1 )

2897 70 CONTINUE 2898 80 CONTINUE 2899 90 CONTINUE 2900 Z(K) = WK 2901 100 CONTINUE 2902 S = 1.0DO/SCASUM(N, Z,1) 2903 CALL CSSCAL(N,S,Z,1) 2904 C 2905 C SOLVE CTRANS(L)*Y = V 2906 C 2907 DO 120 KB = 1, N 2908 K = N + 1 - KB 2909 IF (K.LT. N) Z(K) = Z(K) + CDOTC(N-K,A(K+1,K),1,Z(K+1),1) 2910 IF (CABS1(Z (K)).LE. 1.ODO) GO TO 110 2911 S = 1.ODO/CABSI(Z(K)) 2912 CALL CSSCAL(N,S,Z,1) 2913 110 CONTINUE 2914 L = IPVT(K) 2915 T = Z(L) 2916 Z(L) = Z(K) 2917 Z(K) = T 2918 120 CONTINUE 2919 S = 1.0DO/SCASUM(N,Z,1) 2920 CALL CSSCAL(NS,Z,1) 2921 C 2922 YNORM = 1.ODO 2923 C 2924 C SOLVE L*V = Y 2925 C 2926 DO 140 K = 1, N 2927 L = IPVT(K) 2928 T = Z(L) 2929 Z(L) = Z(K) 2930 Z(K) = T 2931 IF (K.LT. N) CALL CAXPY(N-K,T,A(K+1,K)1,,Z(K+1), 1) 2932 IF (CABS1(Z(K)).LE. 1.ODO) GO TO 130 2933 S = 1.ODO/CABS1(Z(K)) 2934 CALL CSSCAL(N,S,Z,1) 2935 YNORM = S*YNORM 2936 130 CONTINUE 2937 140 CONTINUE 2938 S = 1.ODO/SCASUM(N,Z,1) 2939 CALL CSSCAL(N,S,Z,1) 2940 YNORM = S*YNORM

2941 2942 2943 2944 2945 2946 2947 2948 2949 2950 2951 2952 2953 2954 2955 2956 2957 2958 2959 2960 2961 2962 2963 2964 2965 2966 2967 2968 2969 2970 2971 2972 2973 2974 2975 2976 2977 2978 2979 2980 2981 2982 2983 2984 C C C SOLVE U*Z = V DO 160 KB = 1, N K = N + 1 - KB IF (CABS (Z(K)).LE. CABS1(A(K,K))) GO TO 150 S = CABS1(A(K,K))/CABS1(Z(K)) CALL CSSCAL(N,S,Z,1) YNORM = S*YNORM 150 CONTINUE IF (CABSI(A(K,K)).NE. O.ODO) Z(K) = Z(K)/A(K, IF (CABS1(A(K,K)).EQ. O.ODO) Z(K) = DCMPLX(1 T = -Z(K) CALL CAXPY(K-1,T,A(1,K),1,Z(1),1) 160 CONTINUE C MAKE ZNORM = 1.0 S = 1.ODO/SCASUM(N,Z,1) CALL CSSCAL(N,S,Z,1) YNORM = S*YNORM C IF (ANORM.NE. O.ODO) RCOND = YNORM/ANORM IF (ANORM.EQ. O.ODO) RCOND = O.ODO RETURN END C NAASA 2.1.044 CGESL FTN-A 05-02-78 THE UNIV [ SUBROUTINE CGESL(A,LDA,N,IPVT,B,JOB) IMPLICIT REAL*8(A-H,O-Z) INTEGER LDA,NIPVT(1),JOB COMPLEX*16 A(LDA,1),B(1),K).ODO,O.ODO) c_-..,._ o'~i OF MICH COMP CTR C C C C C C C C C C C C C C C CGESL SOLVES THE COMPLEX SYSTEM A * X = B OR CTRANS(A) * X = B USING THE FACTORS COMPUTED BY CGECO OR CGEFA. ON ENTRY A COMPLEX(LDA, N) THE OUTPUT FROM CGECO OR CGEFA. LDA INTEGER THE LEADING DIMENSION OF THE ARRAY A. N INTEGER THE ORDER OF THE MATRIX A.

2985 C 2986 C IPVT INTEGER (N) 2987 C THE PIVOT VECTOR FROM CGECO OR CGEFA. 2988 C 2989 C B COMPLEX (N) 2990 C THE RIGHT HAND SIDE VECTOR. 2991 C 2992 C JOB INTEGER 2993 C = 0 TO SOLVE A*X = B. 2994 C = NONZERO TO SOLVE CTRANS(A)*X = B WHERE 2995 C CTRANS(A) IS THE CONJUGATE TRANSPOSE. 2996 C 2997 C ON RETURN 2998 C 2999 C B THE SOLUTION VECTOR X 3000 C 3001 C ERROR CONDITION 3002 C 3003 C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A 3004 C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY 3005 C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER 3006 C SETTING OF LDA. IT WILL NOT OCCUR IF THE SUBROUTINES ARE 3007 C CALLED CORRECTLY AND IF CGECO HAS SET RCOND.GT. 0.0 3008 C OR CGEFA HAS SET INFO.EQ. 0. 3009 C 3010 C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX 3011 C WITH P COLUMNS 3012 C CALL CGECO(A,LDA,N,IPVT,RCOND,Z) 3013 C IF (RCOND IS TOO SMALL) GO TO... 3014 C DO 10 J = 1, P 3015 C CALL CGESL(A,LDA,N,IPVT,C(1,J),O) 3016 C 10 CONTINUE 3017 C 3018 C LINPACK. THIS VERSION DATED 07/14/77 3019 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. 3020 C 3021 C SUBROUTINES AND FUNCTIONS 3022 C 3023 C BLAS CAXPY,CDOTC 3024 C FORTRAN DCONJG 3025 C 3026 C INTERNAL VARIABLES 3027 C 3028 COMPLEX 16 CDOTC,T

3029 INTEGER K,KBL,NM1 3030 C 3031 NM1 = N- 1 3032 IF (JOB.NE. O) GO TO 60 3033 C 3034 C JOB = 0, SOLVE A * X = B 3035 C FIRST SOLVE L*Y = B 3036 C 3037 IF (NM1.LT. 1) GO TO 30 3038 DO 20 K = 1, NM1 3039 L = IPVT(K) 3040 T = B(L) 3041 IF (L.EQ. K) GO TO 10 3042 B(L) = B(K) 3043 BK) = T 3044 10 CONTINUE 3045 CALL CAXPY(N-K,T,A(K+1,K)1,,B(K+1),1) 3046 20 CONTINUE 3047 30 CONTINUE 3048 C 3049 C NOW SOLVE U*X = Y 3050 C 3051 DO 40 KB = 1, N 3052 K = N + 1 - KB 3053 B(K) = B(K)/A(K,K) 3054 T = -B(K) 3055 CALL CAXPY(K-1,TA(1,K) 1,B(1),1) 3056 40 CONTINUE 3057 GO TO 100 3058 50 CONTINUE 3059 C 3060 C JOB = NONZERO, SOLVE CTRANS(A) * X = B 3061 C FIRST SOLVE CTRANS(U)*Y = B 3062 C 3063 DO 60 K = 1, N 3064 T = CDOTC(K-1,A(1,K),1,B(l),1) 3065 B(K) = (B(K) - T)/DCONJG(A(K,K)) 3066 60 CONTINUE 3067 C 3068 C NOW SOLVE CTRANS(L)*X = Y 3069 C 3070 IF (NM1.LT. 1) GO TO 90 3071 DO 80 KB = 1, NM1 3072 K N- KB

3073 3074 3075 3076 3077 3078 3079 3080 3081 3082 3083 3084 3085 3086 3087 3088 3089 3090 3091 3092 3093 3094 3095 3096 3097 3098 3099 3100 3101 3102 3103 3104 3105 3106 3107 3108 3109 3110 3111 3112 3113 3114 3115 3116 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) 8(L) = B(K) B (K) = T 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN END C NAASA 2.1.043 CGEFA FTN-A 05-02-78 THE UNIV OF MICH COMP CTR SUBROUTINE CGEFA(A,LDA,N,IPVT,INFO) IMPLICIT REAL*8(A-H, O-Z) INTEGER LDAN,IPVT(1),INFO COMPLEX*16 A(LDA,1) C C C C C C C C C C C C C C C C C C C C C C C C C C C CGEFA FACTORS A COMPLEX MATRIX BY GAUSSIAN ELIMINATION. CGEFA IS USUALLY CALLED BY CGECO, BUT IT CAN BE CALLED DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. (TIME FOR CGECO) = (1 + 9/N)*(TIME FOR CGEFA). ON ENTRY A COMPLEX(LDA, N) THE MATRIX TO BE FACTORED. l 4Oh I LDA N INTEGER THE LEADING DIMENSION OF THE ARRAY A. INTEGER THE ORDER OF THE MATRIX A. ON RETURN A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT. THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE L IS A PRODUCT OF PERMUTATION AND UNIT LOWER TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. IPVT INTEGER(N)

3117 C AN INTEGER VECTOR OF PIVOT INDICES. 3118 C 3119 C INFO INTEGER 3120 C = 0 NORMAL VALUE. 3121 C = K IF U(K,K).EQ. 0.0. THIS IS NOT AN ERROR 3122 C CONDITION FOR THIS SUBROUTINE, BUT IT DOES 3123 C INDICATE THAT CGESL OR CGEDI WILL DIVIDE BY ZERO 3124 C IF CALLED. USE RCOND IN CGECO FOR A RELIABLE 3125 C INDICATION OF SINGULARITY. 3126 C 3127 C LINPACK. THIS VERSION DATED 07/14/77 3128 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. 3129 C 3130 C SUBROUTINES AND FUNCTIONS 3131 C 3132 C BLAS CAXPY,CSCAL,ICAMAX 3133 C FORTRAN DABS,DIMAG,DCMPLX,DREAL 3134 C 3135 C INTERNAL VARIABLES 3136 C 3137 COMPLEX*16 T 3138 INTEGER ICAMAX,J,K,KP1,L,NMI c 3139 C 3140 COMPLEX*16 ZDUM 3141 REAL*8 CABS1 3142 CABS (ZDUM) = DABS(DREAL(ZDUM)) + DABS(DIMAG(ZDUM)) 3143 C 3144 C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING 3145 C 3146 INFO = 0 3147 NM1 = N - 1 3148 IF (NM1.LT. 1) GO TO 70 3149 DO 60 K = 1, NM1 3150 KP1 = K + 1 3151 C 3152 C FIND L = PIVOT INDEX 3153 C 3154 L = ICAMAX(N-K+1,A(K,K),1) + K - 1 3155 IPVT(K) = L 3156 C 3157 C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED 3158 C 3159 IF (CABS1(A(L,K)).EQ. O.ODO) GO TO 40 3160 C nI I

3161 3162 3163 3164 3165 3166 3167 3168 3169 3170 3171 3172 3173 3174 3176 3176 3177 3178 3179 3180 3181 3182 3183 3184 3185 3186 3187 3188 3189 3190 3191 3192 3193 3194 3195 3196 3197 3198 3199 3200 3201 3202 3203 3204 C C C C C C C C C C C C C C INTERCHANGE IF NECESSARY IF (L.EQ. K) GO TO 10 T = A(L,K) A (LK) = A (KK) A(K,K) = T 10 CONTINUE COMPUTE MULTIPLIERS T = -DCMPLX(1.ODO,O.ODO)/A(KK) CALL CSCAL(N-K,T,A(K+1,K),1) ROW ELIMINATION WITH COLUMN INDEXING 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),1,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.ODO) INFO = N RETURN END NAASA 1.1.014 CAXPY FTN-A 05-02-78 THE UNIV OF MICH COMP CTR SUBROUTINE CAXPY(N,CA,CXINCX,CY,INCY) CONSTANT TIMES A VECTOR PLUS A VECTOR. JACK DONGARRA, LINPACK, 6/17/77. IMPLICIT REAL*8(A-H, O-Z) COMPLEX*16 CX(1),CY(1),CA INTEGER I,INCXINCY.IX,IY,N IF (N.LE.O)RETURN

3205 IF (DABS(DREAL(CA)) + DABS(DIMAG(CA)).EQ. O.ODO ) RETURN 3206 IF(INCX.EQ.1.AND. INCY.EQ.1)GOTO 20 3207 C 3208 C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS 3209 C NOT EQUAL TO 1 3210 C 3211 IX = 1 3212 IY = 1 3213 IF(INCX.LT.O)IX = (-N+1)*INCX + 1 3214 IF(INCY.LT.O)IY = (-N+1)*INCY + 1 3215 DO 10 I = 1,N 3216 CY(IY) = CY(IY) + CA*CX(IX) 3217 IX = IX + INCX 3218 IY = IY + INCY 3219 10 CONTINUE 3220 RETURN 3221 C 3222 C CODE FOR BOTH INCREMENTS EQUAL TO 1 3223 C 3224 20 DO 30 I = 1.N 3225 CY(I) = CY(I) + CA*CX(I) 3226 30 CONTINUE o 3227 RETURN 3228 END 3229 C NAASA 1.1.012 CDOTC FTN-A 05-02-78 THE UNIV OF MICH COMP CTR 3230 COMPLEX*16 FUNCTION CDOTC(N,CX,INCXCY,INCY) 3231 C 3232 C FORMS THE DOT PRODUCT OF TWO VECTORS, CONJUGATING THE FIRST 3233 C VECTOR. 3234 C JACK DONGARRA, LINPACK, 6/17/77. 3235 C 3236 IMPLICIT REAL*8(A-H, -Z) 3237 COMPLEX*16 CX(1),CY(1),CTEMP 3238 INTEGER I,INCX,INCY,IX,IY,N 3239 C 3240 CTEMP = DCMPLX(O.ODO,O.ODO) 3241 CDOTC = DCMPLX(O.ODO,O.ODO) 3242 IF (N.LE.0)RETURN 3243 IF(INCX.EQ.1.AND.INCY.EQ.1)GOTO 20 3244 C 3245 C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS 3246 C NOT EQUAL TO 1 3247 C 3248 IX = 1

3249 3250 3251 3252 3253 3254 3255 3256 3257 3258 3259 3260 3261 3262 3263 3264 3265 3266 3267 3268 3269 3270 3271 3272 3273 3274 3275 3276 3277 3278 3279 3280 3281 3282 3283 3284 3285 3286 3287 3288 3289 3290 3291 3292 IY = 1 IF(INCX.LT.O)IX = IF(INCY.LT.O)IY = DO 10 I = 1,N CTEMP = CTEMP + IX = IX + INCX IY = IY + INCY 10 CONTINUE CDOTC = CTEMP RETURN (-N+1)*INCX + 1 (-N+1)*INCY + 1 DCONJG(CX(IX)) CY(IY) C C C CODE FOR BOTH INCREMENTS EQUAL TO 1 20 DO 30 I = 1,N CTEMP = CTEMP + DCONJG(CX(I))*CY(I) 30 CONTINUE CDOTC = CTEMP RETURN END C NAASA 1.1.018 CSSCAL FTN-A 05-02-78 SUBROUTINE CSSCAL(N,SA,CX,INCX) THE UNIV OF MICH COMP CTR C C C C C C C C C C C SCALES A COMPLEX VECTOR BY A REAL CONSTANT. JACK DONGARRA, LINPACK, 6/17/77. IMPLICIT REAL*8 (A-H, O-Z) COMPLEX*16 CX(1) REAL*8 SA INTEGER I,INCX,N,NINCX IF(N.LE.O)RETURN IF(INCX.EQ.1)GOTO 20 CODE FOR INCREMENT NOT EQUAL TO 1 NINCX = N*INCX DO 10 I = 1,NINCX,INCX CX(I) = DCMPLX(SA*DREAL(CX(I)),SA*DIMAG(CX(I))) 10 CONTINUE RETURN CODE FOR INCREMENT EQUAL TO 1 oo CO 20 DO 30 I = 1,N

3293 3294 3295 3296 3297 3298 3299 3300 3301 3302 3303 3304 3305 3306 3307 3308 3309 3310 3311 3312 3313 3314 3315 3316 3317 3318 3319 3320 3321 3322 3323 3324 3325 3326 3327 3328 3329 3330 3331 3332 3333 3334 3335 3336 CX(I) = DCMPLX(SA*DREAL(CX(I)), SA*DIMAG(CX(I))) 30 CONTINUE RETURN END C NAASA 1.1.010 SCASUM FTN-A 05-02-78 THE UNIV OF MICH COMP CTR REAL*8 FUNCTION SCASUM(N,CX,INCX) C C C C C C C C C C C C C C C C C TAKES THE SUM OF THE ABSOLUTE VALUES OF A COMPLEX VECTOR AND RETURNS A SINGLE PRECISION RESULT. JACK DONGARRA, LINPACK, 6/17/77. IMPLICIT REAL*8 (A-H, O-Z) COMPLEX*16 CX(1) REAL*8 STEMP INTEGER I. INCX,N,NINCX SCASUM = O.ODO STEMP = O.ODO IF (N. LE.0) RETURN IF(INCX.EQ.1)GOTO 20 CODE FOR INCREMENT NOT EQUAL TO 1 NINCX = N*INCX DO 10 I = 1,NINCX,INCX STEMP = STEMP + DABS(DREAL(CX(I))) + DABS(DIMAG(CX(I))) 10 CONTINUE SCASUM = STEMP RETURN CODE FOR INCREMENT EQUAL TO 1 20 DO 30 I = 1,N STEMP = STEMP + DABS(DREAL(CX(I))) + DABS(DIMAG(CX(I))) 30 CONTINUE SCASUM = STEMP RETURN END NAASA 1.1.019 CSCAL FTN-A 05-02-78 THE UNIV OF MICH COMP CTR SUBROUTINE CSCAL(N,CA,CX, INCX) SCALES A VECTOR BY A CONSTANT. JACK DONGARRA, LINPACK, 6/17/77.

3337 3338 3339 3340 3341 3342 3343 3344 3345 3346 3347 3348 3349 3350 3351 3352 3353 3354 3355 3356 3357 3358 3359 3360 3361 3362 3363 3364 3365 3366 3367 3368 3369 3370 3371 3372 3373 3374 3375 3376 3377 3378 3379 3380 IMPLICIT REAL*8 (A-H,O-Z) COMPLEX*16 CA,CX(1) INTEGER I,INCX,N. NINCX C IF(N.LE.O)RETURN IF(INCX.EQ.1)GOTO 20 C C 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 C CODE FOR INCREMENT EQUAL TO 1 C 20 DO 30 I = 1,N CX(I) = CA*CX(I) 30 CONTINUE RETURN END C NAASA 1.1.021 ICAMAX FTN-A 05-02-78 THE UNIV OF MICH COMP CTR INTEGER FUNCTION ICAMAX(N,CX,INCX) C C FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. C JACK DONGARRA, LINPACK, 6/17/77. C IMPLICIT REAL*8 (A-H,O-Z) COMPLEX*16 CX(1) REAL*8 SMAX INTEGER I,INCX,IX,N COMPLEX*16 ZDUM REAL*8 CABS1 CABS1(ZDUM) = DABS(DREAL(ZDUM)) + DABS(DIMAG(ZDUM)) C ICAMAX = 1 IF(N.LE.1)RETURN IF(INCX.EQ.1)GOTO 20 C C CODE FOR INCREMENT NOT EQUAL TO 1 C IX = 1 SMAX = CABS1(CX(1)) I I

3381 IX = IX + INCX 3382 DO 10 I = 2,N 3383 IF(CABS1(CX(IX)).LE.SMAX) GO TO 5 3384 ICAMAX = I 3385 SMAX = CABS1(CX(IX)) 3386 5 IX = IX + INCX 3387 10 CONTINUE 3388 RETURN 3389 C 3390 C CODE FOR INCREMENT EQUAL TO 1 3391 C 3392 20 SMAX = CABS1(CX(1)) 3393 DO 30 I = 2,N 3394 IF(CABS1(CX(I)).LE.SMAX) GO TO 30 3395 ICAMAX = I 3396 SMAX = CABS1(CX(I)) 3397 30 CONTINUE 3398 RETURN 3399 END 3400 CCCCCCCCCCCCCCCCCCCCCCCCC 3401 C DREAL DOESN'T SEEM TO WORK, SO THIS FUNCTION IS A SUBSTITUTE 3402 REAL*8 FUNCTION DREAL(X) 3403 COMPLEX*16 X,X2 3404 REAL*8 XA(2) 3405 EQUIVALENCE (X2,XA(1)) 3406 X2=X 3407 DREAL=XA(1) 3408 RETURN 3409 END

BIBLIOGRAPHY Bergman, D. J. (1978), "The Dielectric Constant of a Composite Material - A Problem in Classical Physics," Phys. Rep.,Vol. 43, 377-407. Bohren, C. F. and L. J. Battan (1980), "Radar Backscattering by Inhomogeneous Precipitation Particles," J. Atmos. Sci., Vol. 37, 1821-1827. Gradshteyn, I. S. and I. M. Ryzhik (1980), "Table of Integrals, Series, and Products," Academic Press, New York. Granqvist, C. G. and 0. Hunderi (1978), "Optical Properties of Ag-SiO2 Cermet Films: A Comparison of Effective-Medium Theories," Phys. Rev. B, Vol. 18, 2897-2906. Harrington, R. F. and J. R. Mautz (1975), "An Impedance Sheet Approximation for Thin Dielectric Shells," IEEE Trans. Antennas Propag., Vol. AP-23, 531-534. Harrington, R. F. (1982), "Field Computation by Moment Methods," Krieger, Malabar, FL. Herrick, D. F. (1976), "Analytical Evaluation of Kernels," Radiation Laboratory Memo No. 013714-512-M, University of Michigan. Inspektorov, E. M. (1982), "Using Integral Equations of the Second Kind for Analysis of Diffraction of Thin Shields," Izv. Vyssh, Uchebn. Zaved. Radiofiz., Vol. 25, 1099-1101. -172 -

-173 - Jackson, J. D. (1975), "Classical Electrodynamics," Wiley, New York. Keller, J. B., R. E. Kleinman, and T.B.A. Senior (1972), "Dipole Moments in Rayleigh Scattering," 3. Inst. Math. Its Appl., Vol. 9, 14-22. Kleinman, R. E. (1965), "Low Frequency Solution of Three-Dimensional Scattering Problems," Radiation Laboratory Report #7133-4-T, University of Michigan. Kock, W. E. (1948), "Metallic Delay Lenses," Bell Syst. Tech. J., Vol. 27, 58-83. Maxwell Garnett, J. C. (1904), "Colours in Metal Glasses and in Metallic Films," Philos. Trans. R. Soc. London, Vol. 203, 385-420. Mosotti, 0. F. (1850), Mem. Soc. Ital., Vol. 14, 49. Polder, D. and J. H.Van Santen (1946), "The Effective Permeability of Mixtures of Solids," Physica, Vol. 12, 257-271. Senior, T.B.A. (1975), "Low Frequency Scattering Data for Dielectric Bodies," Radiation Laboratory Memo #013714-505-M, University of Michigan. Senior, T.B.A. (1976), "Low Frequency Scattering by a Dielectric Body,"' Radio Sci., Vol. 11, 477-482. Senior, T.B.A. (1982), "Low-Frequency Scattering by a Perfectly Conducting Body," Radio Sci., Vol. 17, 741-746. Senior, T.B.A. and H. Weil (1982), "On the Validity of Modeling Rayleigh Scatterers by Spheroids," Appl. Phys. B., Vol. 29, 117-124.

-174 - Senior, T.B.A. and T. M. Willis (1982), "Rayleigh Scattering by Dielectric Bodies," IEEE Trans. Antennas and Propag., Vol. AP-30, 1271. Senior, T.B.A. (1983), "Low Frequency Scattering by a Metallic Plate," Electromagnetics, Vol. 3, 131-144. Senior, T.B.A. and D.A. Ksienski (1984), "Determination of a Vector Potential," Radio Sci., Vol. 19, 603-607. Senior, T.B.A. and M. Naor (1984), "Low Frequency Scattering by a Resistive Plate," IEEE Trans. Antennas Propag. Vol. AP-23, 272-275. Stevenson, A. F. (1953), "Solution of Electromagnetic Scattering Problems as Power Series in the Ratio (Dimension of Scatterer)/Wavelength," J. Appl. Phys., Vol. 24, 1134-1142. Stevenson, A.F. (1954), "Note on the Existence and Determination of a Vector Potential," Quart. Appl. Math., Vol. 12, 194-197. Von Hippel, A. R. (1954), "Dielectrics and Waves," Wiley, New York. Weil, H. (1984), private communications. Willis, T. M. (1982), "Low Frequency Scattering by a Thin Dielectric Plate," Radiation Laboratory Memo #01955-502-M, University of Michigan. Wilton, D. R., S. M. Rao, A. W. Glisson, D. H. Schaubert, 0. M. Al-Bundak, and C. M. Butler (1984), "Potential Integrals for Uniform and Linear Source Distributions on Polygonal and Polyhedral Domains," IEEE Trans. Antennas Propag., Vol. AP-32, 276-281.

-175 -Zienkiewicz, 0. C. (1982), "The Finite Element Method," McGrawHill, London.