388919-1-F SCATTERING BY RESISTIVE STRIPS AND PLATES by Department of Thomas B.A. Senior Radiation Laboratory Electrical Engineering and Computer Science The University of Michigan Ann Arbor, Michigan 48109 Final Report on P.O. No. 46-4385 Prepared for Northrop Corporation Aircraft Division Mail Zone 3020/92 3901 West Broadway Hawthorne, CA 90250 July 1985 388919-1-F = RL-2554

CHAPTER I. INTRODUCTION This is the final report of the work carried out for the Northrop Corporation under purchase order 36-4385 during the period ending 31 May 1985. The purpose of the study was to examine the effect of non-zero resistivity on the backscattering cross sections of strips and plates with particular reference to angles of incidence which are close to grazing. In the case of a strip with the incident electric vector parallel to the surface attention is confined to edge-on incidence, and the analysis is presented in Chapter II. The backscattering cross section is determined for a variety of uniform and quadratically-tapered resistiviti;es, and one of the significant conclusions is that tapering is not always better. When the incident magnetic vector is parallel to the surface, the scattering for edge-on incidence is zero regardless of the resistivity, and the incidence which is most important is that corresponding to the traveling wave lobe. The nature of this lobe is examined in Chapter III, and it is shown that even a small resistivity is sufficient to eliminate the traveling wave as such. Unfortunately, this merely bares a far side lobe of the specular flash. A strip is the two-dimensional analogue of a finite plate and in spite of the useful information that can be derived from a study of the simpler two-dimensional structure, there are many circumstances under which the model is irrelevant to the realistic problems of a finite plate. Substantial effort has been devoted to the development -1 -

-2 -of an efficient and effective code to compute the scattering from a finite resistive plate of arbitrary shape and resistivity. The work that has been accomplished is described in Chapter IV, and though there is still much to be done, the program in its present form does compute the current induced in the plate. The development of the program was carried out by Mr. J. W. Burns and Professor D. A. Ksienski (a present and former student, respectively, of the author), and it is a pleasure to acknowledge their assistance.

CHAPTER II. E POLARIZATION One of the more difficult problems in cross section reduction is to reduce the scattering from the edge of an electrically thin structure such as an aircraft wing or tail fin. This is particularly true at angles close to grazing incidence with the electric vector parallel to the surface, and the significant scattering that can occur in this case is dominated by the edges. We shall consider here the problem of a resistive strip or ribbon illuminated by an E-polarized plane wave at edge-on incidence (qo = r) and examine the magnitude of the backscattered field compared with that of the corresponding perfectly conducting strip. 2.1 Uniform Resistivity For a uniform resistive strip with resistivity R ohms per square the backscattered far field amplitude P(rr,r) can be expressed as a sum of front and rear edge contributions as (Senior, 1979a): P(T,r) = pf + pr (2.1) In terms of P the backscattering cross section per unit length is a= 2x P(T,1T)l2 (2.2) If the strip width w is more than about x/2 where x = 2'/k is the wavelength, the front edge return p is identical to that for a half plane of the same resistivity and -3 -

-4 - = - 6 {ZJ(O,n)}2 = - {K(O,n)}2 (2.3) - 16 where n = 2R/Z, J(x,n) is the current on a resistive half plane at a distance x from the edge for the same incident plane wave, and K(O,n) is a "split" function which appears in the Wiener-Hopf solution for a half plane. Z is the intrinsic impedance of free space and a time variation e it has been assumed. K(O,n) is real for real n and its values have been tabulated (Senior, 1979a). For complex n 212 exp {-n/7 [1 - ln(n/2)]} 1|n ~ 1 K(O,n) A- (2.4) -/2 exp {-1/(Trn)} In ~ 1 and K(O,n) can also be expressed in terms of the function i (z) introduced by Maliuzhinets (1958) as K(O,r) = 4 - i(x/2) } (2.5) ( + /1 + n)2 T/2) where cos x = 1/n. In a recent article (Volakis and Senior, 1985), two simple expressions for ~ (z) are derived which, when used in conjunction with known identities, approximate the function to a high degree of accuracy throughout the entire complex z plane. For the rear edge scattering and with the same restriction on w, r = iZJ(w, p iy{ZJ(w,n) }2 (2.6)

-5 - where (Senior, 1979b) y = {4K(O,n)}-2 (2.7) so that pr _ 1 ZJ(w,n) (2.8) r. -(2.8) p 4 ZJ(On) The above formulas are valid for complex n as well as real, and a procedure for computing the exact analytical expression for J(x,n) is described in Senior (1981). In the special case of perfect conductivity, ZJ(x,0) = 2 2 ei(x+/4) (2.9) rrkx If n is real and non-zero, J(x,n) is asymptotic to J(x,0) as kx -+ o, but if n >> 1, kx has to be very large indeed before (2.9) constitutes a good approximation to J(x,n). The current on a resistive half plane never exceeds its value at the corresponding point of a perfectly conducting half plane, and since its magnitude is a monotonically decreasing function of x and n, it follows from (2.8) that prl < - (2.10) for all w and n. For a perfectly conducting strip the front and rear edge contributions have particularly simple expressions. Since K(0,0) = /2, Eq. (2.3) gives

-6 - f i P = - (2.11) and from (2.6), (2.7) and (2.9), r e2ikw P = (2.12) This is the smallest rear edge contribution that can be achieved with any uniform resistive strip whose width w is such that kw > n, and the behavior of Ipr as a function of n is illustrated in Fig. 2.1. To make the edge-on backscattering cross section of a strip as small as possible over a wide band of frequencies, it is necessary to minimize the front and rear edge contributions individually. The smallest value of IPfl is achieved by choosing n as large as possible, and if, for example, n = 4, the resulting front edge contribution is almost 20 dB below that for a perfectly conducting strip. For the rear edge the preceding results suggest that the minimum return is obtained with a perfectly conducting strip. This is certainly true for all except the very narrowest strips, and Fig. 2.2 shows |pr as a function of w/x for five different values of (real) n. Having n t 0 inevitably increases IPr and though, for fixed w/x, the return ultimately decreases as n increases, it remains above that for n = 0 for all reasonable values of kw and n. The obvious solution is to taper the resistivity from a maximum at the front to zero at the rear, and provided this is done smoothly, it is natural to expect a rear edge contribution comparable to that for a perfectly conducting strip with no new source of scattering created.

-7 - -i W0 - 5X wax wu 3X w=6X I-2 103 I0 0 5 X 10 Fig. 2.1: The magnitudes of the front ( ---) and rear ( ) edge contribution of uniform resistive strips as functions of n. (This is an expanded version of Fig. 2 of Senior and Liepa, 1984.)

-8 - 10' i "-12 I0nX I"prl _.17 a' \ \ 0 2 W/X 4 6 Fig. 2.2: The magnitudes of the rear edge contribution ( —) for five uniform resistive strips, compared with the contribution for perfectly conducting strips ( ---) computed using (2.12).

-9 - 2.2 Quadratically-Tapered Resistivity Based on prior experience with resistive tapers, attention was confined to the quadratic form 2 -R(x) = Z (1 - x) (2.13) where x is measured from the front edge. Thus, n specifies the largest value of the resistivity, occurring at the front of the strip. Using program REST-E which solves the integral equation for an E-polarized plane wave incident on a resistive strip, the total induced current and the backscattered far field were computed as a function of w/x for a sequence of real n. From the computed data for P(7,7), the real and imaginary parts of P and Pr were extracted (Ksienski, 1985a), and the resulting values of Pf were almost identical to those for the front edge contribution of a uniform resistive strip having the same n. The magnitudes of the rear edge contributions are plotted in Fig. 2.3, and we observe that they exceed the contribution for a perfectly conducting strip but show a similar dependence on w/X. An empirical expression for Pr is r 2i kw f 5_____ P = -+ 5 /8n (2.14) 47Tkw L (kw )2)3! and since the first term in parentheses in the rear edge contribution for a perfectly conducting strip, the second term is the additive effect of the resistive taper. Under all circumstances, however, the rear edge return exceeds that for perfect conductivity, and for a narrow strip may even exceed

-10 - IP, 0 2 "/A 4 6 Fig. 2.3: The magnitudes of the rear edge contribution for tapered strips with five different values of n compared with the contribution ( ---) for perfectly conducting strips. ^' ^^^^\ \ ^^^^ \^ ^^ ^? \^^ - \\i. contribution ( ---) for perfectly conducting strips. 10 I 4 I v

-11 - the return for the corresponding uniform resistive strip. Comparison of Figs. 2.1 and 2.3 shows that if n = 10 the tapering is only effective if w > 2.4 x, and the analogous conditions for n = 6 and n = 4 are w > 1.4 x and w > 0.9 x respectively. Thus, as the strip width and/or frequency is reduced, we ultimately reach a point at which tapering becomes counter productive. Equations (2.14) and (2.3) are sufficient to provide the edge-on backscattered field of a resistive strip with a quadratic taper, and confirm that at high frequencies for which w >> X the most effective way to reduce the scattering is to increase the resistivity at the front edge to as large a value as possible and to taper the resistivity smoothly to zero at the rear. Although we have considered only the particular taper (2.13), it seems probable that the results for other smooth monotonic tapers will be similar.

CHAPTER III: H POLARIZATION Traveling waves are a major contributor to the backscattering cross section of a long thin body when the magnetic vector is perpendicular to the plane of incidence, and our work with H-polarization was concerned with studying the effect of traveling waves on the backscattering cross section of a strip or ribbon. 3.1 Traveling Wave Considerations* Traveling waves are one of the most important contributors to the scattering from long slender bodies, and in the case of structures such as the fuselage of an aircraft it is possible to produce a reasonably complete description of the backscattering by taking into account the traveling waves and, where appropriate, the specular contributions and the side lobes thereof. When the incident magnetic vector is perpendicular to the plane formed by the direction of incidence and the axis of the body, the fan-shaped pattern characteristic of a traveling wave is a key feature of the overall scattering pattern, and the first lobe is often the major contributor to the scattering at angles close to end-on incidence. The scattering pattern of a thin plate or strips also displays a similar lobe structure at angles close to edge-on incidence when the magetic vector is parallel to the surface, and it is customary to * A shortened version of this section has been published by Senior and Yang (1984). -1 2 -

-13 - attribute at least the first lobe to a traveling wave. In effect, we are thinking of each narrow slice of the plate as a filament or wire supporting a traveling wave and the angle at which the lobe appears is consistent with this argument. Mathematically, however, there is no justification for this reasoning, and it is therefore of interest to compare the high frequency expressions for the backscattered fields of wires and strips. For a wire of length z viewed at an angle e to end-on, an approximate form of the expression obtained by Peters (1958) for the backscattering cross section attributable to traveling waves is a = - S[2 where 2 SI = C^ {|!sin e sinc [ (1 - cos 0)]. (3 1) sinc. sin X sinc X = X, and r is the effective voltage reflection coefficient for the traveling wave. In deriving (3.1) we have assumed that the phase velocity is that of light and that ks >> 1, implying C = 1n(2kz) - (1 - y) where y = 0.5772.... is Euler's constant. The angles e = e n' n = 1,2,3,..., at which the maxima of the traveling wave lobes occur are given by the solutions of the transcendental equation

-14 - tan n = (1 + cos e )+n (3.2) with on 2 (1 - cos n) For large z/X 6 -49.4 (degrees), (3.3) 1 6 98.1 (degrees), (3.4) and comparison with numerical solutions of (3.2) shows that e is 1 accurate to within one percent for t/x > 1.8 and e is accurate to 2 three percent for k/x > 3.6. The results are also in good agreement with measured data. Chang and Liepa (1967) measured the backscattering patterns of a series of 81 wires having lengths varying from 0.3 to 5.42 x and radius 6.27 x 10T3 x, and the values of e and e obtained 1 2 from these patterns are shown in Fig. 3.1, along with the curves corresponding to (3.3) and (3.4). We have also used the measured amplitude of the first lobe to determine the effective reflection coefficient r in (3.1). The amplitudes oscillate in a quite regular manner with a maximum to minimum ratio of about 6 dB and maxima occurring at z/x = 0.4 + 0.5n, n = 0,1,2,... (approx.). It is clear that this is attributable to a resonance behavior of the wire, a feature which is not reproduced by (3.1), but if we incorporate this into the effective reflection coefficient, the resulting values of r appropriate to the

80' K x % x, KK NX I.. -XX x % KI % x N. LK - 60 x deg.) (deg.) 40 20' 0 x 2 -- ~ X \% X Xx K% X- Mk-*A X " 8 - I In I - - - 0 I I 2 3 4 Fig. 3.1: Measured angles (xxx) of the first and second traveling wave lobes for a thin wire. The theoretical curves were computed using (3.3) and (3.4).

-16 - first lobe are shown in Fig. 3.2. The average is about 0.65 and varies little with z/x, and this is consistent with the value generally assumed (Ruck et al., 1970) for a pointed body. We now consider an infinitesimally thin, perfectly conducting strip of width w occupying the region 0 < x < w, -o < T < ~ of the plane y = 0 of a Cartesian coordinate system x,y,z and illuminated by a plane wave having i = z eik(x cos po + y sin o) At large distances the scattered magnetic field can be written as s = e i(ip-w/4) p k e where p,( are cylindrical polar coordinates with x = p cos (, y = p sin (. In terms of P the scattering cross section per unit length in the z direction is = 2x and in the particular case of backscattering, (o =. For kw >> 1 a uniform second order GTD expression for the backscattered far field amplitude is (Senior, 1979b) p( - i 1 + co s s 1 F (/Zkw sin )} 4 cos 2

I'0' I,0 a r 0-5 -0 It II It I I I I II I' I I 1 ' '1I 1 XX i I I x I x X I I 14 *11.X / i1 t, I I 11 1 \ I A IX x I I I I I I I / I / \ I f *. y I x I I I 1 I I it II I \ I I I I I / I I k X I / ~-x X/! N iii ii -- i i _ i i i i ~ i ii 0 I I I 2 I 3 I 4 Fig. 3.2: Effective reflection coefficient r deduced from the thin wire data. curve is only to guide the eye. The dashed

-18 - i 1 - cos - e-2ikw cos ( 1 2 -irr/4 F F cos (3)5) 4 cos 2 2 + O([kw]-) valid for 0< < < r. F(T) is the Fresnel integral 00 F(T) f ei2 du T and (3.5) is identical to the result obtained by asymptotic expansion of the expression obtained by Khaskind and Vainshteyn (1964). A simple program to compute P(c, 0o) has been developed for use on an IBM PC computer. It is designated P-RIB-H and is described in Appendix A. When ( = O or f, IPI = 0 as expected, and the backscattering pattern of a strip is illustrated by the plot of IPI versus ( for w/x = 4 in Fig. 3.3. Since the pattern is symmetrical about the broadside aspect, it is sufficient to consider only 0 < ( < w/2. The lobe centered on ( = 24 degrees is the one usually attributed to a traveling wave, and the locations of this and the adjacent peaks are shown in Fig. 3.4. From the variation as a function of w/x it is evident that all peaks other than the first are most logically associated with the side lobes of the specular return. The first peak is primarily determined by the first term in (3.5), and its magnitude is (2kw (1 2kw) / 1 - C(u) - S(u) + i[C(u) - S(u)] where u = kw(l - cos ()

12 8 I 4 0 90 60 30 1f (degrees) 0 Fig. 3.3: Backscattering pattern of a strip of width w = 4x computed using (3.5).

801 x x K 60 (d eg.) 40* 20 - K x x x K K K K K K x x x x K Kx K K K K K K K "I 'Y-% K- X x -xK x -% x xx - -. —Z. K rN) C0 0 m 0 9 -I 0 I I 2 3 4 W/X Fig. 3.4: Angles at which the maxima in the pattern occur. The theoretical curve was computed using (3.6).

-21 - and C(u) and S(u) are the real integrals defined and tabulated in Abramowitz and Stegun (1964). For kw ~>> 1 the magnitude is approximately {C(u)}2 + {S(u)}2, and its maximum is 0.901 occurring at u = u = 2.29. 1 The resulting value of p is 4 =48.9 (degrees), (3.6) i V w and the corresponding curve is included in Fig. 3.4. To the next order in l/(kw) the peak value is {C(u )}2 + {S(u )}2 + 4 {C(u ) + S(u )} = 0.901 ( + 0.136 ), 1 +u01361 (1 (3.7) and this is plotted in Fig. 3.5 along with the peak values obtained from computations of |P[. The oscillations are attributable to the effect of the second term in (3.5) which was neglected in our analysis, and (3.7) is an excellent approximation to the mean. In spite of the different mathematical formulas which describe the effects of traveling waves on wires and strips, the close agreement between the angles o and q at which the first (dominant) lobe occurs is remarkable. For all practical purposes, it is sufficient to use (3.3) to locate the lobe, thereby lending support to the idea of treating a planar structure as an assemblage of elementary wires, and to using the term "traveling wave" in the case of a strip.

-22 - 1-25 - x K K K x X x KK -- - X x. K K K K KK K -- K K K- ~ IPI 0.75 - 0-5 L 0 - v I I 2 2 3;5 4 4 W/ Fig. 3.5: Peak values (xxx) of the traveling wave lobe for a strip compared with the value ( ) given by (3.7).

-23 - 2.2 Effect of Non-Zero Resistivity As we have seen, a traveling wave is a significant contributor to the backscattering cross section of a perfectly conducting strip at angles close to edge-on. It may be necessary to reduce the magnitude of the resulting peak, and the prevailing wisdom is that a relatively small amount of loss is adequate for this purpose. To determine the reduction as a function of the resistivity, we have used program REST-H or its specialized version STRIP-H (see Memorandum 2500-394-M) to compute the backscattered fields of uniform resistive strips of width w = 1.0(0.1)4.0 X and have examined the magnitudes and locations of the peaks in the backscattered patterns. For perfectly conducting strips the angle ) (measured in degrees from edge-on) at which the maxima are found to occur are shown in Fig. 3.6(a), and the results are almost identical to those obtained using the second order GTD expression (3.5) for the field. We observe that there are two distinct types of maxima; the ones associated with the main lobe of a traveling wave whose location is given approximately by the formula (3.6), and those which are more logically attributable to the side lobes of the specular flash. The latter correspond to the maxima of the pattern factor (sin X)/X with X = kw cos q, and are given by = arc cos (2n + 1) -w (3.8) with n = 1,2,3,.... As seen from Fig. 3.6(a), they lie on a sequence of trajectories which are distinct from the oscillatory curve describing the traveling wave peak. For large w/x the traveling wave dominates the far side lobes, which then show up only as an oscillation in the

801 (a) (b) *. 0. 9 0... * 0 0. 0 *, 6 0 * 4, 4, 0 0 0, 0 0 4 0 0 * 0 0 9 0, 0 0 0 4, 0 0 0 0 9, 0 0 0. 0 0 1011, 40, 0 I-,0 8O1.. 0 0.... 0,. 0 0 0, 0 0. 0... I 0, 4, 0 0 b,. a I 0, 0 0 0 0 0 0 0 0.. 0 16 " *. 0 0 0 0 6 * 0 0, * 0 0, " * e 0 0 0 0, 0 0, 4, 0 0 0 & 0 0 0 - 0 (C) (d) 40. 0 4, 0 0 0 e * 0 0, 0 0 0 0 0 0 0. 0 I I 0 I 1 2 W/,\ -1F I -9 2 W/ x 3 4 Fig. 3.6: Angles ~ (in degrees) at which peaks in the backscattering pattern occur for (a) R = 0, (b R/Z = 0.05, (c) R/Z = 0.10 and (d) R/Z = 0.20. The curves in (a) are plotted using (3.6) and (3.8).

-25 - location and the magnitude of the traveling wave lobe. The corresponding peak values of IPI (in dB) are shown in Fig. 3.7(a). The analogous results for resistive strips having R/Z = 0.05, 0.1 and 0.2 are shown in Figs. 3.6 and 3.7, (b) through (d). The point to be observed is that as R increases from zero the minimum strip width for which a traveling wave exists also increases. Thus, for R/Z = 0.05 (corresponding to 19 ohms), a traveling wave is present only for w/x > 1.9, and for smaller w/x the peak which occurs at almost the expected position (see Fig. 3.6(b)) is actually a side lobe of the specular flash. This is evident from the manner in which the peak location changes with increasing w/x, and the logical explanation is that for w/x < 1.9 the small amount of loss has reduced the traveling wave peak below the level of the side lobe. When R/Z = 0.1 there is no traveling wave lobe unless w/x > 3.4. The way in which the peaks for w/x < 3.4 are "taken over" by the side lobe of the specular flash is graphically illustrated in Fig. 3.7(c), and when R/X = 0.2 there is no evidence of a traveling wave for any w/x < 4. It is apparent that even a small amount of loss is sufficient to eliminate a traveling wave for modest values of w/X. On the other hand, if the purpose of the non-zero resistivity is to reduce the near edge-on cross section, the presence of the peaks regardless of their origin could still be of concern. If we simply take the peak which occurs closest to edge-on, skipping from (side lobe) trajectory to trajectory as necessary, and plot its magnitude as a function of w/x, the results shown in Fig. 3.8 are obtained. For R/Z < 0.15 the points form a reasonably smooth oscillatory curve for all w/x > 1.4, but when R/Z = 0.2 the consequence of changing trajectories is

I0,b a (a).0 0 0 0.*0 0 @0. 0 * 0 0 a 0 0 0 0 0 0 0 0 *o X X 0.,: 0. 0 XX X0XXXX2XX9AXIX @0 @0 0 0 0 0 00 0 (b).0 0 0 00 0 0 0.0.0. 0*0 0 0 0 0 9 0 --I0 10 -Ip 0-.0.0 0 *. 0 0 0. 0 00 0 030 0 OX* XXXX * *iX 0 0 (C) *0 0000 0 (d) 0 0 O 0.00 N).0 0 99000 -I0 -0 0 0 000 40 0 0.0 0 0 0**0 00 0.0.. 0. 00 0..000 0 o. o0 0 0 0 0 0 00 0 0 0 0 0..0 0 0 0 *0 0 X Xx aX 0 0 0.0 ' 0 0 0 a 2 W4k 3 4 0 -TI Wr 2 / I3 -I 4 Fig. 3.7: Magnitudes IPI in dB of the peaks for (a) R = 0, (b) R/Z = 0.05, (c) R/Z = 0.10 and (D) R/Z = 0.2. The curves are merely to guide the eye and the crosses correspond to the traveling wave.

-27 - 1.0. x K (a ) K --- —. --- -! x K K Xx - x x x x X KX 0* 75'0*5 - Kx X K N "I K \Y.., K "I K ~ Kx x ~ - -- X X -— 4~.. —.K K xX -- —, X K -- ----- (C x N — A x x -)L~ t x (d) K ~~-% ---- -X.-~ x K 0.25J 0 2 W/A 3 4 Fig. 3.8: Magnitudes of the actual or nearby traveling wave peak for (a) R = 0, (b) R/Z = 0.05, (c) R/Z = 0.10 and (d) R/Z = 0.20. The curves were obtained using (3.9).

-28 - to produce noticable discontinuities. As expected, a non-zero resistivity decreases the magnitude of the peak, and the reduction below the value for a perfectly conducting strip increases with increasing R/Z and w/X. For R/Z = 0.2 and w/x = 4 the reduction is approximately 12 dB. We have not succeeded in developing a theoretical formula for the traveling wave peak showing the dependence on R/Z and w/x. Though a uniform rigorous second order GTD expression for the bistatic scattered field of a uniform resistive strip is available (Senior 1979b), the result for H polarization is discontinuous in the limit R * 0 and cannot be used to show the effect of resistivities as small as those of interest here. Nevertheless, from an examination of the data in Fig. 3.8 it appears that an expression of the form IP| = A + B (3.9) is adequate to predict the average return as a function of w/x for a given R. By a lease squares fit it is found that R/Z = 0, A = 0.88, B = 0.23 = 0.05, = 0.67 = 0.45 = 0.10, = 0.26 = 0.53 = 0.20, = 0.13 = 0.42 and it is reasonable that these values be used in conjunction with (3.9) to estimate the effect of a small amount of resistivity.

CHAPTER IV. RESISTIVE PLATES The strip or ribbon is a two-dimensional analogue of a finite plate, and there is a considerable amount of useful information about the scattering from a lossy plate that can be obtained by considering the simpler two-dimensional problem. Nevertheless, the information is limited in its applicability, and in cases such as near-grazing incidence when the side edges of a plate play a role it is necessary to consider the plate directly. 4.1 Background A major activity has been the development of an effective and efficient code to compute the scattering from a finite, planar resistive plate of infinitesimal thickness when illuminated by an incident plane wave. We remark that the restriction to a resistive plate is not, in fact, a restriction at all. As pointed out in a recent publication (Senior, 1985: included here as Appendix B), a thin layer whose permittivity and permeability both differ from the free space values of the surrounding medium can be simulated using superposed (electrically) resistive and (magnetically) "conductive" sheets, and these sheets are uncoupled when they are planar. A conductive sheet is the electromagnetic dual of a resistive one, and thus by running the program twice with the polarization and the sheet parameters appropriately defined, the solution for the most general imperfect plate can be obtained by addition. A special case of a combination sheet is the opaque plate having an impedance boundary condition imposed at its surfaces. -29 -

-30 - Several years ago a program was developed (Naor and Senior, 1981) to solve the integro-differential equations for the components of the total electric current induced in a resistive plate when illuminated by a plane wave. The program employs rectangular subdomains and uses simple differencing to carry out the differentiations numerically, and because of this the accuracy and generality are less than what we would like. Based on our recent experience with a program to treat dielectric plates at low frequencies (Ksienski, 1985) where a static analysis is appropriate, it was felt that the numerical differentiation could be avoided, and that a more efficient, accurate and general program would result. The efficiency is important because of our desire to treat plates up to a square wavelength or two in area with a matrix of only modest (e.g. 128 x 128) size, and this was one factor that led to our use of the C language. It is also important that the program accommodate plates of general shape (and this motivated the choice of arbitrary triangular subdomains), having an arbitrarily specified non-uniform resistivity when illuminated by a plane wave incident in any direction with any polarization. It is believed that these objectives are met by the program we have developed. 4.2 Program Description The plate is assumed to lie in the plane z = 0 of a Cartesian coordinate system x,y,z, and is simulated by a resistive sheet whose resistivity R (ohms per square) is R = kt(Ir-1) (4.1) are the thickness and relative permittivity of the plate where t and r are the thickness and relative permittivity of the plate material, k and Z are the propagation constant and intrinsic impedance

-31 - of free space, and a time factor e it is assumed and suppressed. The boundary conditions at the plate are x+E z x = 0 and z x E = Rz x J where R = z x R| is the (EFIE) total electric current, for 3 is then and an electric field integral equation A = i -nc + i~ kZ ^ + 1 S R z x J = z x + - ) J G dS 4=7r k2 (4.2) where G is the free space Green's function G = (e ikr)/r. An equivalent version of (4.2) can be derived using vector and scalar potentials and is R z x J = z x EkZ + z x J G dS' - - z x fv p/E G dS' (4.3) S ith with v* J = P/e i k (4.4) and this is actually the form used by Naor and Senior (1981). In both instances the differentiation was carried out numerically. To avoid this, we first split (4.2) into Cartesian components as follows:

-32 - i kr i kr i'L3Z +3 XJ k'L 4~rs 3, dS' + -i; aa k (4i i 5) R3, E- 4 s i vr iz. VZ i e_-~ dS% + wff VI Inc +I --- y r R3 =EY 4, i kr 2.! e_ d + 3 r y2 S x 3 x a ay the equations bcm wihen discretized~ and~ gut in matrix form (inc, + ~AIOX) + ~Bl (3 YI (4.6) i inc) + Il (3Jx) + PI (3 J (Ey f?I(3,I wher( j n gtinn the vector Cnan i nc an i-eement Column ech subdomaini, i C) is inNedent field on e e X(Y~) x(y) compo "ent Of i nc notaini g the Yen ~lm vector C ) s an ji-elemntcoum at the (J, (l IY mponent Of~ the induced cU sampling pointS5 cntiin h resi'tiv't~e ~R'1is n Nx tA matrix n at the samplin e muts~tual impedances containing th m x t matrices Il"lllZ are 'nsbetween the subdomal ~ le~ ~iix equation as witen asasigeat -fhe equations (.)cnb inc - I Y Ey ztA-element colum~n matrix r - I - 5 I (Di - ~cl I 2%-el ement ZNx N matrix c~lumn matr-Y

-33 - so that ( 31 - - [[A][R L ] l [C] I [D]-[R] \ -En I (4.8) The matrix inversion is carried out using standard IBM FORTRAN library routines employing LU decomposition contained in the MTS NAAS package. The discretization is based on triangular subdomains, with "tent" subsectional basis functions. The latter produce a current which is linearly varying over each subdomain and piecewise continuous over the plate. The concepts are illustrated in Figs. 4.1 and 4.2. If the currents at the vertices Vl, V2 and V3 of a triangular subdomain are Fig. 4.1: Simple example of triangular subdomains for a triangular plate. A12 II VI (,yIY.) V2 (2,,y,) Fig. 4.2: "Tent"V3 (current basis functions.) Fig. 4.2: "Tent" current basis functions.

-34 - I1, I2 and I3 respectively, then I1 = a x + p y + y 1 1 1 1 I2 = c x + y + Y 2 2 2 2 2 I3 = a x + ~ y + Y 3 3 3 3 3 and the current over the entire subdomain is I(x,y) = (a + a + a )x + (3 + X + y )y + X + + 1 2 3 1 2 3 1 2 3 The resistivity is specified at the vertices of the subdomains, and assumed to be linearly varying over the subdomain and piecewise continuous over the plate. Special care is necessary when computing the matrix elements in (4.8) and the novel features of the program are concerned with this. In the progenitor (static) program (Ksienski, 1985b) the vertices of the subdomains were used as the sampling points, and because the singularity of the kernel of the integral equation is integrable, there is no difficulty in computing the field of the current over each subdomain at the vertices. For the dynamic problem, the higher order singularities make this approach impossible, and after examining a variety of generalized function theory methods (e.g., Fikioris, 1965; Lee et al, 1980; Asvestas, 1983; Miron, 1983), an alternative procedure was developed. For each subdomain including the self cell, the kernel was rationalized in the manner shown below and the sampling points were chosen at the centroids of the triangles. The program still computes the current at the vertices, but does so by linearly interpolating the values of the current at the centroids of the adjacent subdomains. Thus,

-35 - in Fig. 4.3, the current at V1 is found by interpolating the known currents at C1, C2, C3 and C4 found from the integral equations. / = Fig. 4.3: Current computation based on the centroids. As much as possible of the integration and differentiation was done analytically to maximize the speed and accuracy of the computation. To this end, the singular parts of the dynamic kernel which are most troublesome numerically were treated analytically, and the integrals over the subdomains were evaluated by conversion to line integrals using the method of Wilton et al (1984). Thus, the kernel was written as eikr 1 k r ei 1 k2r - 1+ ik - +k2r + ik 2 r - r 2 r r 2 and the numerical integration and differentiation was limited to the bracketted terms. The resulting equation for the x component of the current at the jth centroid is as follows: inc ikZ 1~ (1 i kr RJ Enc + 2 Jxi ir + ik -k r)dS' (evaluated analytically) i S. 1 ikZ eVkr 1 k2r + i4 J i (eikr - ik + k2r)dS'(evaluated numerically) i S. 1 (cont.)

-36 - i x 47rk ZJ aX2 J1 4 k ZJ xi i S. 1 xi (1+ ik - k2r)dS' (evaluated analytically) x r 2 (k - - - ik + kr jdS' (evaluated numerically) ~x2 r 2 axr +4k xy s i ( 1 + ik - k2r)dS (evaluated analytically) i S. 1 + cZk 1. 4Tr e 1 k2r (evaluated numerically) yi axy ikr -1 ik + 2 )dS' (evaluated numerically) The equation for Jyj is similar. Note that the derivatives are taken inside the integrals only when the integrands are non-singular. The numerical integration was carried out using the method of Hammer et al (1956), which is exact for a fifth order polynomial. Because of time limitations the program is not yet complete. The final part necessary to compute the scattered field of the plate has not been finished, nor does the program make use of any symmetries that the plate may possess. In its present form the program is limited to the computation of the current induced in a plate of arbitrary size and resistivity. The source code consists of approximately 2200 lines of C language code, plus approximately 1000 lines of IBM FORTRAN NAAS library routines. The program compiles and, to judge from the few runs that have been made on a VAX, appears to produce accurate results. Nevertheless, there is further work that must be done and apart from the

-37 -additions necessary to compute the scattered field we are also aware of some changes that should be made to improve the accuracy. A source listing is included as an attachment to this report.

APPENDIX A: IBM PC PROGRAM P-RIB-H A program has been developed for use on an IBM PC computer to determine the high frequency bista'ic scattered field of a perfectly conductinn strin or ribbon for H nolarization. A perfectly conducting strip of width w occupies the region 0 < x < w, -~ < z < X of the plane y = 0 of a Cartesian coordinate system x,y,z, and is illuminated by an H-polarized plane wave having i A -ik(x cos 0 + y sin o0) H = z e where a time factor e iwt is assumed and suppressed. At large distances the scattered magnetic field can be written as Rs z ei(kp-r/4) p(, where p,p are cylindrical polar coordinates with x = p cos q, y = p sin p. In terms of the two-dimensional far field amplitude P the scattering cross section per unit length in the z direction is Ca('o) = lP(,o) 2 () and in the particular case of backscattering, o0 =:. From symmetry it is sufficient to consider only i7/2 < < fn where P = ir corresponds to grazing incidence. Using a GTD approach, Senior (1979b) developed an expression for the bistatic scattered field of a uniform resistive strip through second order terms, and when specialized to the case of a perfectly conducting strip with p0 = q, the result is -38 -

-39-, = 1 + cos + i 1 - cos e -2ikw cos P( ) -4 cos 4 cos 2 1/2 e-17/4 eikW(l - cos) (2) 2 e(i)/4 e -+ 0([kw]-1) (2),rkw jsin + where the origin of phase is at the left hand edge of the strip. At broadside (p = i/2) the infinities of the first two terms cancel and (2) reduces to Ii, 1 i (2 \1/2 i(kw-7/4) -1 P(2 '2) kw - ) e )+ O([kw] )(3) For ( > 7r/2 the first and second terms on the right hand side of (2) are the contributions of the front and rear edges respectively, with the former vanishing for grazing incidence. The third term is the second order contribution and this clearly fails when p = zr, but by expressing (2) in terms of the half plane current and then using a uniform asymptotic representation of the current, a uniform expression for P(p,f) valid for 7T/2 < < rT is found to be F 1/2 ) = + cos (2 i4+ikw( - cos ) P( - 4 cos - kw 2 + 1 - c e2k cos 2 -/4 sin F (/2-kw cos-) + O([kw]-1) (4)

-40 - where F(T) is the Fresnel integral 00 F() = e du. (5) T For large values of T, i i 2 F() - 2- e and when this is substituted into (4) we recover the asymptotic expression (2). On the other hand, for small T F(T) = - eLvH 4 + 0(r) showing that P(T,7) = 0, as expected. In general F(T) = 1 C( 2)+i -[2 S(T2)] (6) where C and S are cosine and sine integrals whose computation is described by Boersma (1960). Thus 2 e i/4 sin - F(J/2kw cos () = sin {1 - C(u) - S(u) + i[C(u) - S(u)]} f~ 2 2 where u = kw(l + cos ). The introduction of the Fresnel integral to provide a uniform behavior in the vicinity of ( = ff produces a difficulty near broadside

-41 - since the infinities of the first two terms on the right hand side of (4) no longer cancel precisely when p = w/2. For this reason it is desirable to introduce a second Fresnel integral whose main effect is to produce a uniform behavior near ( = 0 in spite of the fact that we are only concerned with ( > i/2. The result is P(('_) = - Tt o cos F ( sin ) 4 cos i 2 \ + i 1- cos ( e-2ikw cos 2 -i4 sin, F (2k cos + 0([kw]-1) (7) valid for 0 < (p < r. The infinities at ( = r/2 now cancel precisely, and we remark that (7) is identical to the result obtained by asymptotic expansion of the uniform expression of Khaskind and Vainshteyn (1964). A program designated P-RIB-H (perfectly conducting-ribbon-H polarization) has been written to compute the function P(p,p) using (7). It is coded in BASIC on our IBM personal computer and computes a/x (in dB), JPJ and arg P (in degrees) as functions of the angle ( for 91 < ( < 179 degrees. The broadside angle ( = 90 degrees is omitted to avoid numerical problems, but in practice, a knowledge of the far field at 91 degrees is sufficient to determine the broadside return. When ( = 180 degrees, I|P = 0, and this angle is also omitted to avoid problems in computing arg P. Some results obtained with P-RIB-H are compared with those of an integral equation solution by the moment method in Figs. 1 through 6.

-42 - The agreement is good for all w/X > 0.5, though we do observe a somewhat larger phase discrepancy than expected for p > 130 degrees with w/x = 2. The traveling wave lobe is faithfully reproduced and this enables us to use (7) or, indeed, (4) for studying traveling wave effects. From the GTD viewpoint the Fresnel integrals were introduced to match the second order expansion for 0 < ) < 7r into the zero values at grazing incidence, but one interesting consequence is that (7), as opposed to (2), is reasonably accurate even for small values of w/x.

-43 - 6 - 4 - 2 - 0 0 \ 0 0 o Oo 90 120 150 180 4 (degrees) Fig. 1: IPI for w/x = 2 computed using P-RIB-H (ooo) compared with the results of an integral equation method ( —).

-44 - 0 00 \ ~~oo o -100- \ 90 0 0 \0 -200 90 120 150 180 0 (degrees) Fig. 2: Arg P for w/x = 2 computed using P-RIB-H (ooo) compared with the results of an integral equation method ( —).

-45 - 3_ 2 - 1 -0 \ 0 ----- I - 90 120 150 180 0 (degrees) Fig. 3: Arg P for w/x = 1 computed using P-RIB-H (ooo) compared with the results of an integral equation method ( ).

-46 - 0000 o O o -100 - s \o \0 bOo - 200 90 120 150 180 0 (degrees) Fig. 4: Arg P for w/x = 1 computed using P-RIB-H (ooo) compared with the results of an integral equation method (-).

-47 - 90 120 150 180 0 (degrees) Fig. 5: Arg P for w/X = 0.5 computed using P-RIB-H (ooo) compared with the results of an integral equation method (-).

-48 - 0 -oo00 bO - 200 - - 90 120 150 180 0 (degrees) Fig. 6: Arg P for w/x = 0.5 computed using P-RIB-H (ooo) compared with the results of an integral equation method ( ).

-49 - iF., -..'......H..:P f" ':N j.. f -'J-' Ir fb 'T"IIj F i' fI i I I I. _. -:, f I.::f. 4:_.;j I 1.1- L VIb....... Il- r % f2; fj.]f I — Sr '.: I U I JE F f::w WIFEE ~, ii I Z1.1 UA - F"C T''l'.... Aj I 1i. I ~ i) I 4::. L~ I3: K f I::.i Ii f2i HT FO VFt T I 'II: % ] I F '114E: 1'I 7 I I ('fI T Tff11N F'I I' IJ F 1-i l -ri T S C.J ili i T~...ff... I.I~IIC~f %I I::IzT D X. V... tI J11 1 1 __ f:.jfi l OlS gr F3 N L D::l.JIE: ii:I 1 ' I Y C ( 1 I IN J iILr L — ) ~ff.i LP1

-50 - ~,.'] i l,...............Z I _' I.......: ':. L::..:.i;:;'..i "::: l riY 41 r LTH F' PIT N T I L 4. -R TS IN::.. 4 1XE DCI NDT S TO AO I J 6 1i..N FA Cl. i!...:F'.]Iv'",..:.:.:,iE....E!i':iGi"...~-:' l 'l [~,q v,''.', -;: F::I l:: L-{' E;::':: Fjj~ Id 1F dI 6 TI DR z../. j':' "i"] j',::':,!!:'- F -:.'.C E( F-x 662 II i,6:,0} F'OR:; X':::::"+F 1/:i. i-3C:' TO F:I]-P::.I/:!ISO.': SE P;iiEF:: DX z~~ I. ('. 710 ''71 1 F, - I i' J 3 ':'9':::' ( I -: C............ T.. T... 3,. ZI I " I X - I N J3 -.J FF"1- CO ( ) '7:2,::.: i:'-i:L:: Ei, 70 S:.2 (.,:,: "':0...... I L( ~.1. T....:... Fi I ]J IT EF O M PFI iT i- F-I i I R E fE1 -.I E F! i- I T F F L LF TV! f C F F7 'E I\-I I (7- 1 T i.. i I:: -I....: l ' I -.:. PFC-, L. `1.E%1T C (::, f:-''A "'I" Ei F;-: F-!"',. I2 i',! I3:J; S30 FTFCX I-....... I..... ~5 FR2R......., F ' l-.... F ),6:5: IF T N E 1.:= r-.-4X2... 67 Cl ' FiR S~~~Z 7:.G SB10: -..; t......-~......:.,...;..,:...,',. 1 77 ".~~~.......:..... a. l 5........ -.~...._..,...7 9i 2 0:':! F: PI...:P. 'Tt'0 le!, X:: 3:; iE A} I"i" F::'S ]: /:-'- MENT EA"" V R 0 F R:::: E iS Si I t E-... N X 0, V YR ~9-.40:1 ALE AN OT.E AD:.O, OPTTO..:,~. 4 6,E,-,.'.:i.... I~.- L~ —`1.J~ 950':O.... 970'!"",:X]IAINR PAR l', THDBV TGA 69:02-:.,;'-;"...,, 016::57 ZY-'-5 -3 050 856.1 ZY-i-i i"' i:.' 07575249!' ZY't 7!"i'UT " —i1 B50663781 c?!5,::::,: 1', (],, ' - *{*5~7: ';{,!/29B^^ZY***5 -":: f::':: (3 E. ^IY *6 -:2;: (.'::, 291 — W ***7 - 033 9 7 ',:: 050:3 SQ (X:Z: Y:;!-''( Ii',I(:AF:~'.Y F: (:: Z; 'i *i']F 'i'HtOS ( Z:;)'; iI,!'EG)~I.,:::, RE,;...

-51 - T P Z fft D........... 'Ii '1....n.........,.~ I-.... ~ i 'I IC;7 16I E.1 F.41.~I ~ ~ fj p4( d ~;v F:H(tSEz:::~.I~Y ~sc~v~i.i:ME'F( )):~ I.8C: I III f~ji ii'Ii i~ ''ijfff:f~f...f:.....If.... e~ (:,Jkj

,A I'LN i;!;, IL. Li Tradrsdctions on Antenrlns dnld P'ropa(dtion, Vol. AP-33, No. 5, May 1985 ufi ted, L- uutlol (1) iiiplhex thlat th1er is no mllagletic current aJId, hence, Il!it the peitfeability of tlie corrctsonding layer ii t(he Sastle as that of tile suIIounduI g (l'lCe pa.c~) IICdIUII. because of this, (3) can be written as x { X ':( ) + )(-)1 = -2^, (4) Combined Resistlve and Conductive Sheets THOMAS B. A. SENIOR, PFLLOW, I.PS, Abslrwf-T iable at W Iblk ayr of merial whoee perilmlvliy aad permeability )bot dilfer fro the vlats for The *srromidig meliiI,. combiailomll rlsivi codiclcv* sheet Is defi eolleadbt Ilr properties decrtibed. 1. INTRODUCTION Thin layers of lossy material are of obvious interest for cross scction reduction purposes. A mathematical model of such a layer is a resistive sheet, and during the last few years the scatter11g froli this type of sheet has been extensively explored [l -.3I. The electromnagnetic dual of an electrically resistive sheet is a "'magnetically conductive" one 141, and to simulate a thin layer whose permittlvity and permeability both differ from the surrounding medium, it may be necessary to include this sheet as well. The properties of a combined sheet consisting of coincident resistive and conductive ones are examined. Although the two sheets are, in general, coupled, with each affecting the scattering wiiin the other, decoupling occurs when the sheets lie in a plaiie. '1his is true regardless of the (planar) configuration, and has irnplications for the development of analytical and numerical methods to predict the scattering from lossy plates. 11, BOUNDARY CONDITIONS An electrically resistive sheet is simply an electric current sheet whose strength is proportional to the tangential electric field at its surface, and in recent years the concept of such a sheet has found niany useful applications. As noted by Levi-Civlta (see 14, p. 191), its electromagnetic properties are completely specified by its resistivity R in ohms per square, and the boundary condiionis at its surface are and in the special case R 0 the sheet is perfectly conducting, whereas itf t 0i, the sheet is no longer present. For a material of targe conductivity o, R = (r)where r is the thickness of the layer, and an alternative expressioii is 151 -^ (5) valid if I Er - 11 I 1. In (5),k andZ(=l/Y) are the propagation constant and iuipedance, respectively, of free space, e, is the relative permittivity of the layer material, and a time factor e-"" hias been assunled. The electrolmlagnetic dual is a magnetically conductive sheet 161 supporting only a magnetic current J*, and by analogy with (1), (2), and (4) tie boundary conditions at its surface are n x I(+) - X I(-) = 0 (6) (-X k('-) = -J with (7) X {n x 1H(+)+ H(-) } = -2R*.* (8) where RI is the conductivity. By virtue of (6), a layer modeled by such a sheet must have its relative permittivity unity, and from (5) an expression for RH is iY = — (9) (Mr - l)0 valid for 1i, - I i 1I where i, is the relative permeability of the layer 11aterIal. 11I. COMBINATION SHEET To model a layer whose permittivity and permeability both differ tromn their free space values a logical approach is to consider a conbination sheet consisting of coincident resistive and conductive sheets at which the boundary conditions are (2)(4) and (6-(8). if s and t are unit vectors in the plane of the sheet such that at every point s- - 0 and X i = n, implying that s, i, and A form a right-handed system, (3), (4), (7), and (8) can be written as L'(+) + k,(-) = -2R{(H(+)- H-)),;(+) - E,(-)= -(2R)- '{Jt(+) + I/I(-)}; ',4(+) + E,(-) = 2R {H,(+) -H,(-)}.,(+) - A,(-) = (2R)-' {,(+) + H,(-)). Ax t(+) - A XE(-) = 0 (I) (2) wlthl (3) Addition and subtraction of the equations in pairs gives where; is a unit vector normal drawn outward to the positive (plus) side of the sheet and J is the total electric current supMaluscript received October 4, 1984. I he ullrt r with the jKdiuaion Laborlury, Department of Ellectrcal higKiLieeCrin nd Computle Science, Univcralty of Michiganu Ann Arbo, Ml 4 I 1.,(t)- - 1 ),(+) + * I ) ' H(-) 4R /\ 4R'.',(t)=N(I H/,(+)-R (I - H(- ) 4 * 4R, ( I (10) -52 -

showing that, in general, the sheet is partially transparent, but if 4~r 1 (11) the conditions (10) reduce to Leontovich boundary conditions 171 on the two sides of the sheet. The combined sheet is then opaque and is simply an impedance (boundary condition) sheet with the same surface impedance - 2R on each side. In view of (1 ) an equivalent formula for the surface impedance is q = (R/ R*) /2, and when the expressions for R and Re are inserted, we rind / /-1 2 1 /2 - r- ) Z a! Z ) (12) tuEr as expected. Under most circumstances the resistive and conductive sheets comprising a combination sheet are coupled inasmuch as the strength (as measured by the current supported) and the scatterng of each sheet are effected by the presence of the other. To illustrate this fact, consider a closed cylindrical sheet illuminated by an t'-polarized wave. In terms of the cylindrical polar coordinates p, 0, z, the sheet is defined as the surface p = a and the incident plane wave is assumed to have P.'/ - ize i k cos By expanding the scattered and interior fields in cylindrical wave functions, the solution for each type of sheet can be obtained by mode tnatching. In particular, the interior field is stn =t i E En(-i)"AnJn(kp) cos no (13) n — O whlere ae = 1 and en = 2, n > 0, and for a resistive sheet nka Z An= + - - Jn(ka)H(L)(ka) (14) 2 R [ior a conductive sheet A,, i. + - -— R k)HI'(ka) (15) where the prime denotes differentiation with respect to ka, but for a combination sheet n 1 - - 1 --- —-- + -- -J. nka)H(0)kg) 4RR* 4RR* 2 R ()n nkas Y 2 R* R If (II) is satisfied the field interior to the combination sheet vanishes, showing that the sheet is then opaque, and we also observe that (14) and (15) can be obtained from (16) by taking the lilits R* - oo and R -* a, respectively. More importantly, however, the coefficient An for the combination sheet is not the sum o,1 the coefticients for the individual sheets, which demonstrates the coupling. IV. PLANAR SHEET Although the sheets comprising a combination sheet are generally coupled, an exception occurs in the special case of a planar sheet. To prove this, let E - 1 + E(,) + Ek2) where?(t) and k(2) are the electric fields radiated by the induced electric and magnetic currents, respectively. If /(' ) and H(2) are the associated magnetic fields, the symmetry about a planar magnetic sheet implies i X {fH(2)(+) - H(2)(_)} = and nX ()(*) =~i*. Hence AX {E(2)(+) + (2)(-)} - O and when these expressions are substituted into the boundary condition (4) for a combination sheet, we find f x {h x [E() + E(')(+)+ E(-)+ E()(-)1} = -Rn x {H )(+) - ~ )(-)} which is simply the condition for the corresponding resistive sheet in isolation. Similarly, (8) reduces to the condition for a conductive sheet by itself, showing that for a planar combination sheet the resistive and conductive parts scatter independently of one another. This is true for a plate of any (planar) configuration and it is therefore sufficient to restrict the development of analytical and/or numerical procedures to the simple case of a resistive plate. By application of the duality principle the solution for the corresponding conductive plate can be deduced, and the solution for the combination plate then follows by addition of the component solutions. Thus, for a planar plate, we can achieve the added generality of a combination sheet without any increase in complexity. The resulting plate is partially transparent unless ( 1) is satisfied, which represents the special case when the impedance boundary condition is applied. REFERENCES (Ii T. B. A. Senior, "Scattering by resistive strips," Radio Sd., vol. 14, pp. 911-924. Sept./Oct. 1979. 12] --- "Backscattering from resistive strips," IEEE Trans. Antennas Proposat., vol. AP-27, pp. 808-813, Nov. 1979. [31 T. B. A. Senior and V. V. LUpa, "Backscattering from tapered resistive stris," IEEE Trans. Antennas Propagat. vol. AP-32, pp. 747-751. July 1984. 14 14. Baemann, Electrical and Optical Wave Motion, Cambridge, England Camhbridge Univ. Pre", 1915. (5) -. F. Haringion and J. R. Mautz, "An impdance sheet approximation for thin dielectric shells," IEEE Trans. Antennas Propogat., vol. AP23, pp. 531-534, July 1975. 61. TB.A. Senior, "Some extensions of Babiae's principle in electromagnetic theo)-y, IEEE Trans. Antennas Propagat., vol. AP-25, pp. 417 -420, May 1977. [71 --- "Impedance boundary conditions for imperfectly condacdt surfaccn, ' App. Sci. Res., vol. 8, sec. B, pp. 418-436, 1960. -53 -

References Abramowitz, M. and I. E. Stegun (1964), Handbook of Mathematical Functions, National Bureau of Standards Appl. Math Series 55, U.S. Government Printing Office, Washington, D.C,; p. 321. Asvestas, J. S. (1983), "Comments on 'Singularity in Green's function and its numerical evaluation'", IEEE Trans. Antennas Propagat. AP-31, 174-177. Boersma, J. (1960), "Computation of Fresnel integrals", Math Tables and Other Aids to Comp. 14, 380. Chang, S. and V. V. Liepa (1967), "Measured backscattering cross section of thin wires", University of Michigan Radiation Laboratory Report No. 8077-4-T. Fikioris, J. G. (1965), "Electromagnetic field inside a current carrying region", J. Math. Phys. 6, 1617-1620. Hammer, P. C., O. J. Marlowe and A. H. Stroud (1956), "Numerical integration over simplexes and cones", Math Tables and Other Aids to Comp. 10, 130-137. Khaskind, M. D. and L. A. Vainshteyn (1964), "Diffraction of plane waves by a slit and a tape", Radio Eng. Electron. 10, 1492-1502. Ksienski, D. A. (1985a), "A method of resolving data into two maximally smooth components", Proc. IEEE 73, 166-168. Ksienski, D. A. (1985b), "Scattering by distributions of small thin particles", Ph.D. dissertation, University of Michigan. Lee, S. W., J. Boersma, C. L. Law and G. A. Deschamps (1980), "Singularity in Green's function and its numerical evaluation", IEEE Trans. Antennas Propagat. AP-28, 311-317. -54 -

Maliuzhinets, G. D. (1958), "Excitation, reflection and emission of surface waves from a wedge with given face impedances", Sov. Phys. Dokl. 3, 752-755. Miron, D. B. (1983), "The singular integral problem in surfaces", IEEE Trans. Antennas Propagat. AP-31, 507-509. Naor, M. and T.B.A. Senior (1981), "Scattering by resistive plates", University of Michigan Radiation Laboratory Report No. 018803-1-T. Peters, L., Jr. (1958), "End-fire echo area of long, thin bodies", IRE Trans. Antennas Propagat. 6, 133-139. Ruck, G. T., D. E. Barrick, W. D. Stuart and C. K. Krichbaum (1970), Radar Cross Section Handbook (Plenum Press, New York); p. 132. Senior, T.B.A. (1979a), "Backscattering from resistive strips", IEEE Trans. Antennas Propagat. AP-27, 808-813. Senior, T.B.A. (1979b), "Scattering by resistive strips", Radio Sci. 14, 911-924. Senior, T.B.A. (1981), "The current induced in a resistive half plane", Radio Sci. 16, 1249-1254. *Senior, T.B.A. (1985), "Combined resistive and conductive sheets", IEEE Trans. Antennas Propagat. AP-33, 577-579. Senior, T.B.A. and V. V. Liepa (1984), "Backscattering from tapered resistive strips", IEEE Trans. Antennas Propagat. AP-32, 747-751. *Senior, T.B.A. and S. J. Yang (1984), "Traveling waves on thin bodies", Electronics Lett. 20, 1050-1051. *Volakis, J. L. and T.B.A. Senior (1985), "Simple expressions for a function occurring in diffraction theory", IEEE Trans. Antennas Propagat. AP-33, 678-680. -55 -

Wilton, D. R., S. M. Rao, A. W. Glisson, D. H. Schaubert, O. M. Al-Bundak and C. M. Butler (1984), "Potential integrals for uniform and linear source distributions on polygonal and polyhedral domains", IEEE Trans. Antennas Propagat. AP-32, 276-281. *S Supported by present contract. -56 -

Jun 6 11:27 1985 Page 1 Attachment 1 2 3 4 5 6 7 8 9 10 11 12 13 14 16 16 17 18 19 20 21 22 23 24 26 26 27 28 29 30 31 32 33 34 36 36 37 38 39 40 41 42 43 44 45 46 47 48 49 60 51 52 63 54 66 566 /* function definitions */ #include <math.h> double darea(),area(), sumang(); /* macro definitions */ #define Darray(Ara,I1,I2,L1,L2) Ara[(I1)+(I2)*(L1)] /* Macro to generate other\ macros to mimic fortran\ arrays */ #define Mnumpoi 100 /* maximum number of points */ #define Mnumtri 100 /* maximum number of triangular patches */ #define 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 */ #define True 1 /* value of logical true */ #define False 0 /* value of logical false */ #define Tri(J1,J2) Darray(tri,J1,J2,Mnumatri,Mnumpoi) /* set up tri as two\ dimensional array with limits (Mnumatri, Mnumpoi) */ #define Poil(Jl,J2) Darray(poil,J1,J2,Mnumatri,Mnumpoi) /* set up poil as two\ dimensional array with limits (Mnumatri, Mnumpoi) */ #define Poi2(J1,J2) Darray(poi2,J1,J2,Mnumatri,Mnumpoi) /* set up poi2 as two\ dimensional array with limits (Mnumatri, Mnumpoi) */ #define Poi(Jl,J2) Darray(poi,Jl,J2,Mnumtri,3) /* set up poi as two dimensional\ array with limits (Mnumtri,3) */ #define Matrix(Jl,J2) Darray(matrix,J1,J2,Mnumpoi,Mnumpoi) /* set up matrix as\ two dimensional array with limits (Mnumpoi,Mnumpoi) */ #define Cmatrix(Jl,J2) Darray(cmatrix,Ji,J2,Mnumpoi,Mnumpoi) /* complex version of Matrix, used with the complex structure */ #define Cvmat(Jl,J2) Darray(cvmat,Jl,J2,2*Mnumtri,2*Mnumpoi) /* define Cvmat, used to store centroid-vertex contributions */ #define Vvmat(Jl,J2) Darray(vvmat,J1,J2,2*Mnumpoi,2*Mnumpoi) /* define Vvmat, used to store vertex-vertex contributions */ struct complex { double real; double imag; }; /* /* /* End of Defir /* /* int flagsym; /* int flagsyms; /* I itions / Beginning of external variable declarations *r/ 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=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:

Jun 5 11:27 1986 Page 2 57 0 - source patch is original patch (presumably first quadrant) 68 1 - source patch is in second quadrant 69 2 - source patch is in fourth quadrant 60 3 - source patch is in third quadrant */ 61 int eof; /* end of file variable (=0 means end of file) */ 62 int potflag; /* flag marks whether or not potentials have been computed 63 0 - potential has not been computed 64 1 - potential has been computed assuming x excitation 65 2 - potential has been computed assuming y excitation 66 3 - potential has been computed assuming z excitation */ 67 struct { 68 double x[Mnumpoi]; 69 double y[Mnumpoi]; 70 struct complex res[Mnumpoi]; 71 struct complex [2*Mnumpoi]; 72 struct complex jtotx[Mnumpoi]; 73 struct complex Jtoty[Mnumpoi]; 74 point; /* point is a structure which contains the locations of all 75 of the points used in defining the triangular patches. A point 76 number 0 is permitted as C defines arrays beginning with 0. 77 X and y are the coordinates of the points, J is the computed 78 current resulting from a single excitation with x components 79 listed first and y components displaced by Mnumpoi. jtotx and 80 Jtoty contain the x and y components of the total current. */ 81 82 struct { 83 int numatri[Mnumpoi]; 84 int Tri(Mnumatri,Mnumpoi-1); 86 int Poil(Mnumatri,Mnumpoi-1); 86 int Poi2(Mnumatri,Mnumpoi-1); 87 } apoint; /* apoint is a structure which contains lists of triangles 88 associated with each point. numatri contains the number of 89 triangles associated with each point; tri contains the list of 90 triangle numbers associated with each point; and poil and poi2 91 contain the other two vertices used in defining the triangle. 92 poll and poi2 are indices to point. tri is an index to 93 triang. 94 struct { 96 struct { 96 double x[Mnumtri]; 97 double y[Mnumtri]; 98 } centroid; 99 int Poi(Mnumtri,3-1); 100 double poten[Mnumtri]; 101 struct complex cpoten[Mnumtri]; /* complex version of poten */ 102 ) triang; /* triang is used to solve the scattering problem with z 103 excitation. centroid.x and.y contain the coordinates 104 of the centroid of each triangle. poi contains the list of 105 points (vertices) associated with each triangle, and is an 106 index to point. poten is the potential of each triangle as 107 determined from solution of the matrix problem. */ 108 double Matrix(Mnumpoi,Mnumpoi-l); /* This is "THE" matrix. Trivial points 109 (ie, those which have zero potential from 110 symmetry considerations) are skipped when 111 reading or filling the matrix, which is 112 always done sequentially. */

Jun 6 1t:27 1985 Page 3 113 struct complex Cmatrix(Mnumpoi,Mnumpoi-1); /* complex version of Matrix */ 114 struct complex Cvmat(2*Mnumtri,2*Mnumpoi-1); 115 struct complex Vvmat(2*Mnumpoi,2*Mnumpoi-1); 116 struct complex ecvec[2*Mnumtri]; 117 struct complex evvec[2*Mnumpoi]; 118 double fvect[Mnumpoi]; /* fvect is the forcing vector for the matrix problem 119 and after solution of the matrix problem contains 120 the solution vector. */ 121 struct complex cfvect[Mnumpoi]; /* complex version of fvect */ 122 int mnumpoi; /* This is the total number of points. Points are assumed to be 123 numbered sequentially from 0. mnumpoi <= Mnumpoi. */ 124 lnt mnumtri; /* This is the total number of triangles. Triangles are assumed 125 to be numbered sequentially from 0. mnumtri <= Mnumtri. */ 126 union { 127 char str[2]; 128 char let; 129 } com; 130 /* This is the first character of each input line, and is 131 interpreted as a command. Valid commands are: 132 d - display linkup of points and triangles 133 and display potentials if defined 134 e - specify direction of exciting field, 135 and solve the resultant matrix problem. 136 g - graph potentials of points/triangles 137 h - enter heading, ie, a descriptive one line title. 138 m - enter material parameters of plate: t and tau. 139 p - enter coordinates of next point 140 r - regenerate matrix problem using finer mesh 141 s - define symmetry to be assumed in solving problem 142 t - enter definition of next triangle */ 143 char hedstr[82]; /* contains one line heading, description of data */ 144 double dYk; 145 double dYZ = 378.7; /* impedance of free space */ 148 double t; /* Thickness of the plate */ 147 double tau; /* permittivity of the plate */ 148 double taui; /* imaginary part of the permittivity. This variable is also 149 used as a flag: if taui is zero, computations are performed 160 assuming *taun is purely real; if taui is non-zero, the routines 151 appropriate for a complex tau are invoked */ 152 double dipmom; /* contains the real part of the computed dipole moment */ 153 double dipmomi; /* contains the imaginary part of the computed dipole moment */ 154 double Epsilon = le-9; /* A very small number, used for approximate equality */ 155 double Pi = 3.1415926535; /* Pi */ 166 double dYxxr,dYxxi,dYxyr, dYxyi,dYyxr,dYyxi,dYyyr,dYyyi; /* used in interface */ 157 double theta, phi, alpha; /* specifies the direction of propagation (theta 158 and phi), and the polarization of the electric 159 vector with alpha. alpha is specified in the 160 same manner as phi. */ 161 double rtheta, rphi, ralpha; /* contains angles in radians */ 162 /*****************************************************************************/ 163 /* 164 /* M A I N P R 0 G RAM */ 165 /* */ 166 /* This is a program which calculates scattering by an arbitrarily shaped, */ 167 /* thin, flat, resistive plate. Excitation may be specified in the x, y, */ 168 /* or z directions. The plate is made up of an arbitrary number (nominally */

Jun 6 11:27 1985 Page 4 169 /* less than 50) of triangular patches, upon which a method of moments */ 170 /* solution is obtained. Contributions from each patch is calculated via */ 171 /* surface integrals which are evaluated analytically, thus contributions */ 172 /* from each patch is obtained exactly (at least to 10 plus digits). The */ 173 /* only approximations which are made are that the potential varies linearly */ 174 /* with z inside the plate, and the division of the plate into triangular */ 175 /* patches, inside of which the field varies linearly with x and y. The */ 176 /* shape and size of the triangular patches are completely arbitrary, and it */ 177 /* is noted that the patches need not be contiguous. */ 178 /* / 179 /***************************************************************************/ 180 main() { 181 hedstr0[] = '\0'; 182 for(eof=scanf("%ls",com.str);eof == 1;) { 183 switch (com.let) { 184 case 'P': 185 case ':p' 186 case ': 187 case 't': 188 rdata(); /* read in the data */ 189 continue; 190 case 'R': 191 case 'r': 192 rres(); /* read in the resistivity of each point */ 193 continue; 194 case 'H': 196 case 'h': 196 gethed(); 197 break; 198 case 'E': 199 case 'e': 200 suplup); /* read in direction of excitation 201 and solve the resulting matrix problem */ 202 default: 203 printf(Cc is an invalid command\n,com.let); 204 while(scanf('Xc",com.str), com.let 1= '\n'); 206 /* flush out bad command */ 208 break; 207 } 208 eof=scanf("%ls",com.str); 209 } 210 } 211 /***************************************************************************** 212 / */ 213 /* ROUTINE: darea(i) */ 214 / */ 216 /* darea is called by suplup to calculate the area of a triangle. darea */ 216 /* accomplishes this with a call to area. The argument i specifies the */ 217 /* triangle number. */ 218 /* */ 219 /***************************************************************************** 220 double darea(i) 221 int i; /* pointer to triangle */ 222 { 223 roturn(area(point.x[triang.Poi(i,0)] point.y[triang.Poi(1,0)], 224 point.x triang. Pol(i,1)], point.y [triang.Poi (i.1) ]

Jun 5 11:27 1985 Page 5 225 point.x[triang.Pol(,2)],point.y[trlang.Pol(i,2)])); 226 } 227 /******************************************************************************/ 228 /* */ 229 /* ROUTINE: suplup */ 230 / */ 231 /* suplup is called by the main program to read in excitation, decompose the */ 232 /* problem into several symmetric problems if appropriate and possible, and */ 233 /* compute the current induced on the resistive plate. */ 234 /*, */ 235 /******************************************************************************/ 236 suplup() { 237 int 1, J, i2; /* loop variables */ 238 scanf ('"le %le %le %le'",dYk,.theta,tphi,kalpha); 239 rtheta=theta*Pi/180.; 240 rphi=phi*Pi/180.; 241 ralpha=alpha*Pi/180.; 242 for (i=O;i<mnumtri;i++) { 243 ecvec[i].real = cos(dYk*(triang.centroid.x[il*sin(rtheta)* 244 cos(rphi)+trlang.centroid. y [i] *sin(rtheta)*sin(rphi))))* 245 sqrt(sin(rtheta)*sin(rtheta)*sin(rph) *sin(rph) +cos(rtheta)* 246 cos (rtheta))* 247 cos(ralpha); 248 ecvec[i].imag = sin(dYk*(triang.centroid.x[i]*sin(rtheta)* 249 cos(rphi)+triang.centroid.y[i] *sin(rtheta)*sin(rphi)))* 250 sqrt(sin(rtheta)*sin(rtheta)*sin(rphi)*sin(rphi)+cos(rtheta)* 251 cos(rtheta))* 252 cos(ralpha); 253 ecvec[i+mnumtri].real =cos(dYk*(triang.centroid.x[i]*sin(rtheta)* 254 cos(rphi)+triang.centroid. y[i *sin(rtheta)*sin(rphi)))* 255 sqrt(sin(rtheta)*sin(rtheta)*cos(rphi)*cos(rphi) +cos (rtheta)* 256 cos (rtheta))* 257 sin(ralpha); 258 ecvec[i+mnumtri].imag =sin(dYk*(triang. centroid. x[i] sin(rtheta)* 259 cos(rphi)+triang.centroid.y[iJ*sin(rtheta)*sin(rphi)))* 260 sqrt(sin(rtheta)*sin(rtheta)*cos(rphi)*cos(rphi)+cos(rtheta)* 261 cos(rtheta))* 262 sin(ralpha); 263 for (i2=O;i2<mnumpoi;i2++) { /* fill 'cmat */ 264 cpcon(i,i2); 265 > 266 } 267 cv2vv(); 268 solvem(2*mnumpoi); /* using Vvmat, evvec, and point.J */ 269 for (i=O;i<mnumpoi;i++) { 270 point.jtotx[i].real=point.J i].real; 271 point.jtotx[i].imag=point. [i].imag; 272 point.jtoty[i].real=point.j[i+mnumpoi].real; 273 oint. J toty[i]. imag=point.J [i+mnumpoi].imag; 274 275 printy(); 276 277 /******************************************************************************/ 278 /* */ 279 /* ROUTINE: cv2vv */ 280 /*

Jun 6 11:27 1985 Page 6 281 /* cv2vv is called by suplup to convert from the centroid vertex matrix to * 282 /* the vertex vertex matrix. 283/*1 284 /*************************************** 285 cv2vvo){ 286 mnt i,J,J2; /*loop variables*/ 287 double sarea,tarea; /* stores area of triangles*/ 288 for (j=0;J<mnumpoi;J++) { /* convert to point-point * 289 for CJ2=0,sarea=0.evvec[J].real~evvec[J].imag= 290 evvec[J+mnumpoi].real= evvec[J+mnumpoiJ.mag-O; 291 J2<apoint,.numatri(J];32++) { 292 tarea=darea(apoint.TriCJ2,J)); 293 sarea += tarea; 294 e vve.c J].real i-= tarea*ecvec[apoint.Tri(j2,J)].real; 295 evvec[J].imag += tarea*ecvec~apoint.TriCJ2,j)].imag; 296 evvec[J+mnumpoiJ.real += tarea* 297 ecvec apoint,.Tri~j2,j)+mnumtri].real; 298 evvec[J +mnumpoiJ. imag += tarea* 299 ecvec apoint. Tri~j2,j) +mnumtri].imag; 300} 301 evvec[JJ real /sarea; 302 evvec j.I~mag /sarea; 303 evvec[J+mnumpoi].real /sarea; 304 evvec[j+mnumpoiJ.Imag 1=sarea; 305 for Ci=0;i<nmnumpoi;i++) { 306 for (J2=Osarea=OVvmat(iJ).real=Vvrnat(i~j).imag307 Vvmat~i j +mnumpoi).ra=vatin~j+mnumpoi).imag308 Vvmat(i+mnumpoi.J).real=Vvmat~i+mnumpoi,J).imag309 Vvmat Ci+mnumpoi j +mnumpoi).real= 310 Vvmat Ci+mnumpoi~j +mnumpoi).imag-o; 311 J2<apoint.numatri~ij;j2++){ 312 tarea~darea (apoint.TrCJ 29 i); 313 sarea += tarea; 314 Vvmat~ij).real += tarea*Cvmat~apoint.TriCJ2,i),J).real; 315 Vvmat~i,J).imag += tarea*Cvmat~apoint.Tri(J2,i),J).imag; 316 Vvrnat(i1j+mnumpoi).rea1 += tarea* 317 Cvmat~apoint.Tri~j2,i).J+mnumpoi).real; 318 Vvmat(j+mnumpoiJi).imag += tarea* 319 Cvmat~apoint.TriCJ2.i).J+mnumpoi).imag; 320 Vvmat~i+mnumpoi.J).real += tarea* 321 Cvmat~apoint.TriCj2,i)+mnumtrl~j).real; 322 Vvmat~i+mnumpoi,j).imag += tarea* 323 Cvmnat(apoint.TriCJ 2,i)+mnumtri,j)A.mag; 324 Vvmat~i+mnumpoi~j+mnumpoi).real += tarea* 325 Cvmat~apoint.Tri(C 2,i)+mnumtri.j+mnumpoi).real; 326 Vvmat(i+mnumpoi~j+mnumpoi) imag += tarea* 327 Cvmat~apoint.TriCj2,i)+mnumtri~j+mnumpoi).imag; 328 } 329 Vvmat~i,J).real /sarea; 330 Vvmat(IlJ).imag /sarea; 331 Vvmat(19J+mnumpoi).real /sarea; 332 Vvmat(iJ+mnumpoi).imag sarea; 333 Vvmat(l+mnumpoi.J).real /sarea; 334 Vvmat~l+mnumpol,j). mag /sarea; 335 Vvmat~l+mnumpo,j +mnumpoi).real1 sarea; 336 Vvrmat~i+mnumpoi,j+mnumpoi).imag 1=sarea;

Jun 6 11:27 1985 Page 7 337 } 338 } 339 ) 340 /******************************************************************************/* 341 /* */ 342 /* ROUTINE: cpcon(i,i2) */ 343 /* */ 344 /* cpcon calculates the elements of the matrix in terms of the contribution */ 345 /* of a source point to a centroid. The source point contribution is broken */ 346 /* up into triangle contributions, which are analyzed by calling contrib. */ 347 /* cpcon is the only routine which calls contrib, and thus is the */ 348 /* interface between the control program, written by Dave Ksienski, and the */ 349 /* matrix element numerical evaluation routines, written by Joe Burns. cpcon*/ 350 /* also calculates the contribution to the self cell from the RJ term. */ 351 /* 352 /******************************************************************************/ 353 cpcon(i,i2) 354 int 1,12; /* field point i, source point i2 */ 355 < 356 struct complex rsum; /* holds values of resistance at centroid */ 357 int iloop; /* loop variable */ 358 for (Cvmat(i,i2).real=Cvmat(i,i2).imaCmagvmat(1,i2+mnumpoi).real= 359 Cvmat(l,i2+mnumpoi).imag=Cvmat(i+mnumtri, 2).real= 360 Cvmat(i+mnumtri,i2).lmag=Cvmat(i+mnumtri, 12+mnumpo).real= 361 Cvmat(i+mnumtri,i2+mnumpoi). ima=O, 362 iloop=O;iloop<apont. numatri [i2; loop++) { 363 contrib(triang.centroid.x[i],triang.centroid.y [i point.x [i2], 364 point.y[i2],point.x[apoint.Poil(iloop,i2) ] 365 point.y[apoint.Pol (iloop,i2)],potpoint.apont.Poi2(iloop,i2)], 366 point.y[apoint.Poi2(iloop,i2)]); 367 /* >> << */ 368 /* >> contrib returns values through the external variables dYabc, where << */ 369 /* >> a=x or y and represents the direction of the field < */ 370 /* >> b=x or y and represents the direction of the current << */ 371 /* >> c=r or i and represents either the real or imaginary part < */ 372 /* >> << */ 373 Cvmat (i,2). real+=dYxxr; 374 Cvmat(i,i2).i mag+=dYxxi; 375 Cvmat(i,i2+mnumpoi).real+=dYxyr; 376 Cvmat (i,i2+mnumpoi). imag+=dYxyi; 377 Cvmat(i+mnumtri,i2).real+=dYyxr; 378 Cvmat(i+mnumtri,i2). imag+=dYyxi; 379 Cvmat(i+mnumtrl. i2+mnumpoi).real+=dYyyr; 380 Cvmat(i+mnumtri, 12+mnumpoi). imag+=dYyyi; 381 if(apoint.Tri(iloop,i2) == i) { /* self cell contribution */ 382 rsum.real=(point.res [12].real+ 383 point.res[apoint.Pol (iloop,i2)].real+ 384 point.res[apoint.Poi2(iloop,i2)].real)/3; 385 rsum.imag=(point.res [2]. imag+ 386 point.res[apoint.Poil(iloop,i2)].imag+ 387 point.res[apoint.Poi2(loop,i2)].imag)/3; 388 Cvmat(i,i2).real+=rsum.real; 389 Cvmat(i,i2). imag+=rsum. mag; 390 Cvmat(i+mnumtri,i2+mnumpoi.real+=rsum.real; 391 Cvmat(i+mnumtri,i2+mnumpoi). imag+=rsum. mag; 392 }

Jun 5 11:27 1986 Page 8 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 416 416 417 418 419 420 421 422 423 424 426 426 427 428 429 430 431 432 433 434 436 436 437 438 439 440 441 442 443 444 445 446 447 448 } /******************************************************************************/ /* ROUTINE: rres /* */ /* rres is called by the main program to read in the resistivity associated */ /* with each point. The resistivity is assumed to be complex and linearly */ /* varying. Upon encountering a non-'r" command rres retruns to the main */ /* program. */ /* */ rres() { int i,j; /* loop variables */ int inumb; /* next point */ do { switch (com.let) { case 'R': case 'r': scanf ("%d,kitnumb); scanf(le %le",&point.res [inumb].real, kpolnt.resinumb].imag); continue; default: return; } } while (scanft (ls, com.str) == 1); } /******************************************************************************/ /* */ /* ROUTINE: ddata /* */ /* ddata is called by the main program to display data. ddata display */ /* locations of points, connections of points, and potentials of points if */ /* they have been computed. ddata also displays status of various flags */ /* which indicate symmetries of plate and direction of excitation. */ /* '/ printy () int i,; /* loop variables */ printf %s \n',edstr); printf('k=%lf, theta=%lf phi=%lf, alpha=%lf\n",dYktheta,phi.alpha); for (i=0; i<mnumpoi; i++). printf (point # %d is located at x = %lf, y = %lf\n',i,point.x[i], point.y[i]); printf(" and has current Jx=%lf +ilf, Jy=Xlf +ilf\n, point.jtotx[i].real,point.jtotx[i].imag,point. toty[i].real, point. toty(i]. imag); printf(" and has resistivity R=%lf +ilf\n", point.res[i].realpoint.res[i].imag); printf("point %d is associated with points and triangles (tri,P1,P2)\n,1); for (j=0; J<apoint.numatri[i];j++) { printf(" (d,%d,%d), ",apoint.Tri(J,i), apoint.Poil(j,i), apoint.Poi2(j,i));

Jun 6 11:27 1985 Page 9 449 if (j%6 == 6 II j == apoint.numatri[i]-l) printf('\n"); 450 } 451 } 452 } 453 /******************************************************************************/ 454 /* 455 /* ROUTINE: rdata */ 456 /* */ 457 /* rdata is called by the main program to read in point and triangle */ 458 /* definitions. Upon encountering a non- "p" or,"t' command, rdata calls */ 459 /* ldata to link the data, and then returns to the main program. 460 /* */ 461 **************************************************************************** 462 rdata() { 463 int i; /* loop variable */ 464 int inumb; /* number of next point or triangle */ 465 do { 466 switch (com.let) { 467.case 'P': 468 case 'p': 469 scanf (C'd",inumb); 470 scanf("%le %le",point.x+inumb,point.y+inumb); 471 mnumpoi=inumb+1; 472 continue; 473 case 'T': 474 case 't': 476 scanf ("d,kinumb); 476 for (i=0;i<3;i++) scanf("%d',ktriang.Poi(inumb,i)); 477 mnumtri=inumb+1; 478 continue; 479 default: 480 ldata(); 481 return; 482 } 483 } 484 while (scanf (C s", com. str) == 1); 485 } 48 /***************************************************************************** 487 /* */ 488 /* ROUTINE: ldata 489 /* */ 490 /* ldata links all information concerning points and triangles. This 491 /* information is needed to precisely define each patch on the plate. */ 492 /* */ 493 /******************************************************************************/ 494 ldata() 495 int i,J; /* loop variables */ 496 int isl; /* scratch variable */ 497 for (i=O;i<mnumpoi;i++) apolnt.numatri[i]=0; 498 for (=O;i<mnuntri;i++) ( 499 for (triang.centroid.x[i]=triang.centroid.y[i]=0,J=O;J<3;j++) { 500 isl=triang.Poi(i,); 601 triang.centroid.x i] += point.x[isl]/3; 602 triang.centroid.y[i] += point.y[isl]/3; 603 apoint.Tri(apoint.numatri[isl].isl) i; 604 apoint.Poil(apoint.numatrl[isl],isl) = triang.Poi(i.(J+1)%3);

Jun 6 11:27 1985 Page 10 605 apoint.Poi2(apoint.numatri[isl]++,isl) = triang.Po(i, (j+2)%3); 506 } 507 } 608 return; 609 } 510 /*****************************************************************************/ 511 /* 512 /* ROUTINE: solvem(N) */ 613 /* */ 614 /* solvem is used to solve the matrix problem. The matrix and forcing vector*/ 615 /* are created in xsolv or ysolv. N is the dimension of the vector. solvem */ 616 /* solves the matrix problem by calling the two fortran routines dgeco and */ 617 /* dgesl which preform a decomposition and back substitution. solvem is the */ 618 /, only c routine which calls a fortran routine. */ 619 /* */ 620 /******************************************************************************/ 621 solvem(N) 622 int N; /* N is the order of the matrix stored in the array matrix. N is less 623 than or equal to isize */ 624 { 625 lnt,j; 626 int isize,job; /* arguments for fortran routines, isize is the size of the 527 array matrix, job (=0) indicates type of matrix problem */ 628 int ipvt[2*Mnumpoi]; /* an array used by the fortran routines to store the 629 pivoting vector */ 630 double rcond; /* the condition number of the matrix. */ 631 double z[2*Mnumpoi*2]; /* scratch vector, lengthened for complex case */ 632 isize=2*Mnumpoi; 533 job=0; 634 for(i=0;i<13;i++){ 635 for(j=0;j<13; ++)( 637 536 printf('{lf %lf\n,Vvmat(i,J)); 538 } 639 for(i=0;i<N+l;i++) < 640 point.j [i].real=evvec[i].real; 641 point. [1].imag=evvec [].imag; 642 > 643 cgeco (vvmat,.isize,.N, ipvt,rcond,z); 644 printT('condition number is %e\n',l/rcond); 645 cgesl_(vvmat,isize,tN,ipvt,point. J,job); 646 647 /****************************************************************************/ 648 /* 649 /* ROUTINE: gethed */ 550 /* */ 6651 /* gethed is called by the main program to get the heading of the data. 652 /* gethed fills hedstr with the remainder of the line (the entire line 653 /* except for the letter H which must be in column 1. */ 654 /* */ 65 /***************************************************************************** 666 gethed() { 657 int 1; /* loop variable */ 558 for (i=0; (hedstr[i]=getchar()) 1= '\n'; i++); 559 hedstr[i] = '\0'; 660 return;

Jun 5 11:27 1985 Page 11 561 } 662 ***************************************************************************** 563 /* */ 664 /* ROUTINE: area */ 665 /* */ 666 /* area is called by graphp to calculate the area of a triangle. This is */ 667 /* needed in the computation of the area coordinates. */ 568 /* */ 569 ************************ ************************************************ 670 double area(xl,yl,x2,y2,x3,y3) 671 double xl,yl,x2,y2,x3,y3; /* These are the coordinates of the three 672 points which delimit the triangle */ 673 { 674 return(fabs(x2*y3-y2*x3-xl*y3+yl*x3+xl*y2-yl*x2)/2); 676 } 576 #include <matb.h> 677 double dYk,dYZ; 578 double Rk re,Rk imRkx reRkx im,Rky re,Rky im; 679 double Dx2rk re7Dx2rk Im,Dx2rRx re,Dx2rkx im,Dx2rky reDx2rky im; 680 double Dy2rk-re Dy2rk-im,Dy2rkx-re, Dy2rkx-im,Dy2rky-re,Dy2rky-im; 681 double Dxyrk-re,Dxyrk-im,Dxyrkx-re,Dxyrkx-im,Dxyrky re,Dxyrkyim; 682 double IT1,IT2,IT3,IT4,IT6,IT6,TT7,IT8,IT~; 683 double DX2IT1,DY2IT1,DXYIT1; 684 double DX2IT2,DY2IT2,DXYIT2; 586 double DX2IT3,DY2IT3,DXYIT3; 686 double DX2IT7,DY2IT7,DXYIT7; 687 double DX2IT8,DY2IT8,DXYIT8; 688 double DX2IT9,DY2IT9,DXYIT9; 689 extern double dYxxr,dYxxl,dYxyr,dYxyl,dYyxr,dYyxi,dYyyr,dYyyi; 690 691 692 contrib(xo,yo,xs 1, ys 1,xs2,ys2,xs3,ys3) 593 694 595 double xo,yo,xsl,ysl,xs2,ys2,xs3,ys3; 596 697 698 { 599 /*************************************************************************/ 600 /* */ 601 /* This routine calculates the finite element contribution to the x */ 602 /* component of the current or the y component for x or y excitation */ 603 /* given the vertices of the source triangle xslysl xs2,ys2 xs3.ys3 and*/ 604 /* the observation point. The finite element approximating the current */ 605 /* varies linearly over the triangular subdomain and is assumed to be */ 606 /* one at xsl,ysl and zero at xs2,ys2 and xs3,ys3. */ 607 /* */ 608 /*************************************************************************/ 609 /* */ 610 /* xo x coordinate of the observation point */ 511 /* yo y coordinate of the observation point */ 612 /* xsl x coordinate of vertice 1 of the source triangle */ 613 /* ysl y coordinate of vertice 1 of the source triangle */ 614 /* xs2 x coordinate of vertice 2 of the source triangle */ 615 /* ys2 y coordinate of vertice 2 of the source triangle */ 616 /* xs3 x coordinate of vertico 3 of the source triangle */

Jun 5 11:27 1985 Page 12 617 /* ys3 y coordinate of vertice 3 of the source triangle */ 618 /* a x coefficient of the finite element */ 619 /* b y coefficient of the finite element */ 620 /* c constant coefficient of the finite element */ 621 /* */ 622 /*************************************************************************/ 623 /* */ 624 /* EXTERNAL VARIABLES */ 625 /* */ 626 /************************************************************************/ 627 /* dYxxr real part of the x component of the current for x excitation */ 628 /* dYxxi imag part of the x component of the current for x excitation */ 629 /* dYxyr real part of the y component of the current for x excitation */ 630 /* dYxyi imag part of the y component of the current for x excitation */ 631 /* dYyxr real part of the x component of the current for y excitation */ 632 /* dYyxi imag part of the x component of the current for y excitation */ 633 /* dYyyr real part of the y component of the current for y excitation */ 634 /* dYyyi imag part of the y component of the current for y excitation */ 635 /* */ 636 /* IT1 is the integral of 1/R over the triangle */ 637 /* IT2 is the integral of x/R over the triangle */ 638 /* IT3 is the integral of y/R over the triangle */ 639 /* IT4 is the integral of 1. over the triangle */ 640 /* ITS is the integral of x over the triangle */ 641 /* IT6 is the integral of y over the triangle */ 642 /* IT7 is the integral of R over the triangle */ 643 /* IT8 is the integral of xR over thr triangle */ 644 /* IT9 is the integral of yR over the triangle */ 645 /* DX2IT1 is second derivative of IT1 with respect to x */ 646 7* DY2IT1 is second derivative of IT1 with respect to y */ 647 /* DXYIT1 is derivative of ITI with respect to x and y */ 648 /* DX2IT2 is second derivative of IT2 with respect to x */ 649 /* DY2IT2 is second derivative of IT2 with respect to y */ 650 /* DXYIT2 is derivative of IT2 with respect to x and y */ 651 /* DX2IT3 is second derivative of IT3 with respect to x */ 652 /* DY2IT3 is second derivative of IT3 with respect to y */ 653 /* DXYIT3 is derivative of IT3 with respect to x and y */ 654 /* DX2IT7 is second derivative of IT7 with respect to x */ 655 /* DY2IT7 is second derivative of IT7 with respect to y */ 656 /* DXYIT7 is derivative of IT7 with respect to x and y */ 657 /* DX2IT8 is second derivative of IT8 with respect to x */ 658 /* DY2IT8 is second derivative of IT8 with respect to y */ 659 /* DXYIT8 is derivative of IT8 with respect to x and y */ 660 /* DX2IT9 is second derivative of IT9 with respect to x */ 661 /* DY2IT9 is second derivative of IT9 with respect to y */ 662 /* DXYIT9 is derivative of IT9 with respect to x and y */ 663 /* Rk re real part of the rational kernel */ 664 /* Rk-im imag part of the rational kernel,/ 665 /* Rkx re real part of x times the rational kernel */ 666 /* Rkxilm imag part of x times the rational kernal */ 667 /* Rky-re real part of y times the rational kernal */ 668 /* Rky im imag part of y times the rational kernel */ 669 /* Dx2fk_re real part of second derivative with respect to x of */ 670 /* the rational kernel */ 671 /i Dx2rk-im imag part of second derivative with respect to x of */ 672 /* the rational kernel */

Jun 5 11:27 1985 Page 13 673 /* Dx2rkx re real part of x times second derivative with respect to x */ 674 /* of the rational kernel 676 /* Dx2rkx im imag part of x times second derivative with respect to x */ 676 /* of the rational kernel */ 677 /* Dx2rkyre real part of y times second derivative with respect to x */ 678 /* of the rational kernel */ 679 /* Dx2rky_im imag part of y times second derivative with respect to x */ 680 /* of the rational kernel */ 681 /* Dy2rkre real part of second derivative with respect to y of */ 682 /* the rational kernel */ 683 /* Dy2rk_im imag part of second derivative with respect to y of */ 684 /* the rational kernel */ 685 /* Dy2rkxre real part of x times second derivative with respect to y */ 686 /* of the rational kernel */ 687 /* Dy2rkxim imag part of x times second derivative with respect to y */ 688 /* of the rational kernel */ 689 /* Dy2rkyre real part of y times second derivative with respect to y */ 690 / of the rational kernel 691 /* Dy2rkyim imag part of y times second derivative with respect to y */ 692 /* of the rational kernel */ 693 /* Dxyrkre real part of derivative with respect to x and y of * 694 /* the rational kernel */ 695 /* Dxyrk im imag part of derivative with respect to x and y of */ 696 /* the rational kernel 697 /* Dxyrkx re real part of x times derivative with respect to x and y */ 698 /* of the rational kernel */ 699 /* Dxyrkx im imag part of x times derivative with respect to x and y */ 700 /* of the rational kernel */ 701 / Dxyrky re real part of y times derivative with respect to x and y */ 702 /* of the rational kernel */ 703 /* Dxyrky_im imag part of y times derivative with respect to x and y */ 704 /* of the rational kernel */ 705 /* dYk wavenumber 706 /* dYZ impedance 707 /* */ 708 *********************************************************************** 709 710 double a,b,c,k,k2.PI; 711 double ikrnl re,ikrnllm,lxyk_re,ixyk im,ix2kre,ix2k lm,iy2kre,iy2kim; 712, 713 k=dYk; 714 k2=k*k/2.; 715 PI=3.1415926535; 716 717 ************************************************************************/ 718 /* */ 719 /* Call routines numer and analyt to calculate external variables */ 720 /* */ 721 ****************************************************************** * 722 723 analyt(xo,yo,xsl,ysl,xs2,ys2,xs3,ys3); 724 725 numer(xo,yo,xsl,ysl,xs2,ys2,xs3,ys3); 726 727 /************************************************************************ 728 /*

Jun 6 11:27 1985 Page 14 729 /* Calculate coefficients of finite element 730 / */ 731 /************************************************************************/ 732 733 a=(ys2-ys3)/ ((xsl-xs2) *(ys2-ys3)-(xs2-xs3)*(ysl-ys2)); 734 735 b=(xs3-xs2)/((xsl-xs2)*(ys2-ys3)-(xs2-xs3)*(ysl-ys2)); 736 737 c= ((s2-xs3)*ys2-(ys2-ys3)*xs2)/((xsl-xs2)*(ys2-ys3)-(xs2-xs3)*(ysl-ys2)); 738 739 /*************************************************************************/ 740 /* */ 741 /* Calculate integrals of current over the triangle 742 /* */ 743 /***********************************************************************/ 744 745 746 ikrnl re=a*(IT2-k2*IT8)+b*(IT3-k2*TIT9)+ l-+c* k2*IT7) +a*Rkxre+b*Rky re 747 +c*Rk re; 748 749 ikrnl im=k*(a*IT5+b*IT6+c*IT4)+a*Rkx im+b*Rky im+c*Rk lm; 750 - - - 750 761 ixyk re=a*(DXYIT2-k2*DXYIT8)+b*(DXYIT3-k2*DXYIT9)+c*(DXYITI-k2*DXYIT7) 752 +a*Dxyrkxre+b*Dxyrkyre+c*Dxyrkre; 753 754 ixyk_im=a*Dxyrkxim+b*Dxyrkyim+c*Dxyrk_im; 755 756 ix2k re=a*(DX2IT2-k2*DX2IT8)+b*(DX2IT3-k2*DX2IT9)+c*(DX2ITI-k2*DX2IT7) 757 +a*Dx2rkx re+b*Dx2rky re+c*Dx2rk re; 758 759 ix2k im=a*Dx2rkx im+b*Dx2rky im+c*Dx2rk im; 760 761 iy2k re=a*(DY2IT2-k2*DY2IT8)+b*(DY2IT3-k2*DY2IT9)+c*(DY2IT1-k2*DY2IT7) 762 +a*Dy2rkx re+b*Dy2rkyre+c*Dy2rkre; 763 764 iy2k im=a*Dy2rkx im+b*Dy2rkyim+c*Dy2rk im; 765 766 766 /*************************************************************************/ 767 /* 768 /* Calculate contributions of the finite element currents 769 /* */ 770 /*************************************************************************/ 771 772 dYxxr=((-k*dYZ*ikrnl im)+(-dYZ*ix2k im/k))/(4.*PI); 773 774 dYxxi=((k*dYZ*ikrnl re)+(dYZ*ix2k re/k))/(4.*PI); 775 776 dYxyr=(-dYZ*ixyk im/k)/(4.*PI); 777 778 dYxyi=(dYZ*ixyk re/k)/(4.*PI); 779 780 dYyyr=((-k*dYZ*ikrnl im)+(-dYZ*iy2k lm/k))/(4.*PI); 781 782 dYyyi= ((k*dYZ*ikrnl re)+(dYZ*iy2k re/k)) / (4. *PI); 783 784 dYyxr=dYxyr;

Jun 5 11:27 1985 Page 15 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 dYyxi=dYxyi; analyt(xo,yo,xs 1, ys 1, xs2,ys2,xs3,ys3) double xo,yo,xsl,ysl,xs2,ys2,xs3,ys3; /* */ /* Calculate the contribution of a particular triangle given its vertices */ /* x[ll,y] x[2],y[2] x[3],y[3] and the field point xo,yo * /* List of Variables: */ /*, */ /* x[n] x coordinates of the vertices of the triangle */ /* y[n] y coordinates of the vertices of the triangle /* xo x coordinate of the field point */ /* yo y coordinate of the field point */ /* ai[n] slopes of lines through segments / alpi n] 1./sqrt(l.+a**2) */ /* gitn] ai[n]*qi[n] / /* bi[n] y-intercept of lines through segments */ /* li[n] length of ith segment */ /* ci[n] constant related to the parameterization /* exi[n] indicates sense of x direction along integration path */ /* eyi[n] indicate sense of y direction along integration path /* Dxin] derivative of x' with respect to path length parameter t */ /* D yin] derivative of y' with respect to path length parameter t /* A[n] term of R /* Bin] term of R */ /* C n] term of R */ /* C2[n] sqrt(C[n]) */ /* Rip[n] distance from observation point to plus end of ith segment */ /* Rim n] distance from observation point to minus end of ith segment /* tip n] path length to plus end of ith segment /* tim[m] path length to minus end of ith segment /* Ili n] integral of 1./R over a segment */ /* 12i[n] integral of t/R over a segment */ /* I3i[n] integral of t**2/R over a segment */ /* I4i[n] integral of R over a segment */ /* I6i[n] integral of tR over a segment */ /* I6i[n] integral of (t**2)R over a segment */ /* sign temporary variable */ /* t temporary variable /* R temporary variable /* xt[n] temporary variable /* yt[n] temporary variable /* ntrmx temporary variable /* ntrmy temporary variable */ /* dtrm tempoarav variable */ I I I I I I I I I /* delta[n] (4AC-BB)/(8C) */

Jun 6 11:27 1986 Page 16 841 /* DXIll is first derivative of Ili with respect to x */ 842 /* DYIli is first derivative of Ili with respect to y */ 843 /* DX2I1i is second derivative of Iii with respect to x */ 844 /* DY2I1i is second derivative of Ili with respect to y */ 845 /* DXYIl1 is derivative of Ili with respect to x and y */ 846 /* DXI2i is first derivative of I2i with respect to x 847 /* DYI2i is first derivative of I2i with respect to y */ 848 /* DX2I21 is second derivative of I2i with respect to x */ 849 /* DY2I21 is second derivative of I21 with respect to y */ 850 /* DXYI21 is derivative of I2i with respect to.x and y / 851 /* DXI31 is first derivative of 13i with respect to x 852 /* DYI3i is first derivative of I31 with respect to y */ 853 /* DX2I31 is second derivative of I31 with respect to x */ 854 /* DY2I31 is second derivative of 131 with respect to y / 856 /* DXYI3i is derivative of I31 with respect to x and y 866 /* DXI41 is first derivative of I4i with respect to x */ 857 /* DYI41 is first derivative of 141 with respect to y */ 868 /* DX2I41 is second derivative of I41 with respect to x */ 859 /* DY2I4i is second derivative of I41 with respect to y */ 860 /* DXYI4i is derivative of 14i with respect to x and y / 861 /* DXI6i is first derivative of I6i with respect to x */ 862 /* DYISi is first derivative of I61 with respect to y */ 863 /* DX2ISi is second derivative of 161 with respect to x */ 864 /* DY2I61 is second derivative of 11 with respect to y 865 /* DXYIS5 is derivative of 161 with respect to x and y * 866 /* DXI61 is first derivative of I61 with respect to x */ 867 /* DYI6i is first derivative of I6i with respect to y 868 /* DX2I61 is second derivative of I6i with respect to x 869 /* DY2I61 is second derivative of 161 with respect to y / 870 /* DXYI61 is derivative of 161 with respect to x and y */ 871 /* DXIT1 is first derivative of IT1 with respect to x 872 /* DYIT1 is first derivative of IT1 with respect to y,/ 873 /* DXIT7 is first derivative of IT7 with respect to x */ 874 /* DYIT7 is first derivative of IT7 with respect to y */ 876 /* temp temporary variable 876 /******************************************************************************/ 877 /* */ 878 /* EXTERNAL VARIABLES */ 879 /* 880 /***************************************************** *** 881 /* IT1 is the integral of 1/R over the triangle 882 /* IT2 is the integral of x/R over the triangle */ 883 /* IT3 is the integral of y/R over the triangle */ 884 /* IT4 is the integral of 1. over the triangle 885 /* IT5 is the integral of x over the triangle */ 886 /* IT6 is the integral of y over the triangle 887 /* IT7 is the integral of R over the triangle 888 /* IT8 is the integral of xR over thr triangle 889 /* IT9 is the integral of yR over the triangle 890 /* DX2IT1 is second derivative of IT1 with respect to x */ 891 /* DY2IT1 is second derivative of IT1 with respect to y 892 /* DXYIT1 is derivative of IT1 with respect to x and y */ 893 /* DX2IT2 is second derivative of IT2 with respect to x */ 894 /* DY2IT2 is second derivative of IT2 with respect to y */ 996 /* DXYIT2 is derivative of IT2 with respect to x and y 896 /* DX2IT3 is second derivative of IT3 with respect to x */

Jun 6 11:27 1985 Page 17 897 /* DY2IT3 is second derivative of IT3 with respect to y */ 898 /* DXYIT3 is derivative of IT3 with respect to x and y */ 899 /* DX2IT7 is second derivative of IT7 with respect to x */ 900 /* DY2IT7 is second derivative of IT7 with respect to y */ 901 /* DXYIT7 is derivative of IT7 with respect to x and y */ 902 /* DX2IT8 is second derivative of IT8 with respect to x */ 903 /* DY2IT8 is second derivative of IT8 with respect to y */ 904 /* DXYIT8 is derivative of IT8 with respect to x and y */ 905 /* DX2IT9 is second derivative of IT9 with respect to x */ 906 /* DY2IT9 is second derivative of IT9 with respect to y */ 907 /* DXYIT9 is derivative of IT9 with respect to x and y */ 908 /******************************************************************************/ 909 910 double ai[4] alpi[4],gi[4],bi[4],li[4],ci[4],exi [4] ey[4] Dxi[4],Dyi[4]; 911 double,Rim[4],tip[4],tim[4],x[5],y[5],Rip[4],A[4] B[4] C[4],C2[4]; 912 double Ili[4],I2i[4],I3i[4],I4i [4],I 5i [4I6i[4],di[4]; 913 double DXIli[4],DYIli[4],DX2Ili[4],DY2Ili[4] DXYIll[4]; 914 double DXI2i[4],DYI2i [4],DX2I2i[4], DY2I2i [4],DXYI2i[4]; 915 double DXI3i[4],DYI3i [4],DX2I3i[4],DY2I3i[4],DXYI3i[4]; 916 double DXI4i[4],DYI4i [4],DX24i[4],DY2I4i[4],DXYI4i[4] 917 double DXI5i[4], DYI51[4], DX2I5i[4], DY2I5i[4], DXYI5i[4]; 918 double DXI6i[4],DYI6i[4],DX2I6i [4],DY2I6[4], DXYI6i [4] 919 double DXIT1,DYIT1,DXIT7, DYIT7; 920 double sign,t[3] R[3], xt [3],yt[3], ntrmx,ntrmy,dtrm,delta[4], temp; 921 double sqrt(),log(),fabs(); 922 int n,m; 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 /******************************************************************** **** ******/ /* Calculate ai [n]. If segment is parallel to the y-axis, ai[n] would go */ /* infinity; consequently, set al[n] in this case to zero. This will not */ /* affect any calculations since this feature of ai[n] has been accounted */ /* for. */ /* */ /************************************************************** ****************/ x[l]=xsl; y[l]=ysl; x[2] =xs2; y[2]=ys2; x[3]=xs3; y[3]=ys3; x[4]=x[1]; y[4]=y[1]; for(n=l;n<=3;n++) { if(x[n+1] == x[n]) al [n]=0.; else ai[n]=(y[n+l]-y[n])/(x[n+l]-x[n]); /N***************************t i******************** ****************************/ /* */ /* Calculate alpi [n]. Note when ai[n] goes to infinity, qi [n] goes to zero. */

Jun 5 11:27 1985 Page 18 955 956 for(n=l;n<=3;n++) 957{ 958 if(x[n+l] == x[n]) 959.alpi~n]=0.; 960 else 961 al pi En]=1.0/sqrt(l1.0+ai En]*ai En]); 962} 963 964 /*************************************** 965 ** 966 / Calculate gl (n]. Note when ai~n] ge to inf inity, gil(n] goes to one. / 967/** 968/***************************************/ 969 970 for (n=1;n<=3;nl+g+) 971 { 972 if~x~n+l] == x~n]) 973 gi~n]=1.; 974 else 975 gi (nJ=ai (n] *alpi [n]; 976} 977 978 /*************************************** 979/*/ 980 /* Cal cu late bi [nb i [n],C (n],B n].A [nJ,ci (n] di (n] tip(n].tim (nJ,Rip(n] * 981 /*Rim En] 982 ** 983 /*************************************** 984 985 for Cn=1l;n<=3;n++) 9868 987 988 bi [nJ=y [nJ -ai[nJ *xn]; 989linqrCxn1-nJ*xn1-n]Cyn]ynJ*yn1-nJ; 990 i[]=qt(( n1I- n x[+Ix(]+( n - n yC+I- n 991 992 C En] =alpi En] *alpi En] +gi En] *gi En]; 993 C2Enl~sqrt(CC n]) 994 995 996 997 ci(Il=x(l]; 998 ciE[2]=xE[2]-alpiE(2]*li l[; 999 ciE(3] xE(3]-alpiE(3 *Cl1i [II +liE(2]) 1000 1001 di(Il=y l]; 1002 diE2]=yE2]-giE2]*1iE1]; 1003 diE[3]=yE3]-glE(3 * (1i (l]+liE(2]) 1004 1005 tip l]=li[l]; 1006 tipE2]=li1]1+liE2] 1007 tipE3]=1i(l]+1iE2]+l1.E3] 1008

Jun 5 11:27 1985 Page 19 1009 tim(l]=0.; 1010 tim[21]tipEI]; 1011 tim3] =tip[2]; 1012 1013 for(n=l;n<=3;n++) 1014 { 1015 1016 B[n]=(-2.0)*Calpi[n]*(xo-ci [n])+gi n]*(yo-di[nD) 1017 1018 A~n]=Cxo-ci n])*Cxo-ci[n])+Cyo-di[n])*Cyo-ditn]); 1019 1020 Rip En] sqrt (C En] *tip En] *tip En] iB En] *tip En] +A En] 1021 1022 Rim En] =sqrt (C En] *tim En] *tlm En] iB En] *tim En] -A En] 1023 1024 } 1025 1026 1027 1* * 1028 /* Calculate Dxi~n],Dyi~n] */ 1029 /* 1030 1031 1032 forCn=1;n<=3;n++) 1033 { 1034 1035 if( x~n+1] >= x~n] ) 1038 exicn]=1.0; 1037 ifC x~n+1] < x~n] ) 1038 exiln]=(-1.0); 1039 if( yEn -1 >= y~n] ) 1040 eyi~n]=1.0; 1041 if( yjn+l] < y~n] ) 1042 eyi~n]C-1.0); 1043 1044 DxiEn]=exiEn]*alpiEn]; 1045 1046 Dyi En] eyi En] *fabs (gi En]); 1047 } 1048 1049 1050 / 1051 /* Calculate integral of 1./R over a segment Iii */ 1052 /* Calculate integral of t/R over a segment 121 */ 1053 /* Calculate integral of t**2/R over a segment 13i 1054 /* Calculate integral of R over a segment 14i 1055 /* Calculate Integral of tR over a segment I i 1056 /* Calculate Integral of (t**2)R over a segment I61 1057 1058 1059 1060 for(n=l;n<=3;n++) 1061 1062 1063 IlinJ=(log(2.*C2En]*Rip [n]+2.*C[~ n*tippn]iBC nD-log(2.*C2En]*Rim~nli2. 1064 *C [n] *tim En] i- En] )[)/C2[n];

Jun 5 11:27 1985 Page 20 1065 1366 MiN]= (Rip (nJ-Rim~n])/C2[n]-B~n] *Iiin]/(2.*C~n1)) 1067 1068 I3i [nJ=((tip~n] -3.*B~n] /(2.*C~n])) *Rip En] -(tim~n] -3.*B~n] /(2. *C~nJ)) 1069 *Rim[n])+(3.*B~n]*B~n]/(4.*CEn]) -A~n])*Ili~n])/(2.*C~n]); 1070 1071 141 En]=((2.*C~n]*tip[n] +B[n])*Rip~n]-(2.*C[n]*tim[nj+B[n])*Rim~n])/(4. 1072 *C [n])+ (4. *A n] *C[n] -B [n] *B [n] )*Ili [n] /(8. *C [n]) 1073 1074 15 En) (powr (Rip En).3) -powr (Rim En].3)) 1 (3. *Ct n] ) -B En] *I41 En] / (2. *C En]) 1075 1076 IMi(n] =(((tip(n] -6.*B [n]/(6. *C n]))* owr (Rip[n],3) -(tim(n] -5.*B (n] /(6.*C[n] 1077 )*powr (Rim (n],3))+ (6.*B [n]*B n] /(4.*C (n] )-A [n] )*I4i[n])/(4. *CE]) 1078 } 1079 1080 /*************************************** 1081 / 1082 /*Calculate integral of 1hR over the triangle ITi 1083 ** 1084 /*************************************** 1085 1086 IT1=0.; 1087 for (n1l;n<=3;n++) 1088 IT1=ITI+(Dyi~n]*(alpi[n]*I2i~n]-(xo-ci~nD)*Iii~n])-Dxi~n]*(gi~n]*I2i~nI 1089 -(yo-di[n])*Ili~n])); 1090 1091 /*************************************** 1092 / 1093 1*Calculate integral of x/R over the triangle 1T2 * 1094 / 10956 **************************4************* 1096 1097 IT2=0.; 1098 for Cn=1l;n<=3;n++) 1099 IT2=IT2+(Dyi~n]*Calpi~nJ*alpi En]*I3i~n]+alpin En*C2.*ci (n]-xo)*I2i~nJ+ci~nJ 1100 * (ci [n] -xo) *I I I n] )-Dxi (n] *(gi (nJ*alpi(n] *I31i[n] +(5i (n]*ci (n] 1101 +alpi[n]*(di~nJ-yo))*I2i~n]+ci~nJ*Cdi~n]-yo)*Ili~n )); 1102 IT2=IT2/2.0+xo*IT1/2.0; 1103 1104 /*************************************** 110561 1106 / Calculate integral of y/R over the triangle 1T3 1107 1 1108 /*************************************** 1109 1110 IT3=0.; 1111 for(n1l;n<=3;n++) 1112 IT3=IT3+(Dyi~n]*(alpi~n]*gi~nJ*I3i~n]+(gi~n]*(ci~n]-xo)i-alpi~n]*di[n]) 1113 *I~i[n]+di~n]*(cifn]-xo)*Ili En] -Dxi~n]*(gil~n]*gi~n]*I3i~n]+gi~n] 1114 *(2.*di~n]-yo)*I2i~n]+di~n]*Cdi~n]-yo)*I1 i n])); 1115 IT3=IT3/2.0+yo*IT1/2.0; 1116 1117 /*************************************** 1118 / 1119 /*Calculate Integral of 1 over the triangle 1T4 1120 1

Jun 5 11:27 1985 Page 21 1121 /*************************************** 1122 1123 IT4=0.; 1124 for(n1l;n<=3;n++) 1125 IT4=IT4+ (Calpi En]*Dyi En] gi En]*Dxi En]) *(C p n] *tip En]-tim En]*tim En])/2.0 1126 + (Dyitn]*ci~n]-Dxi [n]*di En])*(tip~n] -timn En].0 1127 1129 e 1130 /*Calculate integral of x over the triangle ITS 1131/** 1133 1134 IT5=0.; 1135 for(n1l;n<=3;n++) 1136 IT5=IT5 —Dyi [n]*(alpi~n]*alpi En]*(powr~tip~n].3)-powr~tim~n],3))/3.0 1137 +alpiEn] *ci En] *(tip En]*tip En] -tim En]*timEn]) +ci En]*ci En] *(tip En] 1138 -timtn])))/2.0; 1139 1140 /*************************************** 1141 / 1142 /*Calculate integral of y over the triangle 1T6 1143 / 1144 /*************************************** 1145 1146 IT6=0.; 1147 for Cn1l;n<=3;n++) 1148 IT6=IT8+(Dxi En] *(gi En] *i En] * Cowr~tip En) 3) -powr(tim~n].3)) /3.0+ i En] 1149 *di~nJ*Ctip~n]*tip~nT -timEn] *tim En])+di En]*di En)*Ctip~n]-tim~n.i)))/2.0; 1150 1152 / 1153 /*Calculate integral of R over the triangle T 1154 / 1156 1157 IT7=0.; 1158 for~n=1;n<=3;n++) 1159 IT7=IT7+CDyi~n]*(alp i~n]*ISi~nJ-(xo-ci~nj)*I4i~nJ)-Dxi~nJ*Cgi~n]*I5i[nI 1160 -Cyo —di~nJ )*I4i~nJ)/30 1161 1162 /*************************************** 1163 */ 1164 /*Calculate integral of xR over the triangle 1T8 * 1165 / 1166 **************************************** 1167 1168 IT8=0.; 1169 for(n1l;n<=3;n++) 1170 IT8=IT8+ (Dyi (n] *(alpi[n] *alpi (n] *I8i(n] +alpi(n] *C2. *ci (n] -xo) *I5i (n] +ci (nJ 1171 *(ci~n]-xo)*I4i~n])-Dxi[n]*(gi~n]*alpi~n]*I6i[n]+(gi~n]*ci[n] 1172 +alpi~n]*(di~n]-yo))*I~i[n]+ci~n]*(di[n]-yo)*I4i[n)) 1173 IT8=IT8/4.0+xo*IT7/4.0; 1174 1175 /*************************************** 1 176 /*

Jun 6 11:27 1986 Page 22 1177 /* Calculate integral of yR over the triangle 1T9 1178 /* *1 1179 1180 1181 IT9=0.; 1182 for(n=l;n<=3;n++) 1183 IT9=IT9+(Dyi[n]*(alp i ngi~n]*I6i[n]+(gi[n]*(ci[n]-xo)+alpinn]*di~n]) 1184 *I5i[n]+di[nJ*(ci n] -xo)*I4i[n] )-Dxiln]*(gi[n]*gi[n]*I6i[n]+gi[n] 1185 *(2.*di[n]-yo)*I+ i n] +ditn*(di[n]-yo)*I4i [n])) 1186 IT9=IT9/4.0+yo*IT7/4.0; 1187 1188 1189 1* 1190 1* Calculate derivatives of Lii *1 1191 1* DXI1I is first derivative of Iii with respect to x *1 1192 1* DYIlI is first derivative of Iii with respect to y 1193 1* DX2IlI is second derivative of Iii with respect to x 1194 1* DY2I1I is second derivative of Iii with respect to y 1195 /* DXYII is derivatIve of Iii with respect to x and y 1196 1197 1198 1199 for~n=l;n<=3;n++) 1200 ( 1201 DXIlIinJ=O.O; 1202 DYIlin]=0.0; 1203 DX2IIinl=OJ.; 1204 DY2I1in]=O.O; 1205 DXYI1i[n=0.0; 1206 1207 t[i1=ttp[n]; 1208 tt2J=tim En]; 1209 1210 R[1J=Rlp[n J 1211 RE2]=RimnJ; 1212 1213 signl1.0; 1214 1215 for(m~l;m<=2;m++) 1216 { 1217 xt[mn=xo-Calpi[nJ*t[m]+ciEn]); 1218 yt[ml yo-Cgi n]*tmm]+dinn]); 1219 ntrmx=2.0*alpi (n]*C2[n]*Rlm]-2.0*CCn]*xt[m]; 1220 dtrm=2.0*C2[nJ*R[m]+2.0*C n]*t[m]+Bn]; 1221 ntrmy=2.0*ginn)*C2[nJ*R[m]-2.0*CmnJ*yt[mJ; 1222 1223 if(m==2) 1224 sIgn=(-1.0); 1226 1226 DXI1I in] =DXIIi [n] -sign*ntrmx/(C(n]*R(m]*dtrm); 1227 1228 DYIli[n]=DYIi En]-sign*ntrmy/CCn]*R/mJ *dtrn); 1229 1230 DX2IIi[n]=DX2Ili~n]4sign*(-Cntrmx*ntrmx)/(C[n]*sqrt(CCn])l*RE] 1231 *REm] *dtrm*dtrm)+(2.0/dtrm)*(1 (0/R[ml - xt[Em]*xt[Em)/ 1232 powr(R[m],3))):

Jun 6 11:27 1985 Page 23 1233 1234 DY2I11[n]=DY2I1i[n]+sign*(-(ntrmy*ntrmy)/(C[n]*sqrt(C[n])*R[m] 1236 *R[m]*dtrm*dtrm) + (2.0/dtrm)* (1.0/R[m]- (yt[m] *yt[m])/ 1236 powr(R[m],3))); 1237 1238 DXYI11[n]=DXYI1i[n]+sign*(-(ntrmx*ntrmy)/(C[n]*sqrt(C[n])*R[m] 1239 *R[m]*dtrm*dtrm)-(2.0*xt[m]*yt[m])/(dtrm*powr(R[m],3))); 1240 } 1241 1242 /******************************************************************************/ 1243 /* */ 1244 /* Calculate derivatives of 121,/ 1245 /* DXI21 is first derivative of I2i with respect to x */ 1246 /* DYI21 is first derivative of I2i with respect to y */ 1247 /* DX2I2i is second derivative of I2i with respect to x */ 1248 /* DY2I2i is second derivative of I2i with respect to y */ 1249 /* DXYI21 is derivative of I2i with respect to x and y */ 1250 /* */ 1251 /******************************************************************************/ 1252 1253 DXI2i[n]=(xt[l]/R[1]-xt[2]/R[2])/sqrt(C[n])-B[n]*DXIli[n]/(2.*C[n]) 1254 +alpi[n]*Ill[n]/C[n]; 1255 1256 DYI2i[n]=(yt[1]/R[1]-yt[2]/R[2])/sqrt(C[n])-B[n]*DYIi[n]/(2.*C[n]) 1257 +gi[n]*Ili[n]/C[n]; 1258 1259 DX2I2i[n]=((1.0/R[1]-xt[l]*xt[i]/powr(R[l],3))-(1.0/R[2]-xt[2]*xt[2]/ 1260 powr(R[2],3)))/sqrt(C[n])-B[n]*DX2Ili[n]/(2.*C[n])+2. 1261 *alpi[n]*DXIli[n]/C[n]; 1262 1263 DY2I2i[n]=((1.0/R[l]-yt[l]*yt(l]/powr(Rl],3))-(1.0/R[2]-yt[2]*yt[2]/ 1264 powr(R[2],3)))/sqrt(C[n])-B[n]*DY2Ili[n]/(2.*C[n])+2.*gi[n] 1265 *DYI2l[n]/C[n]; 1266 1267 DXYI2i[n]=(xt[2]*yt[2]/powr(R[2],3)-xt[l]*yt([1/powr(R[1],3))/C2[n] 1268 +gi[n]*DXIli[n]/C[n]+alpi[n]*DYIll[n]/C[n]-B[n]*DXYIli[n]/ 1269 (2.*C[n]); 1270 1271 /******************************************************************************/ 1272 /* */ 1273 /* Calculate derivatives of I3i */ 1274 /* DXI3i is first derivative of I3i with respect to x */ 1275 /* DYI31 is first derivative of I3i with respect to y,/ 1276 /* DX2I31 is second derivative of I3i with respect to x */ 1277 /* DY2I31 Is second derivative of 131 with respect to y */ 1278 /* DXYI31 is derivative of I3i with respect to x and y */ 1279 /* 1280 /******************************************************************************/ 1281 1282 DXI3i[n]=(((t[1]-3.*B[n]/(2.*C[n]))*xt[l]/R[1]+3.*alpi[n]*R[1]/C[n]) 1283 -((t[2]-3.*B[n]/(2.*C[n]))*xt[2]/R[2]+3.*alpi[n]*R[2]/C[n])) 1284 /(2.*C[n]); 1285 DXI3i[n]=DXI3i[n] + (3.*B[n]*B[n]/(4.*Cn] )-A[n])*DXIi[n]/(2.*C[n]) 1286 -(3.*B[n]*alpi[n]/(2.*C[n])+(xo-ci[n]))*Ili[n]/C[n]; 1287 1288 DYI3i[n]=(((t[l]-3.*B[n]/(2.*C[n]))*yt[ ]/R[1]+3.*gi [n]*R[1]/Cn])

Jun 5 11:27 1985 Page 24 1289 -C(t[2]-3. *B[n) /(2. *C[n]))*yt[2] /R12]1+3.*gi [n] *R[2] /C[n])) 1290 2.Cn; 1291 DYI3i[n]=DYI3i [n)+C3.*B[n]*B[n]/(4.*C[n))-A[n))*DYIii~n]/C2.*C[n]) 1292 -C3.*BI~n]*gi[n]/(2.*C[n])+Cyo —di~n]))*Ili[n]/C~n]; 1293 1294 DX2I31 [n]=(3.*alpi [n]*xt[1]/(C~n] *C[n] *R[1])-(t[1]-3.*B[n]/(2.*C[n])) 1295 *(1./R[1]-xt[1]*xt[l]/powr(R[l],3))/(2.*C~nD)-(3.*alpi En] 1296 *xt[2]/(C[n]*C[n]*R[2])+(t[2]-3.*B[n]/(2.*C[n]))*(l./R[21 1297 -xt[2]*xt[2]/powr(R[2L,3))/(2.*C[n])); 1298 DX2I3i[nJ=DX2I3i [n]+(3.*B[n]*B[n]/(4.*Ckn]) 1299 -A (n]I)*DX2 II i [n] /(2.*C-[n])- (3. *B~n] *alp i n] /C~n] -2. *Cxo 1300 -ci[n]))*DXIli[n]/C~n]+C3.*alpi-[n]*alpi~n]/C[n]-l.)*Ili~n]/ 1301 C~n]; 1302 1303 DY2I3i En]=(3.*gi n]*yt~i]/(CCn]*C~n]*RE1])+Ct~lJ-3.*B~n]/(2.*C[n])) 1304 *(1./R[1]-yt[1]*yt[1J/powrCR[iJ,3))/C2.*C[n]))-C3.*gi~n] 1305 *yt[2]/CC[n]*C~n] *R[2] +(t(2]-3.*B~n]/(2.*C~n]))*(1./R(2] 1306 -yt[2]*yt[2]/powr(R[2L,3))/C2.*C~n])); 1307 DY2I3i~n]=DY2I3i Ln]+(3.*B~n]*B[n]/(4.*C[n]) 1308 -A[n])*DY2Ili En]/(2.*C~nD)-(3.*B~n]*gi En]/C~nJ-2.*Cyo 1309 -di En]))*DYIli[n]/C[nJ+C3.*gi~n]*gi[n]/C[n]-l.)*Ili LnJ/C~nJ; 1310 1311 DXYI3i [nJ=3.*alpi~n]* t[1]/C2.*C(n]*C~n]*R[1iD-(tE1J-3.*B~n]/(2. 1312 *CtnJ) *ytEI *xt(1]/C2.*C[nJ*powrCR[1L,3))+3.*gi~n]*xt~lJ/(2. 1313 *Ctn]*C(n]4cR(1]); 1314 DXYI3i En)=DXYI3i En]-(3. *al p1 En] *t[2]1/(2.*C En]*C En]*R (2])-Ct [2] 1315 -3.*B[n]/(2.*CEn))) *yt(2] *xt[2]/C2.*C~n]*powr(RE2],3))+3. 1316 *gi[n]*xtE2]/C2.*C n] *C~n]*RE2])); 1317 DXYI3i~nJ=DXYI3i[nJ-(3.*B[n]*alpi[nJ/C2.*C~n]) 1318 +xo-ci [n])*DYIli~n]/C~nJ +C3.*B~nJ*B~n]/(4.*C~nJ)-A~nJ) 1319 *DXYIli~n]/C2.*C(n]); 1320 DXYI3i[nJ=DXYI3i~n]-C3.*B~n]*gi~n]/C2.*C~n])+yvo-di~n]) 1321 *DXI1i(n]/C~n]+3.*gi~n]*gi~n]*Iii~nJ/(C~nJ*C~n]); 1322 1324/** 1325 / Calculate derivatives of 14i 1326 /* DXI4i is first derivative of 14 with respect to x 1327 /* DYI4i is first derivative of 14 with respect to y 1328 /* DX2I4i is second derivative of 14i with respect to x 1329 /* DY2I4i is second derivative of I1i with respect to y 1330 /' DXYI4i is derivative of 141 with respect to x and y 1331 / 1332 /*************************************** 1333 1334 delta~n]=(4.*A~n]*C~n]-B~n)*B~n])/(8.*C[n]); 1335 1336 DXI4i[n]=(C(2.0*C~n]*t[1]+B~n])*xt(i]/R(1]-2.*alpi En]*R[1])-((2.0*C~n] 1337 *t[2]+B[n])*x-tE2]/R[2]-2.O*alpl[n]*R (2]))/(4.0*C[n]) +(xo 1338 -ci[n]+D[n]*alpi~n]/(2.0*C[nJ )) *111En]+delta [n>DXIli[n]; 1339 1340 DYI4i~n]=(((2.0*C[n]*t~l]+B~n])*yt(l]/R(i]-2.*gi[n]*R(1] )-((2.O*C~n] 1341 *t[2]+E3[n])*ytE2]/R[2]-2.O*gi~n]*R[2 )Y)/(4.O*C[n])+(ydin 1342 +B~n]*gi~n]/ (2.0*Ctn]))*Ili~n]+delta[n]*DYI1i[n]; y-in i4313

Jun 5 11:27 1985 Page 25 1345 s-B~n])*xt[1]*xt[l]/powr(Rji].3))-(C2.0*C~n]*t21] 1346 +B [n] -4. 0*alpi In]*xt[2] )/R [2] -(2. 0*C In]*t[2] +B~n])*xt[212 1347 *xtI[2]/powr (R [123)))(4.0O*C n]) 1348 DX2I4i~n]=DX2I4i~n]+(1.o -alpinln*alpiln]/Cnln) 1349 *Iii~n]+(2.0*(xo-ci[n])+B n]*alpiln]/C~n])*DXIli~n]+delta[nI 1350 *DX2Ili~n]; 1351 1352 DY2I4i In)=(((2.0*C~n]*t~l]+B~n]-4.0*ginln*yt~l])/R[1]-(2.0*C~n]*t~l] 1353 +B [n] )*yt (l]1*yt (l]/pOwr (RI[1].3) )- ((2. O*C [n] *t1] 1354 +B In] -4. 0*gi In] *yt [2] )/R 212]- (2.- O*C In] *t, [2] +8 In] ) *yt [2] 1355 *ytI(2]/powr (R[(2]3) ))/C4.0*C (n]); 1356 DY2I4i~n]=DY2I4i[n]+(1.0-ginln*giln]/C~n.]) 1357 *Ili~n]+(2.0*(yo-din [n)+B In]*giln]/C~n])*DYIli~n]+deltaLn] 1358 *DY2Ili[n]; 1359 1360 DXYI4i~nJ=((2.0*(g inln*x-t[2]+al pi[n]*ytI2])/R[2]+(2.0 1361 *C[n]*t212]+B[n])*xt[2]*yt[2 /powr(ft[2L,3))-(2.o*(ginl]*xt~l] 1362 +alpi In]*yt~i])/RI1] +(2.O*ClnJ *t[1]+B~n])*xt[1]*yt~l] 1363 /powrC(R l], 3) )) / (4.0O*C(n] ) 1364 DXYI4 i In] =DXYI4 I. n] -gi. In] *alpi, In] *1I1 i In] IC In] 1365 +(xo-ci[n]+B~n]*alpinln/(2.0*C~n]))*DYIli~n]+(yo-dinln+gilnI 1366 *B~n]/(2.0*C~n]))*DXIli~n]+deltaln]*DXYIli~n]; 1367 1368 /*************************************** 1369/** 1370 /* Calculate derivatives of 151 1371 /* DXI5i is first derivative of 151 with respect to x 1372 /* DYIM is first derivative of INi with respect to y * 1373 /* DX2I6i is second derivative of INi with respect to x * 1374 /* DY2I61 is second derivative of 151 with respect to y 1375 / DXYIM is derivative of 15 with respect to x and y 13768/ 1377 /*************************************** 1378 1379 DXI5i In)=(R 11] *xt11] -R212]*xt212]-B In] *DXI4i In]/2. 0+alpi In]*141[En)) 1380 /C~n]; 1381 1382 DYI~i~nJ=(R[1J*yt[1]-R[2)*yt[2J-B~n]*DYI4i~n]/2.0+giln]*I4i[n])/C[n]; 1383 1384 DX2I~i (n] =(R 11]I+xt (IlI*xt CIlI/R CII-R [12-xt(12*xt (12/R (12-B (n] *DX2I4i (n] 1385 202 al nIDX4in])C n 1386 /.+.*liI]*X~~]/~] 1387 DY2I~i 1nJ(R(11]+yt~l]*yt~l]/RI1]-R212-ytI2]*yt212]/R212]-BIn] *DY2I4i In] 1388 /2.O+2.0*gi In] *DYI4i In])/C~nJ 1389 1390 DXYI~i In]=(xt~l]*yt~l]/R 11]-xt[212]*yt[2]/R[12]+gi In] *DXI4i In] -B In 1391 *DXCYI4i Ln]/2.0+alpiln]*DYI4i~n])/C~n]; 1392 1393 /*************************************** 1394 / 1395 /* Calculate derivatives of 161 13968/ DXI6i is first derivative of 161 with respect to x 1397 /* DYI61 is first derivative of 161 with respect to y * 1398 /* DX2161 is second derivative of 161 with respect to x 1399 /* DY2I~i Is second derivative of 161 with respect to y 1400 / DXYI~i is derivative of 151 with respect to x and y *

Jun 5 11:27 1985 Page 26 1401/*/ 1402 /*************************************** 1403 1404 DXI61[(n] = (5.*alpi (n] *powrC(Rl [1,3) /C12. *C[n] *Cl n])+ (t [l]/ (4. *C fl]) -5. 1405 *B[n]/(24.*C[n]*C n])) *3.*R[1]*xt~l])-(5.*alpi~n]*powr(R[2L,3)/ 1406 C12. *C n] *C n])+C(t 2]/ (4. *C n] )-5. *B[n] /(24. *C[n] *C [n]))*3. *R122 1407 *t2) 1408 DXI6i~n]=DXI6i[n]-(6.*B[n]*alpi~n]/(4.*C[n]*C[n])+(xo-c1[nl])/(2.*C~nl])) 1409 *I4i [n]+(5.*B~n]*B~n] /(16.*C~n]*C~n])-A~n]/C4.*C~n])) 1410 *DXI4i[n]; 1411 1412 DYI61[n] =(5. *gi n] *powrCRli],3) /(12. *C n] *C (n])+C(t[(1]/4. *C[n] )-5. 1413 *B~n]/C24.*C~n) *C[n]))*3.*R[1] *ytIl])-C5.*gi En]*powr(RE2],3)/ 1414 (12.*C~n]*C~n])+Ct[2]/C4.*C~n])-5.*B[n]/(24.*C~n] *C~n]))*3.*RE2] 1415 *ytE2]); 1416 DYI6i[n]=DYI61 En]-(5.*B~n]*gi~n]/(4.*C[n]*C~n])+(yo-di[n])/C2.*C~n])) 1417 *14i [n]+(5. *Btn]*Bfn]/(16.*C~n]*C~nJ)-A~n]/(4.*C~n])) 1418 *DYI41En]; 1419 1420 DX2I61. n]C(5.*alpl [n]*R[1]*xt~l]/C2.*C~n]*C~n])+3.*(t~l]/C4.*C~n])-5. 1421 *B[n]/C24.*C~n]*C~n]))*CRE1]+xt~l]*xt~l]/R[1])); 1422 DX2I61 En] =DX2I61 En]-(5.*al p1En] 1423 *R (2]*xt E2]/C2. *Cn] *C (n] )+3. *(tE(2]/C4.*C [n] )-5. *B [n] /C24. 1424 *C En] *C En] )) * CR [2] +xt [2] *xt [2] /R [2])) 1425 DX2I~i (nJ=DX2I6i(n] -2. *C5. *B n] *alpi (n] / 1426 (4.*C~n]*C[n])4-Cxo-cl[n])/C2.*C~nJ))*DXI4i~n]-CI.-5.*alpi~nI 1427 *alpi(n] /C (n] )*I4i(n] /C2. *C n]); 1428 DX2I6i Cn]4=X2I6i (n] +(6. *B (n]*B [n] /16. *C [nJ*C (nJ) 1429 -A~n]/C4.*CED]))*DX2I41 En]; 1430 1431 DY2IBI~n]C5.*g En]*R[1]*ytEI]/C2.*C~n]*C~n])+3.*(t~l]/(4.*C~n])-5. 1432 *BEnl (/24.*C~n]*C~n]))*(R[1]+yt~l]*yt~l]/RE1])); 1433 DY2I6i En]=DY2I61 En]-(5. *gi En] 1434 *RE(2]*yt[2]/C2.*C[n]*C~n])+3.*(t[2]/C4.*C~n])-5.*B~n]/C24. 1435 *C~n]*C~n]))*(RE2]+ytE2]*ytE2]/RE2])); 1438 DY2I6i~n]=DY2I~i~n]-2.*C5.*B n]*gli~n]1(4. 1437 *C~n]*C~n])+Cyo-dl~n])/(2.*C~n]))*DYI4i~n]-(1.-5.*gi~n]*gi~nI 1438 /C~n])*I4I~n]/C2.*C~n]); 1439 DY2I6i~n]=DY2I6i~n]+C5.*B~n]*B~n]/(18.*C~n]*C~n])-A~n]/ 1440 (4. *C(n])) *DY2I4i En]; 1441 1442 DXYI6i~n]C56.*g En] *R[1]*xt~l]/(4.*C~n]*C~n])+5.*alpi En] *RE1]*ytE1]/ 1443 (4.*Cfn] *CEn])+3.*Ct~l]/(4.*C~n])-5.*B~n]/C24.*C~n]*C n])) 1444 *xt~l]*yt~l]/RE1]); 1445 DXYI8i~n]=DXYI6i n]-C5.*gi~n]*R[2]*xt[2]/(4.*C~n]*C~n])+5. 1446 *a~lpi~n]*RE2]*Tt[2]/(4.*C~n]*C~n])+3.*CtE2]/(4.*C~n])-5. 1447 *B~n]/(24.*C~n *C~n]))*xtE2]*ytE2]/RE2]); 1448 DXYI6i En] =DXYI6I En] -56. *alpl En] *gI En] 1449 *I4I[n]/(2.*C~n]*C~n])-(6.*B~n]*gi~n]/(2.*C~n])+yo-di~n]) 1450 *DXI41 [n / (2. *C n]) 1451 DXYI61 En] =DXYI6i En] - (5. *B En] *alpl En] 1 (2. *Cn )D +xo-ci En]) 1452 *DYI41 En]/(2.*C~n])+(6.*B~n]*B~n]/C16.*C~n]*C~n])-A[n]/C4. 1453 *C [n])) *DXYI4i En]; 1454 1455 1456

Jun 5 11:27 1985 Page 27 1467 1458 /* 1459 /* Calculate derivatives of ITi *1 1460 1* DXIT1 is first derivative of ITI with respect to x 1461 /* DYITi is first derivative of ITI with respect to y *1 1462 1* DX2ITI is second derivative of ITi with respect to x 1463 DY2IT1 is second derivative of ITi with respect to y 1464 1* DXYIT1 is derivative of IT1 with respect to x and y 1465 /* 1466 1467 1468 DXIT1=O.0; 1469 DYITIO.0; 1470 DX2IT1O.0; 1471 DY2IT1=O.0; 1472 DXYIT1O.0; 1473 1474 for(n=l;n<=3;n++) 1475 { 1478 DXIT1=DXIT1+Dyi~n]*(alpi[n]*DXI2i(nJ+Cci[nJ-xo)*DXI1I[n]-Ili[nJ) 1477 -Dxi[nJ*(gilnJ*DXI2i n]C+(diXn]-yo)*DXIliinD 1478 1479 DYIT1=DYIT1+Dyi nl*(alpi[n)*DYI2i[n]+(ci[n]-xo)*DYI1I En]) 1480 -Dxi~nl*(giln] *DYI2i(n]+(din]-yo)*DYIlinJ-IliIn]); 1481 1482 DX2IT1=DX2IT1+Dyin]*C(alpi En]*DX2I2i(n]+(ci[nJ-xo)*DX2IIi+nJ-2 1483 *DXlii[nJ)-Dxi[nJ*CgiinJ*DX2I2i[nJ+(dinn]-yo)*DX2I~i~n]); 1484 1485 DY2IT1=DY2IT1+Dyi[n]*Calpi[n]*DY2I2iEnl+(ci[nJ-xo)*DY2I~ icn]) 1486 -Dxi[nJ*Cgi~nJ*DX2I2i(n]+Cdinn]-yo)*DX2Ili[nJ-2.*DYIlin]); 1487 1488 1489 DXYIT1=DXYIT1+Dyi[n]*Cal i[n]*DXYI2i [nJ+(ci[nl-xo)*DXYIli[nJ-DYIli[nJ) 1490 -Dxi[nJ*Cgi n]*DXYI2innJ+(di(n]-yo)*DXYIli[n]-DXIlinnJ); 1491 } 1492 1493 1494 1* 1495 /* Calculate derivatives of IT2 1496 1* DX2IT2 is second derivative of 1T2 with respect to x 1497 /* DY2IT2 is second derivative of 1T2 with respect to y 1498 /* DXYIT2 is derivative of 1T2 with respect to x and y 1499 1600 1501 1502 DX2IT2=0.0; 1503 DY2IT2=0.0; 1604 DXYIT2=0.0; 1606 1506 for(n=l;n<=3;n++) 1507 { 1608 DX2IT2=DX2iTr2+Dyi [n] * (alpi [n]*al pi [n] *DX2I31 (n] +alpi (n]*(2. *ci (n] -xo) 1609 *DX2I2i[n]-.2.*alpi n]*DXI2iCn]+cinl]*(ci[n]-xo)*DX2Iiln] -2. 1510 *ci[n]*DXIli[n]); 1611 DX2IT2=DX2IT2-Dxi[n]*(alpi[n]l*gin]*DX2I3i[n]+C i[n]*cirn] 1512 +alpi[n]*(dirn]-yo))*DX2I2i[n]+ci[n]*('di nj-yo)*DX2I~i~n];

Jun 5 11:27 1986 Page 28 1513 1514 DY2IT2=DY2IT2+Dyi~n]*(alpi~n]*alpi~n]*DY2I3i En]+alpi~n]*C2.*ci~n]-xo) 1515 *DY2I2i~n]+ci~n],*(ci[n]-xo)*DY2Ili~n]); 1516 DY2IT2=DY2IT2-Dxi En] *(gi En]*alpi En] 1517 DYIin-*ain]DIinCgi[n] *ci (n] +alpi [n]*C(di[n] -yo)) 1518 *DY2I2i~n]-2.*ci n]*DYIli[n]+cinln*Cdi~n]-yo)*DY2I~iIn]); 1519 1520 DXYIT2=DXYIT2+Dyi En] *(alp E[n]*alpi En]*DXYI3i En]+al p i n] *(2. *ci En]-xo) 1521 *DXYI2i[n]-alpi~nJ *DYI2i~n]+ci~n]*Ccin] -xo) *DXYI1i[n]-ci~n] 1522 *DYliN]); 1523 DXYIT2=DXYIT2-Dxi~n]*(alpi~n]*gi~n]*DXYI3i~n]+Cgi En]*ci En]+alpi~n] 1524 * (di En] -yo) )*DXYI2i En] -alpiEn]*DXI-2i En]+ci En] *(di En]-yo) 1525 *DXYI1i En]-cI En]*DXI1iEn]) 1526} 1527 1528 DX2IT2=DX2IT2/2.0+xo*DX2IT1/2.0+DXIT1; 1529 DY2IT2=DY2IT2/2.0+xo*DY2IT1/2.0; 1530 DXYIT2=DXYIT2/2.0+DYIT1/2.0+xo*DXYIT1/2.0; 1531 1533/** 1534 /* Calculate derivatives of 1T3 * 1635 /* DX2IT3 is second derivative of 1T3 with respect to x * 15368/ DY2IT3 is second derivative of 1T3 with respect to y 1537 /* DXYIT3 is derivative of 1T3 with respect to x and y 1538 / 1539 /*************************************** 1540 1541 DX2IT3=0.0; 1542 DY2IT3=0.0; 1543 DXYIT3=0.O; 1544 1545 for~n~l;n<=3;n++) 1546{ 1547 DX2IT3=DX2IT3-~Dyi~n]*Calpi~n]*gi~n]*DX2I3i~n]+(gi En]*(ci n]-xo)+di En]) 1548 *DX2I2i(n]-2.*gi~nJ*DXI i~n]+di~n]*Cci~nJ -xo)*DX2Ili[n]-2.*di~n] 1549 *Ili~n]); 1550 DX2IT3=DX2IT3-DxI[n]*Cgi~n]*gi~n]*DX2I3i[n]+gi~n]*C2.*di~nJ-yo) 1551 *DX2I2i~n]4-di~nJ*Cdi~n]-yo)*DX2I~i En]J 1552 1553 DY2IT3=DY2IT3+Dyi~n]*Cgi[n]*alpi~n]*DY2I3i1[n]+Cgi~n]*Ccl[n]-xo)+di En]) 1554 *DY2I2i[n]+di~nJ *(ci~n]-xo)*DY2Ili~n]); 1556 DY2IT3=DY2IT3-Dxi En] * i(gEn]*gi En] 1556 *DY2I3i~n]+gi~nJ *(2.*di~n]-yo)*DY2I2i~n]-2.*gI~n]*DYI2i~n]+dl~n] 1557 *Cdi~n]-yo)*DY2I~i~n] —2.*dl En]*DYIli~n]); 1558 1559 DXYIT3=DXYIT3+Dyi (n] *Calpi[n] *gl(n] *DXYI3I En] +Cgi En] *(ci En] -xo) +dl En]) 1560 *DXYI2i En]-gi En] *DYI2iEn] +di En] * ci En]-xo) *DXYIli En] -di En] 1561 *DYIII~n]); 1562 DXYIT3=DXYIT3-Dxi En]*Cgi~n]*gi~n]*DXYI3i[n]+gi [n]*(2.*dliEn]-yo) 1563 *DXYI2i~n]-giEn] *DXI2i~n]+di~n] *CdiEn] -yo)*DXYI li~n]-di~n] 1564 *DXI~i~n]); 1565 15668} 1567 1568 DX2IT3=DX2IT3/2.0+yo*DX2IT-1/2.0;

Jun 5 11:27 1985 Page 29 1569 DY2IT3=DY2IT3/2.0~yo*DY2ITI/2.0+DYITI; 1570 DXYIT3=DXYIT3/2.0+DXITI/2.0+yo*DXYIT1/2.0; 1571 1573 /* 1674 /* Calculate derivatives ofI T7 1575 /* DXIT7 is first derivative of 1T7 with respect to x 1576 /* DYIT7 is first, derivative of IT7 with respect to y 1577 /* DX2IT7 is second derivative of IT7 with respect to x 1578 /* DY2IT7 is second derivative of 1T7 with respect to y 1579 /* DXYIT7 is derivative of 1T7 with respect to x and y 1580 */ 1581 /*********************************************************** 1582 1583 DXIT7O.O0; 1684 DYIT7O.0; 1685 DX2IT7=0.0; 1586 DY2IT7=0.0; 1587 DXYIT7O.0; 1588 1589 forCn=1;n<=3;n++) 1590 ( 1591 DXIT7=DXIT7+Dyi n]*(alpi n]*DXIi n]J+(ci[n]-xo)*DXI4i n]-I4i[nJ) 1592 -Dxi EnJ*(gin] n*DXISi En]+(di nJ-yo)*DXI4i nJ); 1593 1694 DYIT7=DYIT7+DyiEn]*CalpiEn]*DYI5in)Ci[nl En]-xo)*DYI4iEn]) 1596 -Dxi [nJ*gi nJ *DYI5i nJ+(di[n]-yo)*DYI4iUn]-I4i nJ); 1596 1697 DX2IT7=DX2IT7+Dyi~nJ*Calpinn]*DX2Ii nJ+(CcinnJ-xo)*DX2I4i~nJ-.2 1598 *DXI4i~nJ)-DxinJ*CgiEn]*DX2I5iEUJ+Cdi~nJ-yo)*DX2I4i~n]); 1599 1600 DY2IT7=DY2IT7+Dyi~nl]*alpinJ*DY2I5i~nJ+Cci~n]-xo)*DY2I4i~nJ) 1601 -Dxi~nJ*(gi~n]*DX2I5i~nJ+Cdinn]-yo)*DX2I4i-n]-2.*DYI4iyn]); 1602 1603 1604 DXYIT7=DXYIT7+Dyi En] * Calpi En] *DXYI5i n] + Cci E[n]) -o) *DXYI4i En) -DYI4i En] 1605 -Dxin]l*Cgitn]*DXYI5iEn)+Cdinn]-yo)*DXYI4i]n]-DXI4i[n]); 1606 1607 1608 DXIT7=DXIT7/3.0; 1609 DYIT7=DYIT7/3.0; 1610 DX2IT7=DX2IT7/3.0; 1611 DY2IT7=DY2IT7/3.0; 1612 DXYIT7=DXYIT7/3.0; 1613 1614 /**************************************/ 1615 /* 1616 /* Calculate derivatives of 1T8 1617 /* DX2IT8 is socond derivative of IT8 with respect to x 1618 /* DY2IT8 is second derivative of ITB with respect to y 1619 DXYIT8 is derivative of IT8 with respect to x and y */ 1620 /* 1621 1622 1623 DX2IT8=0.0; 1624 DY2IT8=0.0;

9MISSIN

1626 dLL-.V 1627 for~n=1;n<=3;n++) 1628{ 1629 DX2IT8=DX2IT8+Dyi [n] *(alpi En] *alpi [n] *DX2I6i [n] +alpi [n] *(2. *ci En] -xo) 1630 *DX2I5i~n]-2.*alpi[n]*DXI~i~n]+cinln*Cci[n] —xo)*DX2I4i~n]-2. 1631 *ci En] *DXM~ En]); 1632 DX2IT8=DX2IT8-Dxi [n] *(alpi En] *gi [n] *DX2I6i En] +(gi [n] *ci En] 1633 +alpi~n]*Cdiln]-yo))*DX2I5i~n]+ci~n]*(di tn]-yo)*DX2I4i~nJ); 1634 1635 DY2IT8=DY2IT8+Dyi En] *Calpi En]I*alpi En] *DY2I6i (n] +alpi [n] *(2. *ci En] -xo) 1636 *DY2I5i~n]+ci~n]*(ci~n]-xo)*DY2I4i~n]); 1637 DY2IT8=DY2IT8-Dxi En] *(gi En]*alpi En] 1638 *DY2IMin]-2.*a1pi~n]*DYI5i n]+(Ki~n]*ci~n]+alpi[nJ*Cdi~n]-yo)) 1639 *DY2I~i~n]-.2.*ci n]*DYIMin]+ci n]*Cdi~n]-yo)*DY2I4i~n]); 1640 1641 DXYIT8=DXYIT8+Dyi~n]*(alpi En]*alpi En] *DXYI6i En]+alp i~n]*(2.*ciEn) -xo) 1642 *DXYI5i En]-alpiEn] *DYI~i En]ci-clE]* Cci En]-xo) *DXYI4i1En] -ci nIn 1643 *DYI4i~nJ); 1644 DXYIT8=DXYIT8-Dxi~nJ*Calpi~n]*gi~nJ*DXYI6i En]-Cgi EnJ*ci~n]+alpi~n] 1645 *(di n]-yo))*DXYI5i~nJ-alpi~n]*DXI~i~n]+ci~n]*Cdi~nJ-yo) 1646 *DXYI4i~n]-ci~n]*DXI4iEn]) 1647} 1648 1649 DX2IT8=DX2IT8/4.0-xo*DX2IT7/4.0+DXIT1/2.0; 1660 DY2IT8=DY2IT8/4.0+xo*DY2IT7/4.0; 1651 DXYIT8=DXYIT8/4.0OiDYIT7/4. 0xo*DXYITI/4.0; 1652 1653 1654 /*************************************** 16556/ 16568/ Calculate derivatives of 1T9 * 1657 1* DX2IT9 is second derivative of 1T9 with respect to x 1658 /* DY2IT9 Is second derivative of 1T9 with respect to y * 1659 /* DXYIT9 is derivative of 1T9 with respect to x and y 1660 / 1661 /*************************************** 1662 1663 DX2IT9=0.0; 1664 DY2IT9=0.0; 1665 DXYIT9=0.0; 1666 1667 for (n=1;n<=3;n+i-) 1668{ 1669 DX2IT9=DX2IT9+Dyi (n]*(alpi (n] *gi (n] *DX2I~i (n] Cgi En]*(ci [n] -xo) +di ND) 1670 *DX2I~i (n] -2. *gi [n] *DXI5I [n] +di[n] *Cci (nJ -xo) *DX2I4i (n] -2. *di (n] 1671 *IMin]); 1672 DX2IT9=DX2IT9-DxiEn]*Cgi En]*gi En]*DX2I~i Enj+i-g n]*(2.*di En]-yo) 1673 *DX2I~i En]'diEn~ * (diEn] -yo) *DX2I.4i En]) 1674 1675 DY2IT9=DY2ITQ-sDyi (nJ *(gi[n] *alpi[n] *DY2I6i En) +C(gi En] *(cl En] -xo) +di EN]) 1676 *DY2I5i[n]+di~tnj *(cl~n]-xo)*DY2I4i~n]); 1677 DY2IT9=DY2IT9-DxiEn]*(g i En]*gi En] 1678 *DY2I61En]i-gt~nl*(2. *di~n]-yo)*DY2I~i~n]-2.*gi~n]*DYI5i~nli-di~nI 1679 *(di~n]-yo)*DY2I4i~n]-2.*d!irn] *DYI4i En]); 1680 1733 DY2IT9= (-1) *DY21IIv; 1734 DXYIT9=(-1) *DXYIT9; 17356 1736

Jun 5 11:27 1986 Page 32 1737 1738 1739 1740 1741 1742 1743 1744 1745 1746 1747 1748 1749 1760 1751 1762 1763 1764 1766 1756 1757 1758 1769 1760 1761 1762 1763 1764 1766 1766 1767 1768 1769 1770 1771 1772 1773 1774 1776 1776 1777 1778 1779 1780 1781 1782 1783 1784 1785 1786 1787 1788 1789 1790 1791 1792 powr(x,n) double x; int n; double p,pov(); if(x == 0.) { p=O.; return(p); else { p=pow(x,n); return(p); } I f numer(xo, yo xsi, ys 1,xs2,ys2,x3,s3) double xo,yo,xs1,ys1,xs2,ys2,xs3,ys3; { / */ /* This routine performs the numerical evaluations over the triangles */ xo x coordinate of the field point / yo y coordinate of the field point */ xsl x coordinate of vertice 1 of source triangle */ / ysl y coordinate of vertice 1 of source triangle xs2 x coordinate of vertice 2 of source triangle / ys2 y coordinate of vertice 2 of source triangle* xs3 x coordinate of vertice 3 of source triangle / ys3 y coordinate of vertice 3 of source triangle /* a weight factor / b weight factor */ c weight factor */ r evaluation point constant */ s evaluation point constant /1* ctrdx x coordinate of the centroid ctrdy y coordinate of the centroid xp[n] x coordinates of the evaluation points yp[n] y coordinates of the evaluation points area area of the triangle k wavenumbar rkr[n] real part of the rational kernel */ */,/

Jun 5 11:27 1985 Page 33 1793 /* rki[n] imag part of the rational kernel */ 1794 /* rkrx[n] real part of x times the rational kernel */ 1795 /* rkix[n] imag part of x times the rational kernal */ 1796 /* rkry[n] real part of y times the rational kernal */ 1797 /* rkiy[n] imag part of y times the rational kernel */ 1798 /* dx2rkr[n] real part of second derivative with respect to x of */ 1799 /* the rational kernel */ 1800 /* dx2rki(n] imag part of second derivative with respect to x of */ 1801 /* the rational kernel */ 1802 /* dx2rkrx(n] real part of x times second derivative with respect to x of*/ 1803 /* the rational kernel */ 1804 /* dx2rkix[n] imag part of x times second derivative with respect to x of*/ 1805 /* the rational kernel */ 1806 /* dx2rkry(n] real part of y times second derivative with respect to x of*/ 1807 /* the rational kernel */ 1808 /* dx2rkiy[n] imag part of y times second derivative with respect to x of*/ 1809 /* the rational kernel */ 1810 /* dy2rkr[n] real part of second derivative with respect to y of */ 1811 /* the rational kernel */ 1812 /* dy2rki(n] imag part of second derivative with respect to y of */ 1813 /* the rational kernel */ 1814 /* dy2rkrx[n] real part of x times second derivative with respect to y of*/ 1815 /* the rational kernel */ 1816 /* dy2rkix[n] imag part of x times second derivative with respect to y of*/ 1817 /* the rational kernel */ 1818 /* dy2rkry[n] real part of y times second derivative with respect to y of*/ 1819 /* the rational kernel */ 1820 /* dy2rkiy[n] imag part of y times second derivative with respect to y of */ 1821 /* the rational kernel */ 1822 /* dxyrkr[n] real part of derivative with respect to x and y of */ 1823 /* the rational kernel */ 1824 /* dxyrki[n] imag part of derivative with respect to x and y of */ 1825 /* the rational kernel */ 1826 /* dxyrkrx[n] real part of x times derivative with respect to x and y of */ 1827 /* the rational kernel */ 1828 /* dxyrkixtn] imag part of x times derivative with respect to x and y of */ 1829 /* the rational kernel */ 1830 /* dxyrkry(n] real part of y times derivative with respect to x and y of */ 1831 /* the rational kernel */ 1832 /* dxyrkiy[n] imag part of y times derivative with respect to x and y of */ 1833 /* the rational kernel */ 1834 /************************************************************************** 1835 /* */ 1836 /* EXTERNAL VARIABLES */ 1837 /* 1838 /**************************************************************************/ 1839 /* Rk re real part of the rational kernel */ 1840 /* Rk-im imag part of the rational kernel */ 1841 /* Rki re real part of x times the rational kernel */ 1842 /* Rkx-im imag part of x times the rational kernal */ 1843 /* Rkyfre real part of y times the rational kernal,/ 1844 /* Rky lm imag part of y times the rational kernel */ 1845 /* Dx2Fk_re real part of second derivative with respect to x of */ 1846 /* the rational kernel */ 1847 /* Dx2rk im imag part of second derivative with respect to x of */ 1848 /* the rational kernel */

Jun 6 11:27 1986 Page 34 1849 /* Dx2rkx re real part of x times second derivative with respect to x of */ 1850 /* the rational kernel */ 1851 /* Dx2rkxim imag part of x times second derivative with respect to x of */ 1852 * the rational kernel * 1853 /* Dx2rkyre real part of y times second derivative with respect to x of */ 1854 /* the rational kernel 1855 /* Dx2rkyim imag part of y times second derivative with respect to x of */ 1856 /* the rational kernel 1857 /* Dy2rk_re real part of second derivative with respect to y of */ 1858 /* the rational kernel */ 1859 /* Dy2rk_im imag part of second derivative with respect to y of */ 1860 /* the rational kernel */ 1861 /* Dy2rkxre real part of x times second derivative with respect to y of */ 1862 /* the rational kernel */ 1863 /* Dy2rkxim imag part of x times second derivative with respect to y of */ 1864 /* the rational kernel 1865 /* Dy2rky re real part of y times second derivative witn respect to y of */ 1866 /* the rational kernel */ 1867 /* Dy2rkyim imag part of y times second derivative with respect to y of */ 1868 /* the rational kernel */ 1869 /* Dxyrkre real part of derivative with respect to x and y of */ 1870 /* the rational kernel */ 1871 /* Dxyrkim imag part of derivative with respect to x and y of */ 1872 /* the rational kernel */ 1873 /* Dxyrkxre real part of x times derivative with respect to x and y of */ 1874 /* the rational kernel */ 1875 /* Dxyrkxim imag part of x times derivative with respect to x and y of */ 1876 /* the rational kernel */ 1877 /* Dxyrky re real part of y times derivative with respect to x and y of */ 1878 /* the rational kernel */ 1879 /* Dxyrkyim imag part of y times derivative with respect to x and y of */ 1880 /* the rational kernel */ 1881 /* dYk wavenumber; */ 1882 /**************************************************************************/ 1883 1884 double x[4],r([4],a.b,.cr.sk.ctrdx.ctrdy,xp[8],yp[8],area,R[8]; 1885 double rkr[8],rkik[8]rkrx[8],rkix[8] rkry[8] rkiy[8]; 1886 double dx2rkr[8],dx2rki[8],dx2rkrx[8],dx2rkix[8],dx2rkry[8],dx2rkiy[8]; 1887 double dy2rkr[8],dy2rki[8] dy2rkrx[8] dy2rkix[8],dy2rkry[8],dy2rkiy[8] 1888 double dxyrkr[8] dxyrki[8] dxyrkrx[8],dxyrkix[8],dxyrkry[8],dxyrkiy[8] 1889 double sqrt(),pow(); 1890 int n; 1891 1892 area=IT4; 1893 k=dYk; 1894 1896 x[l]=xsl; 1896 y[1]=ysl; 1897 x[2]=xs2; 1898 y[2]=ys2; 1899 x[3]=xs3; 1900 y[3]=ys3; 1901 1902 ************************************************************************** 1903 /* */ 1904 /* Calculate the centroid of the triangle */

Jun 5 11:27 1985 Page 35 19056/ 1907 1908 ctrdx —r( l]+x [2]+x [3]) /3.o0; 1909 ctrdy=C(y 1l]+y (2j+y(3]) /3. 0; 1910 1911/**************1**********************/ 1912 / 1913 /* If the field point and the centroid are coincident use the numerical*/ 1914 / method described in Hammer for a quadratic function to avoid the * 1915 /* singularity problem 1916 / 1917 /************************************* 1918 1919 if( Cxo == ctrdx US yo ==ctrdy) 1920 { 1921 r=0.5; 1922 1923 for Cn1l;n<=3;n++) 1924{ 1925 xp[n] =r*x [n] +C.0-r) *ctrdx; 1926 ~ p n]=r*y (n] +C (i0-r) *ctrdy; 1928 for (n=1;n<=3;n++) 1929 R (n] =sqrt (Cro-xp(nJ)*Cxo-xp (n])+C(yo-yp (n])*Cyo-yp EnJ)) 1930 1931 for(n=1;n<=3;n++) 1932 { 1933 1934 rkr~n]=cos~k*R[n])/R~n]-1.0/R~n]+Ic*k*R~n]/2.0; 1935 rki En]=sin Ck*R En] )/R En]-k; 1936 rkrx En] [ n] *(cos(k*REn]) /R [n]-i. /R En]+k*k*R En]/2.0); 1937 rkix n]=xpn]*Csin~k*R~n])/R~n]-k); 1938 rkry En]=yp En] *Ccos Ck*REn]) /REnJ -1. 0/R En]+k*k*R En)/2.0); 1939 rciy En]=yp En] *(sin Ck*REnl))/REn] -10; 1940 1941 1942 dx2rkr~n]=C3. *Cxo-xp En])*Cxo- En] )/CR En] *R~n]) -1.0-k*k 1943 *C(xo-xp[n]I) *Cxo- En]) *cos~k*R~n])/powCR~n].3) 1944 -CI.-3.*Cxo-xp n] *(xo-xp~n])/CR~n]*R~n]))*kc 1945 *sin (k*R (n])/CR n] *R n]) 1946 1947 dx2rki En]=C(I. 0-3. * Cxo-xp En] )*(xo-p En])/CR En] *R~n]) )*k*cos Ck*R En]) 1948 / (R(n] *R n] ) +(3. *Cxo-x En])*(xo-xp En])/CR~n]*R~n])-l.0 1949 -k*kc*Cxo-xp~n])*Cxo-xpn] ))*sin~k*R~n])/powCR~n].3); 1950 1951 dx2rkcrx~nl~xp~n]*C(3.*Cxo-xp~n])*(xo-xp~n])/CR~n]*R~n])-1l.0-k*k 1952 * Cxo-xp En])* Cxo-xp~n]))*cos~lk*R~n])/powCR~n]3) 1953 -C1.-3.*Cxo-xp~n] )*(xo-p En])/CR,[n]*t&n]))k 1954 *sin Ck*R En])/ CREn]*REn])) 1955 1956 dx2rkix~n]=xp En]*CC1.0-3.*Cxo-xp~n])*Cxo-xp~n])/CR~n]*R~n)))*'k*cos~k 1957 *R~n])/ CR~n]*R~n] )+C3.*Cxo-xp~n])*Cxo-xp~n])/CR[n]*R~n])-1.0 1958 -k*k*Cxo-xp~n])*Cxo-xp~n]))*sin~k*R~n])/powCR~n],3)); 1950 1960 dx2rkry En]=ypEn] * C(3.* Cxo-xp En])* Cxo-xpEn]) / CREn] *REn]) -1. 0-k*k

Jun 5 11:27 1985 Page 36 1961 *(xo-xpcn])*Cxo-xp [n]))*cos~k*R[n])/pow(R~nJ.3) 1962 -C1.-3.*(xo —xp[nJ)*(xo-xp En])/CR~n]*R[n]))*k 1963 *sin~k*R~n])/(Rtn]*R~n])) 1964 1965 dx2rkly[n]=y [n]*CC1.0-3.*(xo-xp[n])*(xo-xp[n])/(R~n]*R[n]))*k*cos(k 1966 *RI[h )/CR~n]*R~n])+(3.*(xo-xp[n])*Cxo-xp[n])/(R[n]*R[n])-l.0 1967 -k*k* Cia-xp n]) *(xo-xpEn])) *sin (k*R En] )/pow(R [n],3)); 19C8 1969 19~70 1971 dy2rkr~n]=(3.*Cyo-y En])*(yo-yp En])/(REn] *R En])-l.0-k*k* Cyo-p En]) 1972 *(yo-ypn] ))*cos Ck*R[n] )/pow(R[n].3)-(1.-3.*(yo-yp [nJ 1973 Cyo-y In])/(R~n]*R[n]))*k*sin~k*R~n])/(R[n]*R n]) 1974 1975 dy2rkl~n]=(1.0-3 *(yo-yp~n])*Cyo-y [n])/CR~n]*R[n]))*k*cos~k*R~n]) 1976 /(R~nJ*Rtn])+(3. *(yo-yp EnJ)*yo-yp En])/CR~n]*R[n] )-1.0 1977 -k*k*(yo-yp~n])*Cyo-ypnln]))*sin(k*R~n])/powCR~n].3); 1978 1979 dy2rkrx~n]=xp En]*((3.*(yo-yp~nj)*Cyo-ypn])/(R~n]*R[n])-l.0-c*Ik*Cyo1980 yp[n])*(yo-yp~n] ))*cos~k*R n])/powCR~n].3) -C1.-3.*(yo-yp[n])* 1981 1y-pn)(~]Rn)*~i(*~]/Rn*~]) 1982 (oy~j)C~]Rn)*~i~*~]/Rn*~]) 1983 dy2rkix~n]=x En]*CCI.0-3.*Cyo-ypn])*Cyo-yp-Cn])/CR~n]*R~n]))*kc*cos~k 1984 *R[n])/CR~n]*R~n])+C3. (o-yp~n])*(yo-y En])/CR~n]*R~n])-l.0 1985 -k*Ic*Cyo-yp~n])*(yo-yp n] ))*sin~k*R n])/powCR~n].3)); 1986 1987 dy2rkry~n]=ypEn) *CC3.*Cyo-yp~n])*Cyo-y En])/(R~n]*R~n])-1.0-kc*k*Cyo1988 ypE~] )*( o-p(n) ))*cos Ck*R n])/pow(R~n].3)-(l.-3 *Cyo-ypln])* 1989 (yo-yp~n))A En] *R En]))*k*sin Ck*R En])/(R En]*R En))); 1990 1991 dy2rkiy~n]=yp [n]*((1.0-3.*Cyo-'yp~n])*Cyo-yp~n])/CR~n]*R~n]))*k*cos(k 1992 *R[n])/CR~n]*R~n])+C3.*(yo-yp~n])*(yo-y En])/CR[n]*R~n])-1.0 1993 -Ic*k*(yo-yp~n])*(yo-ypEn) ))*sln~k*R Ln])/powCR~n].3)); 1994 1995 1996 1997 1998 dxyrkr~n]=CC3.-k*k*R~n]*R~n])*cos(k*R[n])+3.*kc*R~n]*sln~k*R~n])) 1999 *CXo-'xp~n])*Cyo-yp~n])/pow(R~n].5); 2000 2001 dxyrki~n]=C(3.-k*k*R~n]*R~n])*sin~k*R~nJ)-3.*k*R~n]*cos~k*R~n])) 2002 (ox n y-p(])/o R(),5 2003 *x-pn)Coy~]/oC~].) 2004 dxyrkrx~n]=x En]*C(C3.-k*k*R~n]*R~n])*cos~lc*R~n])+3.*k*R~n]*sinC~c 2005 *R(n ))*(xo-xp~n])*(yo-yp En])/pow(R En] 5)); 2006 2007 dxyrkixfn]=x En]*CCC3.-k*k*R~n]*R~n])*sin~k*R~n])-3.*k*R~n]*cos~k 2008 *R[n])) * (o-xp En])*(yo-yp En])/pow(REn].5)); 2009 2010 dxyrkry~n]=y En]*CCC3.-kc*k*R~n]*R~n])*cos~k*R~nJ)+3.*k*R~n]*sin~k 2011 *RCn ))*Cxo-xp~n])*(yo-yp~n])/pow(R[n].5)); 2012 2013 dxyrkiy~n]=yp En]*CCC3.-k*k*R~n]*R~n])*sin(k*R~n])-3.*k*R~n]*cos(k 2014 *R[-n]))*(xo-xp~n])*Cyo-yp~n])/pow(R[nI.5)); 2015 2016 }

Jun 5 11:27 1985 Page 37 2017 Rkc-re=(rkr[1]+rkr[2]1-rkr[3J)*area/3.0; 2018 2019 Rk-im=(rki[(1]+rki[(2ji-rki [3])*area/3.0; 2020 2021 Rkx re=Crkrx[1]+rkrx[2]+rkrx[3])*area/3.0; 2022 2023 Rkx im=Crkix[l]+rkix[2]+rkix[3])*area/3.0; 2024 -2025 Rky~re=Crkry[l]+rkry[2]+rkry[3] )*area/3.0; 2026 2027 Rky-im=(rkiy[1]+rkiy,[2]+rkiy[3])*area/3.0; 2028 2029 2030 Dx2rk-re=Cdx2rkr[1] +dx2rkr[2] +dx2rkr[3] )*area/3.0; 2031 2032 Dx2rk-im=Cdx2rki [1]+dx2rki [2]+dx2rki [3])*area/3.0; 2033_ 2034 Dx2rkcx re=Cdx2rlkrx[1]+dx2rkrx [2] +dx2rkrx [3]) *area/3.0; 2035 2038 Dx2rkx-im=Cdx2rkcix[l] 'dx2rkix[2] +dx2rkix[3])*area/3.O; 2037 2038 Dx2rkyVre=(dx2rlkry[1] +dx2rkry [2) +dx2rkry [3]) *area/3.0; 2039 2040 Dx2rky~im=Cdx2rkiy[1] idx2rkiy[2] i-dx2rkly [3]) *area/3.0; 2041 2042 2043 Dy2rk-re=(dy2rkr [1]+dy2rkr [2] +dy2rkr [3]) *area/3.0; 2044 2045 Dy2rk-im=Cdy2rki [1]+dy2rki [2) +dy2rki [3])*area/3.0; 2046 2047 Dy2rkx-re=(dy2rkrx[1] +dy2rkrx [2] +dy2rkrx [3]) *area/3.0; 2048 2049 Dy2rkx-im=Cdy2rkix[1] +dy2rlkix[2] +dy2rkix[3])*area/3.0; 2050 2051 Dy2rky-re=Cdy2rlkry[1] +dy2rlkry[2]+dy2rkry [3] )*area/3.0; 2052 2053 Dy2rky~im=Cdy2rkiy[1]+dy2rkiy[2]+dy2rkiy,[3])*area/3.0; 2054 2055 2056 Dxyrk-re=Cdxyrkr[1] +dxyrkr [2] +dxyrkr [3] )*area/3.0; 2057 2058 Dxyrk-im=Cdxyrki [1]+dxyrki [2]+dxyrki [3])*area/3.0; 2059 2060 Dxyrkx-re=Cdxyrkrx[1] +dxyrkcrx[2] +dxyrkrx[3] )*area/3.0; 2061 2062 Dxy~rkxim=(dxyrkix[l]+dxyrkix[2]+dxy-rklx[3])*area/3.0; 2083 2064 Dxyrky~re=(dxyrkry[1] +dxyrlkry[2] +dxyr-kry[3] )*area/3.0; 2065 2083 Dxyrky~ir=(dxyrkiy[1]+dxy~rkiy[2]+dxyrkiy[3])*area/3.0; 2067 2068 } 2069 2070 /************************************* 2071 / 2072 /* I tefeld point and the centroid are not coincident use the

Jun 6 11:27 1985 Page 38 2073 /* numerical method described in Hammer for a quintic function */ 2074 /* * 2076 /**************************************************************************/* 2076 2077 else 2078 { 2079 r=(1.0+sqrt(15.0))/7.0; 2080 s=(1.0-sqrt(15.0))/7.0; 2081 2082 a=(165.-sqrt(15.O))*area/1200.; 2083 b=(165.+sqrt(16.0))*area/1200.; 2084 c=9.0*area/40.0; 2085 2086 xp[l]=ctrdx;. 2087 yp[1]=ctrdy; 2088 2089 for(n=l;n<=3;n++) 2090 { 2091 xp[n+1]=r*x n]+(1.0O-r)*ctrdx; 2092 yp[n+1]=r*y [n] +(1.0-r)*ctrdy; 2093 xp[n+4]=s*x[n] + (1.0-s)*ctrdx; 2094 yp[n+4]=s*y in] + (1.0-s)*ctrdy; 2096 5 2096 2097 for(n=l;n<=77;n++) 2098 R[n]=sqrt((xo-xp[n])*(xo-xp[n])+(yo-yp En])*(yo-yp[n])); 2099 2100 for(n=i;n<=7;n++) 2101 { ' 2102 2103 2104 rkr[n]=cos(k*RIn])/R[n]-1.o/R[n]-+k*kc*R[n]/2.0; 2106 rki [n]=sin (k*R [n]) /R [n] -k; 2106 rkrx [n] =xp [n]*(cos(k*R[n])/R[n]-1.0/R[n]+k*k*R[n]/2.0); 2107 rkix[n] =xpn] * (sin (k*R [n]) /R n]-k); 2108 rkry[n]=yp[n]*(cos (k*R [n]) /R [n] -1.0/R[n] +k*k*R[n]/2.0); 2109 rkiy[n]=yp[n]*(sin(k*R[n])/R[n]-k); 2110 2111 2112 dx2rkr[n]= (3.*(xo-xp[n])*(xo-xp [n])/(R[n]*R[n] ) -1. O-*k 2113 *(xo-xp[n])*(xo-xp [n]))*cos(k*R[n])/pow(R[n],3) 2114 - (1.-3.*(xo-xp[n] )* (xo-xp [n])/(R[n]*R[n]))*k 2116 *sin (k*R [n] ) / (R [n]*R [n]); 2116 2117 dx2rki[n]=(1.0-3.*(xo-xp[n])*(xo-xp [n])/(R[n]*R[n]))*k*cos(k*R[n]) 2118 /(R[n]*R[n])+(3.*(xo-xp[n])*(xo-xp[n])/(R[n]*R[n])-1.0 2119 -k*k* (xo-xp [n]) * (xo-xpln]))*sin (k*R [n] )/pow(R [n],.3); 2120 2121 dx2rkrx[n]=xp[n]*((3.*(xo-xp[n])*(xo-xp [n])/(R[n]Rn]*Rin])-1.0-k*k 2122 *(xo-xp n])*(xo-xp [n]))*cos(k* [n])/pow(R[n].3) 2123 - (1.-3.*(xo-xp [n] )*(xo-xp [n])/(R[n] *[n]))*k 2124 *sin(k*R[n])/(R[n]*R[n])); 2125 2126 dx2rkix[n]=xp[n]*((1.0-3.*(xo-xp[n])*(xo-xp[n])/(R[n]*R[n]))*k*cos(k 2127 *R[n])/(R[n]*R[n])+(3.*(xo-xp[n])*(xo-xp[n])/(R[n]*R[ni])-1.0 2128 -k*k* (x o-xp n])*(xo-xpDn]))*sin(k*R n])/pow(R n],3));

Jun 5 11:27 1985 Page 39 2129 2130 dx2rkry~n]=yp~n]*((3.*Cxo-xp~n])*Cxo —xp~n])/CR~n]*R[n])-l.0-k*k 2131 *(xo-xptn])*(xo-xp [n]))*cos~k*R[n])/pow(R[n].3) 2132 -C1.-3.*(xo-xp[n] )*(xo-xp En])/CR~n]*R[n]))*k 2133 *sin(k*R~n])/(R~n]*R~n])) 2134 2135 dx2rkly~n)=yp [n]*((1.0-3.*Cxo-xp[n])*(xo-xp~n])/(R~n]*R~n]))*k*cos(k 2136 *R~n])/(R[n]*R[n])+(3.*(xo-xp[n])*(xo-.xp[,n])/(R[n]*R[n])-1.0 2137 -k*k*(xo-xp~nD)*(xo-xpnln))*sin(k*R~n])/powCR~n].3)); 2138 2139 2140 2141 dy2rkr~n]=C3.*(yo-y [n)*(yo-yp En])/(R~n] *R~n])-1.0-k*k*(yo-r En]) 2142 *Cyo-yp n]) *cos (k*R~n] )/pow(R~n].3P-C1.-3.*Cpo- p [nl 2143 (yo-yVpnln)/(R~n]*R~nJ))*k*sin(k*R[nJ)/CR~n] *R In); 2144 2145 dy2rkl~nJ=C1.0-3.*Cyo-yp~n])*(yo-y (nJ)/CR~n]*R~n]))*k*cos~k*R~nJ) 2146 /CR~n]*Rtn E]+C3.*Cyo-ypEn] )*Cyo-yp En])/(R[nJ*R[nJ)-1.O 2147 -Ic*k*Cyo-yp~n])*Cyo-ypfn] ))*sin~k*R~n])/pow(R~n].3); 2148 2149 dy2rkrx~n]~xp (nJ*((3.*(yo-yp~nJ)*Cyo-ypnJ)/CR~n]*R~nJ)-1.0-k*k*Cyo2150 y-p In )* ( yoypn] ))*cos~k*RtnJ)/powC(R[n].3) -CI. -3.(yo-yp (n]) 2151 (yo-yp[n) )/k nJ*R~nJ))*k*sin(k*R~n])/CR~n]*R~n])) 2152 2153 dy2rkix~n]=x (nJ*((1.0-3.*Cyo-yp~nJ)*Cyo.-yp[nJ)/CR En] *R~n]))*k*cos (k 2154 *R[n])/(R~nJ*R[nJ)+C3.*Cyo-yp~n])*CyO-y En])/CR~n]*R~n])-1.0 2155 -Ic*k*Cyo-y-p~n])*Cyo-ypEn] ))*sin~k*R ljn)powCR~n].3)); 2158 2157 dy2rkry[n]=yp [n]*C(3.*(yo-yp~nJ)*Cyoyp n])/CR~n]*R~n])-1.0-k*k*Cyo-' 2158 yp[nJ )*(yo-p(n) ))*cosC~k*R n)pow(R[n].3)-(1.-3.*Cyo-yp~nJ)* 2159 (yo-yp(n]) (fEn]*R~nJ))*k*sin(k*R~n])/CR[nJ*R[nJ)); 2160 2161 dy2rkly[nJly [nJ*CCI.0-3.*(yo-ypn])*Cyo-yp~n])/CR~n]*R~n]))*k*cos~k 2162 *RCn )/(R~nJ*R[nJ)+C3.*(yo-yp~nj)*Cyo-y EnJ)/CR~n]*R~nJ)-l.0 2163 *-k*k* (yo-y-p n])*(yo —yp[n) ))*sin~k*Rfn])/pow(R~n].3)); 2164 2165 2166 2167 2168 dxyrkr [n]C=C3. -k*k*R (nJ*R [n] )*cos (k*R(n] )+3. *k*R(nJ*sin (k*R (n]) 2169 (ox n y-p[])/o R[],5 2170 *x-pE]*y-pn)pwRE].) 2171 dxyrki [nJ=CC3. -k*k*R (nJ*R (n] )*sin (k*R (n ) -3. *k*R (nJ *cos (k*R (n]) 2172)(ox n y-p[])/o R(],5 2173 *x-pn)Coy~]/oC~J.) 2174 dxyrkrx [n]r —p [n]*C(3. -k*k*R(n] *R (nJ )*cos (k*R (n] )+3. *1*R [n]*sin (k 2175 *R~n) ))*(xo-xp~nJ)*Cyoyp~n])/powCR~n],5)); 2176 2177 dxyrklx (n] — p [n]*CC3. -k*k*R(n] *R (n]*sin (1*R (n] )-3. *k*R [n] *cos (k 2178 *R~n])) *(xo-xp En])*(yo-yp~n])/pow(R~nIs)); 2179 2180 dxyrkry[nl-yp En]*CCC3.-k*k*R~n]*R~n])*cos~k*R~n])+3.*k*R~n]*sln~k 2181 *R~n] ))*(xo-xp~n])*(yo-yp~n])/pow(R~n].5)); 2182 2183 dxyrkly~n]:y [n]*(CC3.-k~*k*R~n]*R~n])*sin(k*R~n])-3.*k*R~n]*cos(k 2O.1 84 *RCn]))*(xo-xp~n])*(yo-yp~n])/pow(CR~nJ 5));

Jun 5 11:27 1986 Page 40 2185 2186 } 2187 2188 Rkc re=c*rlcr[1] +a*Crkr[2] +rkr [3] +rkr[4] )+b*(rkr[5] +rkr [6] +rkr [7]); 2189 Rk lzn=c*rkil[1+a*Crki[2]+rki[3]+rki[4])+b*Crki[5]+rkl[6]+rki[7]); 2190 Rki re=c*rkrx[1]+a*(rkrx[2].,+rkrx[3]+rkrx[4])+b*(rkrx[5]+rkrx[6]+rkrx[7]); 2191 Rkxim~c*rkix[l]+a*(rkix[2]+rkix[3]+rkix[4])+b*(rkix[5]+~rkix[6]+rkix[7]); 2192 Rkf-re=c*rkry[1]+a*(rkry[2]+rkry[3]+rkry[l4])+b*(rkry,[5]+rkry[6] +rkry[7]); 2193 2194 2195 Dx2rk-re=c*dx2rkr[l] ia*(dx2rkr[2] +dx2rkr[3] +dx2rlcr[4] )+b 2198 *(dx2rkr[5]+dx2rkr[6] +dx2rkr[7]); 2197 2198 Dx2rk-im=c*dx2rki [l]+a*(dx2rki [2] +dx2rki [3] +dx2rki [4])+b 2199 *Cdx2rki615]+dx2rkl [6] +dx2rkl [7]); 2200 2201 Dx2rkx-re=c*dx2rkrx[1] +a*Cdx2rkrx[2] edx2rlkrx[3]+dx2rkrx[4] )+b 2202 *Cdx2rkrx[5]+dx2rkrx[6]+dx2rkrx[7]); 2203 2204 Dx2rkx im~c*dx2rkix[1] +a*(dx2rklx[2]+dx2rkix[3J+dx2rkix[4])+b 2205 *Cdx2rkix[5]+dx2rkcix[8]+dx2rklx[7J); 2206 2207 Dx2rky-re=c*dx2rkry [1]+a* Cdx2rkry [2]+dx2rkry [3] +dx2rkry [4]) +b 2208 *Cdx2rkryl 5] +dx2rkry [6] +dx2rkry [7]) 2209 2210 Dx2rky-im~c*dx2rki l y[]+a*Cdx2rkiy[2]+dx2rky [3] +dx2rkiy [4] )+b 2211 *Cdx2rkiy[5] +dx2rkiy[6]+dx2rkiy[7]) 2212 2213 2214 2215 Dy2rk re=c*dy2rkr [1]+a* (dy2rkr [2] +dy2rkr [3] +dy2rkr [4]) +b 2216 *Cdy2rkr[6]+dy2rkr[6]+dy2rkr[7]); 2217 2218 Dy2rkim~c*dy2rki [1]+a*Cdy2rki (2]+dy2rki [3]+dy2rki [4])+b 2219 *Cdy2rki[5]+dy2rlci[6]+dy2rki[7]); 2220 2221 Dy2rkx-re=c*dy2rkrx [1]+a* (dy2rkrx [2] +dy2rkrx [3] +dy2rkrx [4]) +b 2222 *Cdy2rkrx[5J+dy2rkrx[6]+dy2rkrx[7]); 2223 2224 Dy2rkx-im~c*dy2rkix[1] +a* Cdy2rkix[2] +dy2rklx[3] +dy2rkix[4] )+b 2225 *(dy2rkix[5]+dy2rkix[6]+dy2rkix[7]); 2226 2227 Dy2rky~re=c*dy2rkry [1] +a*(dy2rkry[2J +dy2rkry [3] +dy2rkry [4]) +b 2228 *(dy2rkry[5]+dy2rkry[6]+dy2rlkry[7]) 2229 22130 Dy2rky im=c*dy2rki y [1]+a* Cdy2rkiy[2] +dy2rky [3] +dy2rkly[4])+b 2231 *Cdy2rkiy[5] +dy2rk iy[6] +dy2rkiy [7] 2232 2233 2234 2235 DxyrkC re=c*dxyrkr li]+a* Cdxyrkr [2] +dxyrkr [3]+dxyrkr [4]) +b 2236 *Cdxyrcr [5] +dxyrkr [6] +dxyrkr [7]); 22137 2238 Dxyrk-im~c*dxyrki [1]+a*Cdxyrkl [2] +dxy-rkl[3] +dxyrki [4])+b 2239 * (dxyrkl [5] +dxyrki [6] +dxyrki [7]); 2240 Dxylrkx re~c*dxyrkrx[1]4a9*(dxyrkcrx [2]+dxyrkrx(3] +dxyrkrx[4])+b

Jun 5 11:27 1985 Page 41 2241 *(dxyrkrx [5] +dxyrkrx[6]+dxyrkrx [7]); 2242 2243 Dxyrkxim=c*dxyrkx [1 ] +a* (dxyrk [dxyr kxrkx [ 3+dxrkix dxyrkix [4]) +b 2244 *(dxyrklx[6]+dxyrkix[6]+dxyrklx[7]); 2245 2246 Dxyrky re=c*d+a*kry +a(dxyrkry [2] +dxyrkry [3]+dxyrkry[4]) +b 2247 *(dxyrkry 5]+dxyrkry[6]+dxyrkry[7]); 2248 2249 Dxyrkylm=c*dxyrky[+a*(dxyrk dxyrkly [2]+dxyrk [3] iy [4])+b 2260 *(dxyrkiy[] +dxyrkiy[6]+dxyrkiy [7] 2251 } 2252 } 2263 C NAASA 2.1.042 CGECO FTN-A 05-02-78 THE UNIV OF MICH COMP CTR 2264 SUBROUTINE CGECO(A,LDA,N,IPVt,RCONDZ) 2266 IMPLICIT REAL*8(A-H,O-Z) 2266 INTEGER LDA,N,IPVT(1) 2267 COMPLEX*16 A(LDA,1),Z(1) 2268 REAL*8 RCOND 2269 C 2260 C CGECO FACTORS A COMPLEX MATRIX BY GAUSSIAN ELIMINATION 2261 C AND ESTIMATES THE CONDITION OF THE MATRIX. 2262 C 2263 C IF RCOND IS NOT NEEDED, CGEFA IS SLIGHTLY FASTER. 2264 C TO SOLVE A*X = B, FOLLOW CGECO BY CGESL. 2266 C TO COMPUTE INVERSE(A)*C, FOLLOW CGECO BY CGESL. 2266 C TO COMPUTE DETERMINANT(A), FOLLOW CGECO BY CGEDI. 2267 C TO COMPUTE INVERSE(A), FOLLOW CGECO BY CGEDI. 2268 C 2269 C ON ENTRY 2270 C 2271 C A COMPLEX(LDA, N) 2272 C THE MATRIX TO BE FACTORED. 2273 C 2274 C LDA INTEGER 2276 C THE LEADING DIMENSION OF THE ARRAY A 2276 C 2277 C N INTEGER 2278 C THE ORDER OF THE MATRIX A 2279 C 2280 C ON RETURN 2281 C 2282 C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS 2283 C WHICH WERE USED TO OBTAIN IT. 2284 C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE 2285 C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER 2286 C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. 2287 C 2288 C IPVT INTEGER (N) 2289 C AN INTEGER VECTOR OF PIVOT INDICES. 2290 C 2291 C RCOND REAL 2292 C AN ESTIMATE OF THE RECIPROCAL CONDITION OF A 2293 C FOR THE SYSTEM A*X = B, RELATIVE PERTURBATIONS 2294 C IN A AND B OF SIZE EPSILON MAY CAUSE 2296 C RELATIVE PERTURBATIONS IN X OF SIZE EPSILON/RCOND 2296 C IF RCOND IS SO SMALL THAT THE LOGICAL EXPRESSION

Jun 6 11:27 1985 Page 42 2297 C 1.0 + RCOND.EQ. 1.0 2298 C IS TRUE, THEN A MAY BE SINGULAR TO WORKING 2299 C PRECISION. IN PARTICULAR, RCOND IS ZERO IF 2300 C EXACT SINGULARITY IS DETECTED OR THE ESTIMATE 2301 C UNDERFLOWS. 2302 C 2303 C Z COMPLEX (N) 2304 C A WORK VECTOR WHOSE CONTENTS ARE USUALLY UNIMPORTANT. 2305 C IF A IS CLOSE TO A SINGULAR MATRIX, THEN Z IS 2306 C AN APPROXIMATE NULL VECTOR IN THE SENSE THAT 2307 C NORM(A*Z) = RCOND*NORM(A)*NORM(Z) 2308 C 2309 C LINPACK. THIS VERSION DATED 07/14/77 2310 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. 2311 C 2312 C SUBROUTINES AND FUNCTIONS 2313 C 2314 C LINPACK CGEFA 2316 C BLAS CAXPY,CDOTC,CSSCAL,SCASUM 2316 C FORTRAN DABS,DIMAG,DMAX1,DCMPLX,DCONJG,REAL 2317 C 2318 C INTERNAL VARIABLES 2319 C 2320 IMPLICIT REAL*8(A-H,O-Z) 2321 COMPLEX*16 CDOTC,EK,T,WK,WKM 2322 REAL*8 ANORM,S,SCASUM,SM,YNORM 2323 INTEGER INFO.J,K,KB,KP1,L 2324 C 2325 COMPLEX*16 ZDUM,ZDUM1,ZDUM2,CSIGN1 2326 REAL*8 CABS1 2327 CABS1(ZDUM) = DABS(DREAL(ZDUM)) + DABS(DIMAG(ZDUM)) 2328 CSIGN1(ZDUM1,ZDUM2) = CABS1(ZDUM1)*(ZDUM2/CABS1(ZDUM2)) 2329 C 2330 C COMPUTE 1-NORM OF A 2331 C 2332 ANORM = O.ODO 2333 DO 10 J = 1, N 2334 ANORM = DMAX (ANORM,SCASUM(N,A(1,J),1)) 2335 10 CONTINUE 2336 C 2337 C FACTOR 2338 C 2339 CALL CGEFA(A,LDA,N,IPVT,INFO) 2340 C 2341 C RCOND = 1/(NORM(A)*(ESTIMATE OF NORM(INVERSE(A)))) 2342 C ESTIMATE = NORM(Z)/NORM(Y) WHERE A*Z = Y AND CTRANS(A)*Y = E 2343 C CTRANS(A) IS THE CONJUGATE TRANSPOSE OF A 2344 C THE COMPONENTS OF E ARE CHOSEN TO CAUSE MAXIMUM LOCAL 2345 C GROWTH IN THE ELEMENTS OF W WHERE CTRANS(U)*W = E 2346 C THE VECTORS ARE FREQUENTLY RESCALED TO AVOID OVERFLOW. 2347 C 2348 C SOLVE CTRANS(U)*W = E 2349 C 2350 EK = DCMPLX(1.ODO, O.ODO) 2351 DO 20 J = 1, N 2352 Z(J) = DCMPLX(0.ODO,O.ODO)

Jun 6 11:27 1985 Page 43 2353 20 CONTINUE 2354 DO 100 K = 1, N 2355 IF (CABS1(Z(K)).NE. 0.ODO) EK = CSIGN (EK,-Z(K)) 2356 IF (CABS1(EK-Z(K)).LE. CABS1(A(K,K))) GO TO 30 2357 S = CABS1(A(K,K))/CABS1(EK-Z(K)) 2358 CALL CSSCAL(N,S,Z,1) 2359 EK = DCMPLX(S,O.ODO)*EK 2360 30 CONTINUE 2361 WK = EK - Z(K) 2362 WKM = -EK - Z(K) 2363 S = CABS1(WK) 2364 SM = CABS1(WKM) 2365 IF (CABS1(A(K,K)).EQ. O.ODO) GO TO 40 2366 WK = WK/DCONJG(A(K,K)) 2367 WKM = WKM/DCONJG(A(K,K)) 2368 GO TO 60 2369 40 CONTINUE 2370 WK = DCMPLX(1.ODO,O.ODO) 2371 WKM = DCMPLX(1.ODO,O.ODO) 2372 50 CONTINUE 2373 KP1 = K + 1 2374 IF (KP1.GT. N) GO TO 90 2375 DO 60 J = KP1, N 2376 SM = SM + CABS1(Z(J)+WKM*DCONJG(A(K,J))) 2377 Z(J) = Z(J) + WK*DCONJG(A(K,J)) 2378 S = S + CABS1(Z(J)) 2379 60 CONTINUE 2380 IF (S.GE. SM) GO TO 80 2381 T = WKM - WK 2382 WK = WKM 2383 DO 70 J = KP1, N 2384 Z(J) = Z(J) + T*DCONJG(A(K,J)) 2385 70 CONTINUE 2386 80 CONTINUE 2387 90 CONTINUE 2388 Z(K) = WK 2389 100 CONTINUE 2390 S = 1.ODO/SCASUM(N,Z,1) 2391 CALL CSSCAL(N,S,Z,1) 2392 C 2393 C SOLVE CTRANS(L)*Y = V 2394 C 2395 DO 120 KB = 1, N 2396 K = N + 1 - KB 2397 IF (K.LT. N) Z(K) = Z(K) + CDOTC(N-K,A(K+1,K),1,Z(K+1),1) 2398 IF (CABS1(Z(K)).LE. 1.ODO) GO TO 110 2399 S = 1.ODO/CABSI(Z(K)) 2400 CALL CSSCAL(N,SZ,1) 2401 110 CONTINUE 2402 L = IPVT(K) 2403 T = Z(L) 2404 Z(L) = Z(K) 2405 Z(K) = T 2406 120 CONTINUE 2407 S = 1.0DO/SCASUM(N,Z,1) 2408 CALL CSSCAL(N,S,Z,1)

Jun 5 11:27 1985 Page 44 2409 2410 2411 2412 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 2439 2440 2441 2442 2443 2444 2445 2446 2447 2448 2449 2450 2451 2.4 52 2453 2454 24655 2456 2457 2458 2459 2460 2451 2462 2463 2464 C C C C C C C YNORM = 1.000 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 CAXPY(N-K,TDA(K+1,K),1,Z(K+1),1) IF CCABS1CZ(K)).LE. 1.000) GO TO 130 S = 1.ODO/CABS1(Z(K)) CALL CSSCAL(NDS,Z,1) YNORM = S*YNORM 130 CONTINUE 140 CONTINUE S = 1.ODO/SCASUMCNZ,1) CALL CSS'CALCN,SZ,1) YNORM = S*YNORM SOLVE U*Z = V DO 160 KB = D N K =N + 1 - KB IF CCABS1CZCK)).LE. CABS1CA(KK))) GO S = CABS1(A(KK))/CABS1CZCK)) CALL CSSCALCNSZD1) YNORM = S*YNORM 150 CONTINUE IF CCABS1(A(KK)) N&E. 0.000) ZCK) = IF CCABS1CACKDK)).EQ. 0.000) Z(K) = D T =-ZC(K) CALL CAXPYCK-1,T,AC1DK),1DZC1),l) 160 CONTINUE C MAKE ZNORM = 1.0 S = 1.ODO/SCASUM(NZ,1) CALL CSSCAL(N,S,Z,1) YNORM = S*YNORM C IF (ANORM.NE. 0.ODO) RCOND = YNORM/ANORM IF CANOM4.EQ. 0.000) RCOND = 0.ODO RETURIN END C NAASA 2.1.044 CGESL FTN-A 05-02-78 THE SUBROUTINE CGESL(A,LDAND IPVT.DD JOB) IMPLICIT REAL*8CA-H,O-Z) INTEGER LDAD N, IPVT (1), JOB COM?-PLEX* 16 A (LDA, 1),B (1) TO 150 (K) /A (K, K) CM.PLX (1. ODO, 0. ODO) EUNIV OF MICH COMP CTR 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

Jun 6 11:27 1985 Page 46 2465 C A COMPLEX(LDA, N) 2466 C THE OUTPUT FROM CGECO OR CGEFA. 2467 C 2468 C LDA INTEGER 2469 C THE LEADING DIMENSION OF THE ARRAY A 2470 C 2471 C N INTEGER 2472 C THE ORDER OF THE MATRIX A 2473 C 2474 C IPVT INTEGER (N) 2475 C THE PIVOT VECTOR FROM CGECO OR CGEFA. 2476 C 2477 C B COMPLEX (N) 2478 C THE RIGHT HAND SIDE VECTOR. 2479 C 2480 C JOB INTEGER 2481 C = 0 TO SOLVE A*X = B, 2482 C = NONZERO TO SOLVE. CTRANS(A)*X = B WHERE 2483 C CTRANS (A) IS THE CONJUGATE TRANSPOSE. 2484 C 2485 C ON RETURN 2486 C 2487 C B THE SOLUTION VECTOR X 2488 C 2489 C ERROR CONDITION 2490 C 2491 C A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A 2492 C ZERO ON THE DIAGONAL. TECHNICALLY THIS INDICATES SINGULARITY 2493 C BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER 2494 C SETTING OF LDA. IT WILL NOT OCCUR IF THE SUBROUTINES ARE 2495 C CALLED CORRECTLY AND IF CGECO HAS SET RCOND.GT. 0.0 2496 C OR CGEFA HAS SET INFO.EQ. 0 2497 C 2498 C TO COMPUTE INVERSE(A) * C WHERE C IS A MATRIX 2499 C WITH P COLUMNS 2500 C CALL CGECO (A,LDA,N, IPVT,RCOND,Z) 2501 C IF (RCOND IS TOO SMALL) GO TO... 2502 C DO 10 J = 1, P 2503 C CALL CGESL(A,LDA,N,IPVT,C(1,J),O) 2504 C 10 CONTINUE 2505 C 2506 C LINPACK. THIS VERSION DATED 07/14/77 2507 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. 2508 C 2509 C SUBROUTINES AND FUNCTIONS 2510 C 2511 C BLAS CAXPY,CDOTC 2512 C FORTRAN DCONJG 2513 C 2514 C INTERNAL VARIABLES 2515 C 2516 COMPLEX*16 CDOTC,T 2517 INTEGER K, KL,NM 1 2518 C 2519 NM.1 = N - 1 2520 IF (JOB.NE. 0) GO TO 50

Jun 5 11:27 1985 Page 46 2521 C 2522 C JOB = 0 SOLVE A * X = B 2523 C FIRST SOLVE L*Y = B 2524 C 2525 IF (NM1.LT. 1) GO TO 30 2525 DO 20 K = 1, NM1 2527 L = IPVT(K) 2528 T = B(L) 2529 IF (L.EQ. K) GO TO 10 2530 B(L) = B(K) 2531 B(K) = T 2532 10 CONTINUE 2533 CALL CAXPY(N-K,T,A(K+1,K),1,B(K+1),1) 2534 20 CONTINUE 2536 30 CONTINUE 2536 C 2537 C NOW SOLVE U*X = Y 2538 C 2539 DO 40 KB = 1, N 2540 K = N + 1 - KB 2541 B(K) = B(K)/A(K,K) 2542 T = -B(K) 2543 CALL CAXPY(K-1,T,A(1,K),1,B(1), 1) 2544 40 CONTINUE 2545 GO TO 100 2546 50 CONTINUE 2547 C 2548 C JOB = NONZERO, SOLVE CTRANS(A) * X = B 2549 C FIRST SOLVE CTRANS(U)*Y = B 2650 C 2551 DO 60 K = 1, N 2552 T = CDOTC(K-1,A(1,K),1,B(1),1) 2553 B(K) = (B(K) - T)/DCONJG(A(K,K)) 2554 60 CONTINUE 266555 C 2656 C NOW SOLVE CTRANS(L)*X = Y 2557 2558 2559 2560 2661 2662 26563 2664 2566 2566 2667 2568 2569 2570 2571 2572 2673 2574 2575 2576 C IF (NM1.LT. 1) GO TO 90 DO 80 KB = 1, NM1 K = N - KB B(K) = B(K) + CDOTC(N-K,A(K+1,K) L = IPVT(K) IF (L.EQ. K) GO TO 70 T = B(L) B(L) = B(K) B (K) = T 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE RETURN END C NAASA 2.1.043 CGEFA FTN-A 05-02-78 SUBROUTINE CGEFA(A, LDA,N,IPVT, INFO) IMPLICIT REAL*8 (A-H, O-Z) INTEGER LDA,N,IPVT(1), INFO,1,B(K+I),1) THE UNIV OF MICH COMP CTR

Jun 6 11:27 1985 Page 47 2577 COMPLEX*16 A(LDA, 1) 2578 C 2579 C CGEFA FACTORS A COMPLEX MATRIX BY GAUSSIAN ELIMINATION. 2580 C 2581 C CGEFA IS USUALLY CALLED BY CGECO, BUT IT CAN BE CALLED 2582 C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED. 2583 C (TIME FOR CGECO) = (1 + 9/N)*(TIME FOR CGEFA) 2584 C 2585 C ON ENTRY 2586 C 2587 C A COMPLEX(LDA, N) 2588 C THE MATRIX TO BE FACTORED. 2589 C 2590 C LDA INTEGER 2591 C THE LEADING DIMENSION OF THE ARRAY A. 2592 C 2593 C N INTEGER 2594 C THE ORDER OF THE MATRIX A. 2595 C 2596 C ON RETURN 2597 C 2598 C A AN UPPER TRIANGULAR MATRIX AND THE MULTIPLIERS 2599 C WHICH WERE USED TO OBTAIN IT. 2600 C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE 2601 C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER 2602 C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR. 2603 C 2604 C IPVT INTEGER (N) 2605 C AN INTEGER VECTOR OF PIVOT INDICES. 2606 C 2607 C INFO INTEGER 2608 C = 0 NORMAL VALUE. 2609 C = K IF U(K,K).EQ. 0.0. THIS IS NOT AN ERROR 2610 C CONDITION FOR THIS SUBROUTINE, BUT IT DOES 2611 C INDICATE THAT CGESL OR CGEDI WILL DIVIDE BY ZERO 2612 C IF CALLED. USE RCOND IN CGECO FOR A RELIABLE 2613 C INDICATION OF SINGULARITY. 2614 C 2615 C LINPACK. THIS VERSION DATED 07/14/77. 2616 C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LABS. 2617 C 2618 C SUBROUTINES AND FUNCTIONS 2619 C 2620 C BLAS CAXPY,CSCAL, ICAMAX 2621 C FORTRAN DABS,DIMAG,DCMPLX,DREAL 2622 C 2623 C INTERNAL VARIABLES 2624 C 2625 COMPLEX*16 T 2626 INTEGER ICAMAX,J,K,KP1,L,NM1 2627 C 2628 COMPLEX*16 ZDUM 2629 REAL*8 CABS1 2630 CABS1 (ZDUM) = DABS(DREAL(ZDUM)) + DABS(DIMAG(ZDUM)) 2631 C 2632 C GAUS3IAN ELIMINATION WITH PARTIAL PIVOTING

Jun 6 11:27 1985 Page 48 2633 2634 2636 2636 2637 2638 2639 2640 2641 2642 2643 2644 2645 2646 2647 2648 2649 2650 2651 2652 2653 2664 2665 2666 2657 2668 2669 2660 2661 2662 2663 2664 2666 2666 2667 2668 2669 2670 2671 2672 2673 2674 2675 2676 2677 2678 2679 2680 2681 2682 2683 2684 2686 2686 2687 2688 C C C C C C C C C C C C C C C C C C C C C INFO = 0 NM = N - 1 IF (NM1.LT. 1) GO TO 70 DO 60 K = 1, NM1 KP1 = K + 1 FIND L = PIVOT INDEX L = ICAMAX(N-K+1,A(K,K),1) + K - 1 IPVT(K) = L ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED IF (CABS1(A(L,K)).EQ. O.ODO) GO TO 40 INTERCHANGE IF NECESSARY IF (L.EQ. K) GO TO 10 T = A(L,K) A(L,K) = A(K,K) A(K,K) = T 10 CONTINUE COMPUTE MULTIPLIERS T = -DCMPLX(1.ODO,O.ODO)/A(K,K) 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+I,K),1,A(K+1,J),1) 30 CONTINUE GO TO 60 40 CONTINUE INFO = K 60 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,CX, INCX,CY, INCY) CONSTANT TIMES A VECTOR PLUS A VECTOR. JACK DONGARRA, LINPACK, 6/17/77. IMPLICIT REAL*8(A-H,O-Z)

Jun 5 11:27 1986 Page 49 2689 COMPLEX*16 CX(),CY(1),CA 2690 INTEGER I,INCX,INCY,IX,IY,N 2691 C 2692 IF(N.LE.O)RETURN 2693 IF (DABS(DREAL(CA)) + DABS(DIMAG(CA)).EQ. O.ODO ) RETURN 2694 IF(INCX.EQ.1.AND.INCY.EQ.1)GOTO 20 2695 C 2696 C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS 2697 C NOT EQUAL TO 1 2698 C 2699 IX = 1 2700 IY = 1 2701 IF(INCX.LT.O)IX = (-N+I)*INCX + 1 2702 IF(INCY.LT.O)IY = (-N+1)*INCY + 1 2703 DO 10 I = 1,N 2704 CY(IY) = CY(IY) + CA*CX(IX) 2705 IX = IX + INCX 2706 IY = IY + INCY 2707 10 CONTINUE 2708 RETURN 2709 C 2710 C CODE FOR BOTH INCREMENTS EQUAL TO 1 2711 C 2712 20 DO 30 I = 1,N 2713 CY(I) = CY(I) + CA*CX(I) 2714 30 CONTINUE 2716 RETURN 2716 END 2717 C NAASA 1.1.012 CDOTC FTN-A 06-02-78 THE UNIV OF MICH COMP CTR 2718 COMPLEX*16 FUNCTION CDOTC(N.CX,INCX,CY,INCY) 2719 C 2720 C FORMS THE DOT PRODUCT OF TWO VECTORS, CONJUGATING THE FIRST 2721 C VECTOR. 2722 C JACK DONGARRA, LINPACK, 6/17/77. 2723 C 2724 IMPLICIT REAL*8(A-H,O-Z) 2725 COMPLEX*16 CX(1),CY(1),CTEMP 2726 INTEGER I,INCX,INCY,IX,IY,N 2727 C 2728 CTEMP = DCMPLX (O0.ODO,O.ODO) 2729 CDOTC = DCMPLX(O.ODO,O.ODO) 2730 IF (N. LE. O)RETURN 2731 IF(INCX.EQ.1.AND.INCY.EQ. 1)GOTO 20 2732 C 2733 C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMFNTS 2734 C NOT EQUAL TO 1 2735 C 2736 IX = 1 2737 IY = 1 2738 IF(INCX.LT.O)IX = (-N+1)*INCX + 1 2739 IF(INCY.LT.O)IY = (-N+1)*INCY + 1 2740 DO 10 I = 1,N 2741 CTEMP = CTEMP + DCONJG(CX(IX))*CY(IY) 2742 IX = IX + INCX 2743 IY = IY + INCY 2744 10 CONTINUE

Jun 6 11:27 1985 Page 60 2745 2746 2747 2748 2749 2760 2751 2752 2763 2754 2755 2756 2757 2758 2769 2760 2761 2762 2763 2764 2765 2766 2767 2768 2769 2770 2771 2772 2773 2774 2776 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 CDOTC = CTEMP RETURN CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C 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 THE UNIV OF MICH COMP SUBROUTINE CSSCAL(N,SA,CX,INCX) C C SCALES A COMPLEX VECTOR BY A REAL CONSTANT. C JACK DONGARRA, LINPACK, 6/17/77. C IMPLICIT REAL*8(A-H, O-Z) COMPLEX*16 CX (1) REAL*8 SA 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) = DCMPLX(SA*DREAL(CX(I)),SA*DIMAG(CX(I))) 10 CONTINUE RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 DO 30 I = 1,N 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 REAL*8 FUNCTION SCASUM(N,CX,INCX) C C TAKES THE SUM OF THE ABSOLUTE VALUES OF A COMPLEX VECTOR AND C RETURNS A SINGLE PRECISION RESULT. C JACK DONGARRA, LINPACK, 6/17/77. C IMPLICIT REAL*8(A-H,O-Z) COMPiLEX*16 CX(1) REAL*3 STEMP INTEGER I,INCX,N,NINCX C SCASUM = O.ODO STE4P = O.ODO IF (N. LE. 0) RETURN IF(INCX.EQ.1)GOTO 20 CTR CTR

Jun 5 11:27 1985 Page 51 2801 2802 2803 2804 2805 2806 2807 2808 2809 2810 2811 2812 2813 2814 2816 2816 2817 2818 2819 2820 2821 2822 2823 2824 2825 2826 2827 2828 2829 2830 2831 2832 2833 2834 2835 2836 2837 2838 2839 2840 2841 2842 2843 2844 2845 2846 2847 2848 2849 2860 2861 2852 2853 2864 2855 2856 C C C CODE FOR INCREMENT NOT EQUAL TO 1 NINCX = N*INCX DO 10 I = 1,NINCX,INCX STEMP = STEMP + DABS(DREAL(CX(I))) + 10 CONTINUE SCASUM = STEM? RETURN C C CODE FOR INCREMENT EQUAL TO 1 C 20 DO 30 I = 1,N STEMP = STEMP + DABS (DREAL(CX (I))) + 30 CONTINUE SCASUM = STEMP RETURN END C NAASA 1.1.019 CSCAL FTN-A 05-02-78 SUBROUTINE CSCAL(N,CA,CX, INCX) C C SCALES A VECTOR BY A CONSTANT. C JACK DONGARRA, LINPACK, 6/17/77. C IMPLICIT REAL*8 (A-H, O-Z) COMPLEX*18 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 INTEGER FUNCTION ICAMAX(N,CX, INCX) C C FINDS THE INDEX OF ELEMENT HAVING MAX. 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 DABS (DIMAG(CX(I))) DABS (DIMAG (CX(I))) THE UNIV OF MICH COMP CTR THE UNIV OF MICH COMP CTR ABSOLUTE VALUE.

Jun 6 11:27 1985 Page 62 2857 COMPLEX*16 ZDUM 2858 REAL*8 CABS1 2859 CABS1(ZDUM) = DABS(DREAL(ZDUM)) + DABS(DIMAG(ZDUM)) 2860 C 2861 ICAMAX = 1 2862 IF(N.LE. 1)RETURN 2863 IF(INCX.EQ.1)GOTO 20 2864 C 2865 C CODE FOR INCREMENT NOT EQUAL TO 1 2866 C 2867 IX = 1 2868 SMAX = CABS1 (CX(1)) 2869 IX = IX + INCX 2870 DO 10 I = 2,N 2871 IF(CABSI(CX(IX)).LE.SMAX) GO TO 6 2872 ICAMAX = I 2873 SMAX = CABS1(CX(IX)) 2874 6 IX = IX + INCX 2876 10 CONTINUE 2876 RETURN 2877 C 2878 C CODE FOR INCREMENT EQUAL TO 1 2879 C 2880 20 SMAX = CABS (CX(1)) 2881 DO 30 I = 2,N 2882 IF(CABS1(CX(I)).LE.SMAX) GO TO 30 2883 ICAMAX = I 2884 SMAX = CABS1 (CX(I)) 2885 30 CONTINUE 2836 RETURN 2887 END 2888 CCCCCCCCCCCCCCCCCCCCCCCCC 2889 C DREAL DOESN'T SEEM TO WORK, SO THIS FUNCTION IS A SUBSTITUTE 2890 REAL*8 FUNCTION DREAL(X) 2891 COMPLEX*16 X,X2 2892 REAL*8 XA(2) 2893 EQUIVALENCE (X2,XA(1)) 2894 X2=X 2895 DREAL=XA(1) 2895 RETURN 2897 END