AFCRL-72-0162 013630-9-T 1363-9-T = RL-2030 THE UNIVERSITY OF MICHIGAN COLLEGE OF ENGNtRING DEPARTMENT OF ELECTRICAL AND COMPUTER ENGINEERING Radiation Laboratory THE NUMERICAL SOLUTION OF LOW FREQUENCY SCATTERING PROBLEMS By THOMAS B.A. SENIOR AND DAVID J. AHLGREN February 1972 Scientific Report No. 12 Contract No. F19628-68-C-0071 Project No. 5635 Task No. 563502 Work Unit No. 56350201 Approved for public release; distribution unlimited. r\ Contract Monitor: JOHN K. SCHINDLER ) MICROWAVE PHYSICS LABORATORY Prepared for: Air Force Cambridg Research Laboratories Air Force Systems Command Laurence G. Hanocom Field Bedford, Masachusetts 01730 Ann Arbor, Michigan 48105

Qualified requestors may obtain additional copies from the Defense Documlllentatiorl Center. All others should applY1 to the National I r. l iii ic.l Information Service.

013630-9-T ABSTRACT The low frequency scattering of electromagnetic and acoustic waves by rotationally symmetric bodies is considered. By concentrating on certain quantities such as the normalised component of the induced electric and magnetic dipole moments, it is shown how the first one or two terms in the far zone scattered fields can be expressed in terms of quantities which are functions only of the geometry of the body. Each of these is the xeighted integral of an elemientary potential function which can be found by solving an integral equation. A computer program has been written to solve the appropriate equations by the moment method, and for calculating the dipole moments, the electrostatic capacity, and a further quantity related to the capacity. The program is described and related data are presented. ii

01 3f 30-9-T TABLE OF CONTENTS ABSTRACT ii ACKNOWLEDGEMENTS iv I INTRODUCTION 1 II PERFECTLY CONDUCTING BODIES 3 2.1 Formnulation 3 2.2 Procedure for P1 9 2.3 Procedure for P33 and C 11 2.4 Procedure for M11 18 2.5 Procedure for M33 21 2.6 Disjoint Surfaces 22 II ACOUSTICALLY SOFT OR HARD BODIES 33 3.1 General Procedure 33 3.2 Soft Bodies 35 3.3 Hard Bodies 40 IV THE COMPUTATIONAL TASK 52 4.1 Integral Equations 52 4. 2 The Kernels and Their Singularities 55 4. 3 The Body and Its Volume 62 V NUMERICAL SOLUTION OF THE INTEGRAL EQUATIONS 69 5.1 P11/V Computation 69 5.2 P33/V0 Computation 73 5.3 M11/V0 Computation 75 5.4 Sampling Rate 78 VI CONCLUDING REMARKS 81 REFERENCES 86 APPENDIX: THE COMPUTER PROGRAM 87 iii

013630-9-T 1. INTRODUCTION When a plane electromagnetic wave is incident on a finite perfectly conducting body, or a plane acoustic wave incident on a finite acoustically soft or hard body, the scattered field in the far zone can be ext:rl.lpanded(. in a power series in the wave number k if k is sufficiently small. The determination of the first few terms in these series requires the solution of certain elementary potential problems. We here consider the potential problems associated with the first (Rayleigh) term in the electromagnetic expansion and the first two terms in each of the acoustic expansions, and show how in the case of a singly connected body of revolution all of these terms can be deduced from the solutions of just five potential problems. If the body is not singly connected, only the axial component of the induced electric dipole moment is affected, and for a body consisting of two separate parts, an expression for the modified component is obtained. Each potential satisfies a simple integral equation. Computer piroglrams are described for solving the equations by the moment method, and since most of the equations are of first order type, the computational procedures are rather similar to those of Mautz and Harrington (1970). The appropriate elements of the electric and magnetic polarisability tensors (Keller et al, 1972) are then computed, along with the electrostatic capacity and a quantity y related to this, and these are sufficient to specify the electromagnetic and acoustic.ac:ltctering for any direction of plane wave incidence and any direction of scattering. For relatively simple geometries, the entire computation takes about 3 seconds on an IBM 360 computer. In our presentation we first examine (Section 2) the problem of a plane electromagnetic wave of arbitrary polarisation and incidence direction, and

AFCRL-72-0162 013630-9-T The Numerical Solution of Low Frequency Scattering Problems By Thomas B. A. Senior and David J. Ahlgren\ The University of Michigan Radiation Laboratory 2455 Hayward Street Ann Arbor, Michigan 48105 February 1972 Scientific Report No. 12 Contract F19628-68-C-0071 Project No. 5635 Task No. 563502 Work Unit No. 56350201 Approved for public release; distribution unlimited. Contract Monitor: John K. Schindler Microwave Physics Laboratory Prepared For Air Force Cambridge Research Laboratories Air Force Systems Command Laurence G. Hanscom Field Bedford, Massachusetts 01730

013630-9-T 2. PERFECTLY CONDUCTING BODIES 2.1 FORMULATION Let B be a finite, closed, perfectly conducting body of revolution about the z axis of a rectangular Cartesian coordinate system (x,y, z). In terms of the cyli ii(d ic1l polar coordinates (p, 0, z ) where p + y, arctan -. x the surface will be described by the equation P - P(z) where p can be a rnultivnlued function of z as, for example, in the case of a disk or a re-entrant shape, but is never infinite and is zero outside some interval in z. Let r be the radius vector to an arbitrary point in the domain 2 exterior to B and let n be a unit vector normal to the surface drawn into V. A linearly polarised electromagnetic wave is incident with electric and magnetic vectors A ^ ik k. r 1 - ae A ^ ikk.r H = Ybe A A AA where k, a and b are mutually perpendicular unit vectors such that b = kAa; Y is the intrinsic admittance of the homogeneous isotropic medium (of permittivity c) exterior to B and a time factor e has been suppressed. 3

013630-9-T For k.smnllll but kr large, the resulting scattered field E, H can be written as ( Kleinman, 1965) ikr s e E,,~ - - - 47 r ikr s e - 47r r k2 1 ^(rA ) + I (rAm) k (rP) - r^(r m) and magnetic dipole moments respectively. (2) where p and m are the electric: As shown by Keller et al ( 1972), ( 11 (P33 - ll) ( = -Y b + (M -M )(b.z)z 11 M 3^1 (3) where Pll, P33, Mll and M33 are functions only of the geometry of the body. For a given body, Pll P33' Mll and M33 are constants whose values are as follows: (i) 11 B a x an (x - 1) dS ds (4) where 1 is such that V2 = 0 1 = X in 2 on B (5) 4

013630-9-T ~1 O(r-2) = O(r as r - o0 ft (:ii) P33 J B where 3 is such that I a z a (z~ 3) dS (6) V2 ~3 = 0 3 3 — z+7 in U on B; (7) -y is a constant chosen to make IJ ~ 3n dS - 0, B implying zero total induced charge on B, and ensuring that (8) -2 3 = O(r ) M11 ff B as r- o. (iii).^) n. x (x - 1) dS (9) where l is such that V 2 I= 0 /1 in V on B an an (10) ax an 5

013630-9-T -2 1. J ((r ) as r - oo M33 j= n (z -t) ds B (iv) (11) where 3 3 is such that v2 q - 0 ~3 az an 3n 3 = 0 (r-2 ~3 ~O(r ) in V on B (12) as r - oo. Although the values assumed by the potential function f3 on B are quite distinct from those of ir, a 1 and 3, nevertheless, as shown by Karp an an (1956) and Payne (1956), 1 M = P. 33 2 11 (13) This obviates the need for solving the potential problem (iv) if the only purpose for finding r3 is to calculate M33. There is one other electromagnetic quantity of interest and this is the electrostatic capacity C of the body in isolation. If the body is raised to the potential unity, the surface charge density is ado p e an s an (14) 6

013630-9-T where (0 is an exterior potential function satisfying the boundary condition =0 1 onB. (15) The electrostatic capacity is then equal to the total charge induced on the surface and is C = - ~ dS. (16) an B Note, however, that if all portions of the surface are not in electrical contact with one another, charge can no longer flow freely over the entire surface, and additional (mutual) capacities can be defined. In particular, such electrical separation has a profound effect on the calculation of P33 and the modifications that result when the surface is disjoint are discussed in Section 2. The five quantities listed in eqs. (4), ( 6 ), (9 ), (11 ), and ( 16) can be computed by solving five separate potential problems of a rather stanlda Xd nature, and the manner in which this is done is as follows. 2 Let V be some potential function satisfying V V = 0 outside and on B, and let V be the regular part of V. V is therefore an exterior potential and we can regard V - V s V1 (17) V- V =V as an incident potential. Green's theorem applied to the function V in the region K then yields V (_) 4 a(r') n' t — r dS 4 - R an B(18) 7

013630-9-T where R r - r'. If the boundary condition on the potential V is V(r) 0, r on B, (19) eq. (18) reduces to V(r) V (r) n V(r) dS. (20) B and since the integral exists for all r including points on B, we can allow r to lie on B and apply the boundary condition (19) to obtain _ 1 1r ar V (r) J R ' t V(r') dS', (21) B 3V which is an integral equation of the first kind for 3-. an If, on the other hand, the boundary condition on the potential V is - V(r) 0, r on B, (22) an - - eq. (18) reduces to v( r)-vi(r)+ 41 f () ( d' (23) B a /I \ and because of the non-integlrable sinrutllnl'ity of n-r (-) at r = r eq. (23) is valid as it stands only if r is not on B. To obtain an integral equation for V, we allow r to approach a point on B in the direction of the I 8

013630-9-T inward normal, in which case it can be shown that |r V(rt) a ()dS '= 27rV(r) + V(rt) (- I djS B B where the bar across the integral signs denotes the Cauchy principal value. Hence V(r) = 2V i(r)+ 4r 4V(r') -a(R) dS (24) 27 JJB for r on B, which is an integral equation of the second kind for V. 2.2 PROCEDURE FOR P 1 The solution of problem (i) and, hence, the computation of P is a straightforward application of the integral equation (21). If V x, Vs = then V = x- - V1 (say) (25) with V 1 0 on B. Since x = p cos 0 and the equation of the surface is independent of 0, the potential V1 must everywhere have the same 0 dependence as V i, implying Vl(r) = Vl(p, z) cos 0. (26) av1 This is true also of n, and we can therefore write 9

013630-9-T av (r) = - Vl(p, z ) Cos0 T (s) cos 0 where s' is arc length along a profile of the body. Moreover, wt 2 1/2 R x p-p') 2 + (z- z) + 2pp '(-cos ) } with (27) (28) =0 - 0' and since dS' p' d 'ds', the integral equation (21) now takes the 27r pcos 0 1 4 1 r 47r 0 s 1 2 r 0 f ~ 'orm s T1(') Cos - ptd0 ds 1 R p'K T (s )ds cos0 1 1 where the kernel is K1 K (p, z; p', z) = 1 K1(p, 7T ~ cos' dt COS R (29) 10

013630-9-T and the integration is along the profile of the body. The integral equation for T1 (s) is therefore p'K T (s')ds' = 2 p P 1 1 (30) which can be solved to determine Pll= /I B B T ( s ). In terms of this quantity x -n V(r)dS an ~~ 27T s /0 10 p2 s2 T(S)dds p cos 0 T1 (s)d0ds which reduces to Pll = 7 2 p T (s) ds 1 (31) 2.3 PROCEDURE FOR P33 AND C The solution of problem (ii), leading to the calculation of P33' involves two successive applications of the integral equation (21). In the first case we consider an incident potential i i V V2 =z (32) and seek the corresponding total potential V2 satisfying the boundary condition 2 11

013630-9-T i V = 0 on B. Since V is everywhere independent of 0, it follows that 2 2 av2 V2 and an are likewise 0 independent. We can therefore write aV an = T2 (st) and the integral equation (21) now becomes (33) 27r s 1 z = 4-r 7 f s 1 1 T (s ) p d ds' 'K T (st)ds' 0 2 where the kernel is K0(pz;p'z) K0 - Ko<Pz;ptzt) 0 | ^0 d R (34) The integral equation from which to determine T (s) is therefore 2 s J p' K T2(s)ds' = 27rz (35) 12

( 1 ':..t.l.1-9-T The second of the two basic problems is that in \\hicli the incident potential is i i V V3 = 1 3 (36) We again seek the total potential V satisfying the boundary condition V3 = 0 on B, and writing av an' n T3(s ) (37) the integral equation (21) takes the form s f pK T (st) ds = 27r 0 (38) from which T (s) can be found. 3 In problem (ii), however, Vi z+ y and V = -3 Thus, i i V = V2 + TV3 2 3 (39) implying V V 2 +,v3, (40) 13

013630-9-T and if we write a, = T(st) (41) then T(s) = T2( s ) + 'Y~T3(s) (42) where T ( s ) and T3 ( s ) are the solutions of the integral equations ( 35) and ( 38 ) respectively. The constant 7 is determined by the condition ( 8 ) for zero total induced charge on B, viz. * fJ a (vi-V) dS = 0. (43) B an But AJJn dS = JJn B B (z+ y) dS = 0 (44) as can be seen by application of the divergence theorem; moreover 2 T s VdS - B = 27r T(s) pd$ dS pT(s) ds * We are here assuming that the surface is not disjoint. 14

013630-9-T s I = 27 p T2(s)ds + 27r p T3(s)ds and hence, by virtue of eqs. (43) and (44), s J P T2(s)ds 0 Since T2(s) and T (s) can befoundfrom theintegral equations (35) and (38), the constant T given in eq. (45) now completes the specification of the surface field T(s), and in terms of T( s) P3= V (r) dS B s 27r zpT(s)ds I0 Hence s s P33 2 n zP T2(s)ds + 2 7 zpT3(s)ds. (46) A valuable by-product of the above analysis is the electrostatic capacity A valuable by-product of the above analysis is the electrostatic capacity 15

013630-9-T C defined in eq. (16). This fact is apparent on recalling that the determination of C requires us to find the exterior potential 0 satisfying the boundary condition ( 15) on B, and this can be accomplished using the integral equation (21) with V = 1 = 1 - o (s = - B) ( = 0 on B). so that The problem is therefore identical to the second of the two basic ones considered above, and indeed 0 3 implying ____ av3 aI0 OV3 = - = - T3(s) an an 3 where T3 (s) is the solution of the integral equation (38). Hence, when the body is at unit potential, the surface charge density as a function of arc length is P, = eT3(s) 3 (47) (48) (49) and the capacity C is s C = 2re 0 pT (s)ds (50) We observe that the denominator of the expression (45) for y is simply 16

013630-9-T C /(2 vT E), which ensures that ' can never be infinite. s Some simplification of the preceding results is possible. Since V2 and V s are both exterior potentials, Green's theorem can be applied to the region 3 2) exterior to B to yield the reciprocity relation V2 B avs dS 3 _dS = an ffv B av2 dS an (51) (Van Bladel, 1968). But av2 an A A T2(s) - n. z and S av3 3an T3(s) an and from the boundary conditions on V2 and V3 2i 3 V = -z 2 V = -1 3 on B. Substituting these into eq. ( 51), we have ff zT3(s)dS J (T2(s)-. dS B B which reduces to s 0 pT2(s)ds = zpT3 (s) ds 3 (52) 17

013630-9-T With the aid of this result, the constant 7 of eq. (45) can be expressed in terms of the surface field T ( s ) alone: 3 5 zpT3(s)d s -( - (53) s pT3(s) ds 0 i.e. s pY1s 2 zpT3(s)ds, (54) but whilst this reduces from four to three the number of separate surface field integrations involved in the calculation of P33' there is no way to avoidE entirely the determination of the surface field T ( s ). Indeed, the simplest expression for P33 is s P 2T J zp T2(s)ds - 2 (55) 33 E -0 2.4 PROCEDURE FOR M11 The solution of problem (iii) leads to the calculation of M and is a 11 straightforward application of the integral equation (24). If i x and V s -~ I 18

013630-9-T then V = x- t = V4 (say) (56) and aV4 /3 n 0 on B as a consequence of the boundary condition (10). Since x p cos 0 and the equation of the surface is independent of 0, V4 must have the same 0 dependence as V, namely, cos 0. In particular, on the surface V4(r) = V4 (s)os0 (57) which enables us to write eq. (24) as V (s) cos 0 = 2p cos 0 + 2 4 2 7T 2r s I f V(s' )cos0' — a( )p'dtdds on' 0 But n = s c ( C os ' + ysin ') - z sin a' ac"' tan-1 ap az (58) where (59) so that an_ (R) ~ (R) - - os a' (p -p cos ) - sin'z - z) (60) R3 19

013630-9-T with i = 0 - 0' as before. Hence cos '0 a )d = cos 0 cos27 1 d 0' = cos 0 cos 0 a an' 1 do d = 2 cos 00 p os' (z' - z) sina -p' cos a where -1 1 (p, z;p', z') 2 = 2 (p,z;p ' z') 2 2 7r 7T * cos ' d R3 3 d R ((31) (62) The integral equation from which to determine V( s ) is now V4(s)= 2p +F1 V4(s ) (Pcos' + 0 - z) sin a' - p cos a'] Q P p ds (63) and in terms of V (s): 4 M = ffn. VDdS 11 Y4(r) B 27r s = I cosaV4(s)cos 0pdods o0 J0 20

013630-9-T s i.e. M11 7r P V4(s) cos.ds (64) ^0 2.5 PROCEDURE FOR M33 Although it is not necessary to compute M33 directly because of the relation (13), the integral equation which the corresponding potential satisfies must be solved if the second term in the low frequency expansion for an acoustically hard body is to be evaluated. It is therefore appropriate to describe the determination of this potential function here. Once again we have a straightforward application of the integral equation (24). If V z and 3 then V z- 3 - V5 (say) (65) av_ with 5/an 0 on B. Since V must be independent of 0, eq. (24) -~- -- — p --- implies 27 S V5(S) 2z+ 2 V5(S ) anR P' ds But 2 27r 2 a- n t d o an' R d/ a do 0 00 21

013630-9-T 2 p cos, + [(Z z) sin at p coso Q where Q is as defined in eq. (61) and 7r ( z; ) - * (66) R R The integral equation from which to determine V5 (s) is therefore S s r V5(s)=2z+ f V5s ()pcosa + (z -z)sinat -ptcosa] Qi )P'dst (67) and we note in passing that M33 n. z V5(s)dS B 2 w s = - sin cV5 (s)p d ds f0 #0 i.e. M33 -2 PV (s)sinds. (68) 0 2.6 DISJOINT SURFACES So far it has been assumed that all portions of the surface are in electrical contact with one another, and if this requirement is not met, the analysis is no 22

013630-9-T longer valid. Thus, for example, an application of the above procedures to a body consisting of two separate spheres leads instead to the solution for the two spheres joined by an infinitesimal wire along the axis of symmetry, and though the presence of the wire (producing electrical contact) does not affect the values of M1 and P1l (and hence M33, by virtue of eq. 13 ), it does have a profound effect on P33 This is not unnatural since P is proportional to the longitudinal (z) component of the induced electric dipole moment. The breakdown in our formulation when B has several distinct parts stems from the imposition of the zero induced charge criterion (8). If charge cannot flow freely between the n parts B1, B2...., B, eq. ( 8) must be replaced by the n equations dS = 0, i1, 2,... n. (69) B i Since this obviously affects only the potential 1 3and leaves the procedure (and results) for Pll, M11 and M33 unchanged, our efforts will be directed at P33 alone with the objective of finding an approach which is applicable when B consists of just two electrically isolated portions B1 and B2. So that we may use to the fullest extent the work that we have already done, it is desirable to have this new appioachi as similar as possible to that appropriate when the two portions are electrically connected. By analogy with problem (ii) of Section 2. 1, the task is to find an exterior potent i:tl < 3 satisfying the equation V 2 3 - 0 in the domain V exterior to B, together with the boundary conditions z+ onB (70) 3 on B (71) 3 = 2 2 23

013(;3(-9-T where the constants 71 and y2are such that f 3n 3ff 3 dS = (72) B an 1 J J, dS. (73) 2 The quantity P3 is then given by eq. (6) as before. 33 Because the boundary conditions on B1 and B2 differ, it is no longer convenient to think in terms of incident and total potentials, with the difference representing the desired exterior potential. Let us therefore consider the basic potential problem in which ) is an exterior potential satisfying the boundary 3 condition on B1 (74) 0 on B2 By ipplic:i/iol of Green's theorem to the domain J, we have (r) 1 a (1) (rt )= S 3 ()J 47T 3 an' R R an' 4r a n R 47T R an 3 B 1 B 24

013630-9-T and the first integi:ll is identically zero since B1 is itself a closed surface. If, now, r is allowed to approach B, application of the boundary condition (74) gives 1 1 4 R B a, ()(r)d r' on B1 0, r on B2 (75) a 3(1) / an' which is an integral equation from which to determine. It can (1) /an c be simplified somewhat by observing that 3( and, hence, are are independent of the azimuthal coordinate 0. When the 0 integration is performed, eq. (75) reduces to T 3 ( W27r, It K0T3(1(s)ds ( 0 0 r on B1 r on B2 (76) c.f. eq. (38), where 3(1) - a (1)( T3 = an' 3 (77) and K is the kernel defined in eq. (34). It will be noted that the integration in (76) is over the entire profile of the body B B + B. Similarly, if p (2) is an exterior potential function satisfying the boundary condition 4(2)_ on B on B2 2 (78) 25

013630-9-T then s 0 r on B1 f PK0T3(2)( )ds 2r p'K T (s')ds 1 (79) 3 2r, r on B2 where T(2 (s) a (2 t) (80) Comparison of eqs. (76) and (79) with (38) shows that T( (s) + T()(s) T3(s) (81) 3 3 3 where T3 ( s) is that surface field quantity which is appropriate when B1 and B2 are electrically connected. If T (s) has already been computed, it is (1) (2) clearly necessary to compute only one of T( (s) and T (s). 3 3 Let us now return to the potential problem set forth in eqs. (70) through (73). As regal ds the boundary conditions (70) and (71), an exterior potential satisfying them is (1) V3 z - V2 13 +Y2 + 3(2 (82) where V2 is the total potential considered in Section 2. 3. Hence _ 3_ az (1) (2) an an 3 2 and since T2( s) is given as the solution of the integral equation (35), it only remains to specify the constants 71 and 72 From the zero charge condition (72) and using the fact that 26

()13(;:3 1-9-T B1 or B2 aZ dS 0 an we have 7 f pT3 ((s)ds + Y2 (1) (1) pT3()(s)ds - 3 f pT2(s) ds (1) (84) where the symbol (1) below the integral signs shows that the integrations are carried out over the profile of the portion B alone. Similarly, from eq. (73), (2) (2) pT3 (s)ds + 72 (2 (2) pT (s)ds = - 3 (2 (2) pT2(s)ds (85) where the integrations are over the profile of B2 alone, and if we now define C = 27rc 11 C2 = 2 r eqs. (84) and pT3(1) (s)ds (1) p T3(1)(s)ds (2) [ (85) take on the more C12 = 2 7re f p T3()(s)ds, (1) (86) C22 = 27re 22pt compact form P T3 2(s)ds (2) I 71Cll +' 2C12 = -2Cre I pT2(s)ds (1) (87 a) 27

013630-9-T 71 C21 + 2 C22 = -2 7 I pT2(s)ds. (87 b) (2) It will be observed that the quantities C11l, etc. all have the dimensions of capacity, and by virtue of eqs. ( 50) and (81 ), C11 + C12 + C21 + C22 = C (88) where C is the capacity when electrical contact is maintained. Rather than solve the eqs. (87) directly, it is more convenient to first eliminate the surface field quantity T2 (s) from the expressions. That this is possible can be shown by application of reciprocity to the exterior potential functions z - V and 3(2) From the pair z -V2 and 3(1) B B Hence A s r T J T2(s) dS = J zT3(1)(s)dS B B implying f s g pT2(s)ds zpT3 (s) ds. (89) (1) #0 Similarly, s pT2(s)ds = zpT3(2)(s)ds (90) (2) 0 and we note that by:dlit Lioin of the last two equations we recover eq. (52 ). 28

013630-9-T Finally, from the function pair 3() and () ff - B1 B1 (2) T3 (s)dS 3 I2 B2 implying (1) T3 (s)dS (1) pT3 (s)ds 3 (1) pT (2)(s)ds = 3 (2) (91) i.e. C12 = C21 (92) as expected. Using eqs. (89) and (90), T2(s) can be eliminated from the eqs. (87) (2)(s ) using eq. (81), we obtain and if we also eliminate T (s) using eq. (81), we obtain 3 ('1 -'2) c+1 - 2E (1 2 )C 21 +'2 27 Is p T3(s)ds = - 2 27 (1) 0 J pT3(s)ds - 27T (2) 0 (1) zpT (s)ds zT3)(s zpT3 (s)ds 3 s -2~re f 0 zpT3(s) ds 291

013630-9-T These can be solved to give s s s s (f zpT3 ds J p(s)ds- r zpT3)(s)ds fpT (s)ds} 1 - 2 ( zp T(s)ds J pT3sds - zp T( ds) s pT ds (1)0 where p J pT3(s(s)ds ds - )())ds. 6 T (3 d(s)ds - zpT 3 T3 (s)ds f (95) ((( 9 )(2) We can now proceed to the calculation of P. If we write this quantity * Pr 33T (S as P33 to distinguish it from the PQQ of eq. ( 55 ) for B1 and B2 in electrical contact, we have, from eqs. (8) and (83), 33 ~ T2(s)+ 7Y1 T3 ( )+ T 3 (2) S B - 2 / zpT TP 2(s)+ y2T3(S)+(7-.1 72)T3 (s)ds (96) 30

013630-9-T But p = 33 2 7r zp (T2(s) + Y T3(s)) ds ( see eq. 46 ) where T is given by eq. ( 53), and thus 33 33+ 2 2 2 ) s rO zpT3 (s) ds+ 27r (y1 -72) zpT ( )(s)ds 3 (97) Moreover, from eqs. (53), (94) and (95 ), after some manipulation, s [ pT3( )(s) ds 2 -7 -. (, -Y -l. (98) s r. pT3 (s)ds which enables us to write eq. (97) as P33 = P33 + 27r(y1 s 0 (z +7y) pT 3()(s) ds (99) The factor (7 Y1 - y2 ) is defined in eq. (93) and invoking yet again the expression (53) for T together with the identity (91), we have 1 2 - s if pT3(s) ds (1) p s pT (s)ds+ + zpT3 )(s) ds) 31

013630-9-T s 5 2(z+,pT )(s)ds 0 giving _( C P33 p33- E L|J ^)pT( sds} (100) where C and A are defined in eqs. (50) 'and (95) respectively. This is our final expression for P33. Compared to the situation when B and B2 are electrically connected, the only additional field quantity that must now be found is T (s), which is given as the solution of the integral 3 equation (76); and since C/c and A are both positive, electrical separation decreases the longitudinal component of the induced electric dipole moment. 32

013630-9-T 3. ACOUSTICALLY SOFT OR HARD BODIES 3.1 GENERAL PROCEDURE Let B now be a finite, closed acoustically soft or hard body of revolution about the z axis of a Cartesian coordinate system (x, y, z). It is of no concern whether B is disjoint or not. A plane acoustic wave is incident and its velocity potential is written as ui A ikk. r e - (101) A s where k is again a unit vector in the direction of propagation. If U is the scattered field that is produced, then U5 satisfies (V2+ k2) Us = 0 in V (102) (103) ( aur - s r ar )k U 0 as r-4oo, and the boundary condition Us on B (104) if B is soft, or au an aui an on B (105) if B is hard. Eqs. (104) and (105) are equivalent to * To avoid any possible confusion, we shall henceforth refer to U as a field. 33

013630-9-T U 0 on B (106) au a U 0 on B (107) an i s respectively, where U - U + U is the total field. A general expression for U(r) at an arbitrary point in 2\ is provided by the Helmholtz representation: 1, k Re ikR U(r) =U (r)+ i fU (r) a - - i R d B (108) where R = s where R s- |- r't as before. For sufficiently small k, U, U and, hence, U can be expanded as power series in i k of the form Ui(r) = (ik)m U (r), (109) m=0 and when these are inserted into eq. (108), the coefficients of like powers of i k on both sides of the equation can be set equal to give - - iU (r)=U(r - JJ7-m I I...-2(r)+ a m - m- 4' (m-t)' ' m 1) RnM S R.{=0 B - a' U(r') dS'? (110) for m 0, 1, 2,... By allowing r to lie on B, an integral equation is obtained from which U (r) can be found; and as is seen by substituting m — 34

013630-9-T the power series for U(r) into eq. (102), 2 2 V U U -= 0 0 1 (111) 2 V U = U m m-2 m > 2, showing that U0(r) and U (r) are potential functions, but U (r) is not unless U (r) - 0. In the far zone (r -- o ) the low frequency expansion of the scattered field deduced from eq. (108) is ikr U (r)-4 e 4 4r r OD m (m-f)' ( * r') m- L L 2 ^) -+1 J J ( 1.r ) m=0 fa0 B X n. r) U _(r) + U(r') dS (112) (Kleinman, 1965, with the correction of a sign error), provided U (r') is taken to be zero. Our objective is to calculate the first few terms in this series. 3.2 SOFT BODIES We now specialise the above results to the case of a soft rotationally symmetric body illuminated by the plane wave (101), and seek the first two terms in the low frequency expansion of the far zone scattered field. By 35

013630-9-T invoking the boundary condition U (r) = 0 m - on B, (113) m =0, 1, 2,...., we have u (r) ikr e (-1) 47 r U0(r'LdS' -ik B B a- U ( ant ' 0 B r') (114) a') dS'+ (k) U (rl) dST+ o(k2 anI1 -JP showing that only the potential functions U (r) and U (r) are required. From eqs. (101) and (107) it follows that i U0 (r) = I 0 i ~ u (r) 1 - A = k.r (115) and by inserting the boundary condition (113) into (110), the latter becomes m. ui (r),1 i 1 Sf U m-1 a.. U (r) = U (r) 4n ( - FR U(r')dS' m 4w ( —~* &nt B which, for r on B, reduces to i U (r) 5 m - m 1 r 4n 7 1 (m-I)! mn-i —1 a ffR m 1 U1(r')dS' B (116) 36

013630-9-T When m 0, eq. (116) gives B a U0(rL')dSt (117) This is identical to the integral equation satisfied by the potential V3 (r) of Section 2. 3, and hence U0(r) = V (r), U (r) - T3(s) an 0 - 3 (118) We note that U0 (r)dS s 2r J pT3(s) ds = 0 C E (1.19) B (see eq. 50), where C is the electrostatic capacity. Fromeq. (116) with m=1, A k. r ffJ Uo(r')dS f' =a () dS B B (120) and using eq. ( 119 ), the left hand side can be written as (k. x)pcos0 + (k. y)psin + (k. z)z 4 -Since the surface of the body is independent of 0, it follows that U 1( r) must have the form 37

013630-9-T (^AAA)(1) + (A (A A (2) C U(3) r U (r) k.x) cos0+(k.y)sin U) (r)+(k, z)U )(r) - 4U ( (121) \ where the individual U )(r), j 1, 1 2, 3, satisfy 1 P = 4r B1 R B a an' U )(r )dS, (122) 1 1 Z " 4 J J R B 1 = 4 J R a 1 (2) -a5; U (r) dS' an' 1 ( 123) a U(3)(r)dS an' U1 (rdS (124) B Comparison of eqs. (124) and (117) shows U(3) (r) - U(r) = a U3(r) = T3(r) -5 U3 - 9 (125) Similarly, U1( (r) cos J is identical to the potential V1 (r ) implying of Section 2. 2, U1) ) ( Vl(p, z) (126) so that a Uv(1 (r) Tl(s) an (127) 38

013630-9-T (2) and U1 (r) is identical to the potential V (r) of Section 2.3, so that a (2) a U (r) = T2(s) (128) an 1 2 It is now a trivial matter to evaluate the right hand side of eq. (114). The first integral is clearly C/e, and the second can be written as r.r ) T3(s )-.xco +(k.y )sin' T1(s') AA C +k. z) T2s ) 4- T3( dS' S. J )z + ] T3(st) - (k. z)T2(s )) ds But S s P'T2(s') ds' t= z''T3(s') ds' (52) 0 J0 and hence the second integral on the right hand side of (114) is s 2- - zr P T3(s) ds' e 47r e r k.j ) where y is as defined in eqs. (53) and ( 54). 39

013630-9-T The low frequency expansion of the far field is therefore r 4 - ikr U (r)- (- i - - y(r ]k) z + O(k ), (129) 4w r kL4E showing that a knowledge of C and y alone is sufficient to specify the first two terms. As demonstrated by Van Bladel ( 1968), a similar result obtains even for a body which is not rotationally synlnletric. 3.3 HARD BODIES The final case to be considered is that in which B is a hard rotationally symmetric body. The boundary condition on U ( r), m = 0, 1, 2,.... mis then n U (r) 0 on B (130) &n mand when this is inserted into eq. (112 ), the low frequency expansion of the far zone scattered field becomes ikr ) r (-1) 'ikJn. r)U (r)dS' +k2. r') U0(r') 47r r' - - - B B ft ^, 2 A,A -U(r ') (. r) dS - ik3 (r.r)2 U0(r ) B -(r.r')U (r + U2(r ) (. r) dS' + O(k4) (131) 40

,!...)-9-T As we shall see later, the first term O(k) is identically zero, and we therefore need U (r), U (r) and U (r) to compute two non-zero terms in the expansion. From eq. (110). and the boundary condition (130), an expression for U (r) at an arbitrary point r in V is m "- -" m Ui (r) U (r) + ( )! (1-m+\)R m - ~ - anti ' o-0 B (132) and in particular, when m = 0, U( r) i i, Ua( r ~ ): 0(r) = U0 (r) + 47 (r) an' d B (133) s i Clearly U0 (r) = U(r) - U0 (r) is an exterior potential function a s a i and - U ( r) = 0 on B since - U0 ( r) = 0. In addition, Dn 0 - an o s -1 U (r) vanishes more rapidly than r as r-) oo since there is no term O(k ) present in the expansion (131), and hence Us(r) =0 (134) implying U(r) U i) Uo(r) = U0(r) = i (135) From eq. (132) with m - 1, we have i 1 -U (r) = U (r) + (rf) ant dR) Bl ( dS' 1 an B (136) 41

013630-9-T which can be converted into an integral equation for U1(r) by allowing r to approach B. Because of the non-integrable singularity of the kernel for r on B, it is necessary to apply a limiting process, and if a bar across an integral sign is again used to denote the Cauchy principal value, we obtain U (r) 2 Ui (r) + 1 l(r) an' dS (137) B for r on B, where U (r) is given byeq. (115). In terms of the cylindrical polar coordinates (p, ~, z), U (r) (k. x) cos 0 + (k. y) sin p + (k.z) z and since the surface of the body is independent of 4, it follows that U1 (r) can be split up into three parts each of which has the i dependence of that part of U1 (r) giving rise to it. In particular, on the surface, U{(r) (k.)cosO + (I.~)sin) V(s) + (k.z )V (s) (138) where V ( s) and V5 ( s ) are the potentials introduced in Sections 2.4 and 2. 5 respectively and satisfying the integral equations ( 63 ) and ( 67 ). For the remaining function U (r) an expression at an arbitrary point r in 2 is given by eq. (132) with m = 2 and is U (r) - U (r) + 41 J U2(r) ar R dS (139) B 42

013630-9-T where (see eqs. 101 and 109) U (r) 1 ( r^)2 (140) 2 2 An integral equation for U (r) can be obtained by allowing r to approach B, but it proves unnecessary to determine U (r) explicitly if the only pur3 2 pose is to calculate the term O(k ) in eq. (131). S To see this we first note that since U (r) 0, the eqs. (111) imply S 0 0 V2 U2 S showing that U2 is a potential function. Moreover, from eq. (139), 2. U2 U - U2 is an exterior potential, being of double-layer type, and 2 2 2 since a i A AA a U2 (r) (k. r) (k.), (141) the boundary condition on US (r) is 2 - a A AA U (r) - (k.r)(k. n) (142) an 2 - 5 for r on B. U2 clearly depends on the direction of incidence as well as 2 that of the normal to the surface, and in principle nine separate but elementary potential problems must be solved to find U In terms of these potentials, 3 3 U2 (r) k.k.G..(ij) (143) 2 m - i j j il1 j=1 where, for convenience, we have put 43

013630-9-T x=x 1 y x2, z =X and the potential functions Gi j ( r) are such that a ( A a Gj(r) - -(n. i)(r.x) an i- - j for r on B. In like manner we can write 3 U -( r) =11 k. F.(r) 1 - If| 1 - i, jl, 2, 3 (144) (145) where the functions F. (r) are such that IL a A a F (r)= -n.x. an i- 1 (146) on B, and comparison of (145) with (115) and (138) shows that on the surface F 1(r) (S) -P - 4 F2(r) = (V (s)-p) F3(r) = V5(SL-z cos 0 sin 0 (147) Following Van Bladel (1968) we now apply reciprocity to the exterior potentials Ft (r) and Gij(r), I, i, j 1, 2, 3, in the region 2) to get B Gij ( r) B an F ) = JF( Fr') ann Gij(r ) ds B B 44

0 1 3 G 3 0-9 -T which reduces to IA (A(r9( A ij 1)'t. ) d~ = ~ I~r)( 1) rt j )d (148) B B when the boundary conditions ( 144 ) and ( 146 ) are employed. Hence B a i I 1kj fGij (rt)(n. XI ) d Sr B 14 44 A t A A ikj F, (rl)( i)(rl.x.) d S I 3 i B I A A T (r.1) (k.nl) (k. ) dS' A S B implying Jf U - )n. r)dS' B f U(t) )(At A )ds u 2-rn.r S B I B and A B )(A. A A A FI(r' r x,) (k. nl) (L rl) dSt +1 (1' (A 2 (k. rt)(. ) dS' B ( 149) 45

013630-9-T This integral is the only form in which U 2(r) enters the far field expansion 32 through terms O (k ), and since the F (r') are known by virtue of the eqs. (147), the integral can be computed without the explicit determination of U (r) itself. We are now in a position to evaluate the individual terms shown in eq. (131). Since U (r) = 1, we have 0 - ( '.r) Uo(r') dS' =0 (150) r0 -B 2 verifying that the leading term in the far field expansion is O (k ), and r. r') U0 r')('. r) = r. ( r.r ') dS B B r. V (r. r')d ' Vo V0 (151) where V0 is the volume of the body. Also, from eq. (138), ffU (r')(n'. r )dS' = r. cos cos 0' x + cos a' sin 0'y-sin 'z } B X \(k. x)cos0' + (k. y)sin O V4(s')+ (k. z)V5(s')p d 0'ds' r. {x x( k(. x)+y(k. p' (s )cos c' ds' r... / 46 -- ----— W I

013630-9-T s A A A. A -( r. z)(k. z) 2 I 0 p'V5(st) sin ae ds' where aCt is the angle defined in eq. ( 59). Hence, from eqs. ( 64) and (68), l U ) (n.)dS. rM - (k. z)(Ml-M33) 1 - 3 B ( 152) where M and M are the elements of the magnetic polarisability tensor discussed in Sections 2.4 and 2. 5 respectively. As we have previously noted, for a body of revolution M33 is related to P1l (see eq. 13). When the results of eqs. ( 150) through ( 152) are substituted into eq. (131), the low frequency expansion of the far zone scattered field is found to be US(r)-" e4 7r (k M( (k ^) (M - M33 V + O (k3) (153) where the actual term involving k is i/ ik3 ( r.r) U 0(r')- (. r') U1( r ' ) + U2(r'] (n. r) dS 2). B (154) Unfortunately, the evaluation of this is rather a messy task. Since U (r) = 1 (see eq. 135), (r. U2) U 0 ) =. B B n' (r. r') dS' a B B 47

013630-9-T Vo r. _ I V0 V ((r. r I)Tt dr' (155) To simplify the treatment of the next two integrals in (154), write F (r) FL(r) + r. x (156) so that (see eq. 145) U1 (r) I k F (r) (157) Using eq. /B B (149 ) we then have (-( r.r) U (r')+ U (r 1.__ (n'. r) dS' A k x ( r')(n. dS'm )^ (n k) - (k.)(r.r)( r)} dS'l W -L I. r A. A r.xf (k.r') 1 (A + (k. B r') (n. r) - (r. r) (n. k). rdS (k. r) dSt (158) But B 1 A )2t ^S 2(Jk r' ( n r~)dS' = ( k.r) k fr' d 7' V0 48

013630-9-T and r r.( n. r k. dS r +(k. r)k r 'd B V 0 as may be shown by analyses similar to that performed above. Hence J r') r(n'. r) - (rA ')('. (k. rt')dS' -. r J' dr' B V 0 (159) which cancels the contribution (155) of the first term in the integrand of (154). The complete integral (154) is therefore 3 A AA A ik3 (r) ( r.) ( k. r n.) - (k.x )(r. rt)(n '.r dS' (160) I B and to simplify this we now invoke the rotational symmetry of the body. From eqs. (147) and (156) we have Fl(r') = V4(s')cos 0, F2 (r') V (s')sin0', F (r) = V (s') 1- 4 2 - 4 3 5 (161) When these are inserted into (160) and the azimuthal integration performed, the contribution of the first term in the integrand is (r - (k.z)(r.z)(zcos psin (s)ds. -.z)(r.z) p (z cos a-p sinoa) V (s)ds 4 49

01:3 (63 0-9-T + r(r.z) 1 6(k.z J p V (s)cos ads p zV5 pzV5(s) sin ads. A A A^ 2 -27r (r. z)(k. z) s fo The contribution of the second term in the integrand of (160) differs only in having r and k interchanged, and when the two are subtracted, the final 3 expression for the term in k in the far field expression (153) is i3 AA ) A) P ik (k-).z r-(k. z p (z cos a p sina) V(s) ds - Jo4 s +(k. z)(r. z) 0 s -2 (k. z)(r. z) 2 p V5(s) cos ads 5 )z V5(s) sina ds 5 (162) Although this is only the 'second non-zero term in the low frequency expansion, it is much more complicated than the second term in the expansion for a soft body. The surface field quantities involved are the same as those associated with Mll and M33, but there is now no simple relationship analogous to (13) which enables us to dispense with V ( s ). If the direction of incidence or 5 50

013630-9-T A A A A observation is parallel to the axis of symmetry, i. e. k - + z or r = + z, the integral containing V4 (s) disappears, but there is no comnparable situation where the integrals containing V5 ( s) are absent except for the special case A A of forward scatter, r = k, when the entire expression (162) vanishes. 51

01:3;30-9-T 4. THE COMPUTATIONAL TASK When this study was first undertaken the main objective was to develop an effective program for computing the quantities Plp P33 and M 1 specifying the low frequency scattering behavior of perfectly conducting rotationally symmetric bodies. The realisation that the calculation of P33 produces as a by product the electrostatic capacity led us to add this to the list of quantities considered, but it was only later that the question of acoustic scattering came up. Since the first two terms in the low frequency expansion for a soft body are expressible in terms of C / c and zy, and y is implicit in the P33 computation, it was only natural to add this to our list, and for a hard body the first term involves no additional work. But the second term, ( 162 ), is another matter. In particular, it requires the explicit calculation of the surface field V ( s) that had hitherto been avoided by virtue of the relation ( 13 ), and even if this were done, the nature of the k term is almost such as to preclude any physical understanding of the data. For these reasons it was decided not to implement the computation of V5 ( s) and, hence, to ignore the second term ( 162 ) in the hard body expansion. The quantities which we are now left with are all ones which are needed for the electromagnetic problem. 4. 1 INTEGRAL EQUATIONS It is convenient to begin by listing the integral equations which have to be solved and the quantities to be computed from their solutions. Assuming that the profile p p ( z) of the finite, closed, rotationally symmetric body has been specified in some manner and its volume V0 computed as a preliminary step, then: 52

013 f; 30-9-T (i) solve I pt K1T1 (S) ds' - 2rp where the kernel K1 is defined in eq. (29); compute (163) p11 VO S 2 p T (s)ds 1 (164) (ii) solve s l p'K T2 (s')ds' = 27z (165) p'K T3 (s')ds' = 2r (166) where the kernel K is defined in eq. (34); retain the option to print out T (s); compute 3 C E "2.r p T3 (s) ds zpT (s)ds, 3 (167) E y= - 2 C (168) 53

013630-9-T 33 V0 ^ 27r r 2 v0 zpT2(s)ds - C V 0 0 0V (169) (iii) if and only if B consists of two separate closed parts solve s '21 P K0 T3(1) (st) ds' 0 B and B2 r on B1; (170) r on B2 compute S V0 0 2 (z+ y) pT3( )(s)ds 3 6 33 VO (171) - I p(1)(s)ds 2 p T (s)ds 27 3 c pT3 (s)ds (1) where the symbol (1) below the integral sign means that the integration is carried out over the profile of B alone (iv) solve V4(s) s ' 2 + -z)si n -Ppds f v4(st) cos aI 2 + (z.- z) sin aI -pt cosa] 1)Pt ds5 0 J T V4(S) - 2p (172) 54

013630-9-T where Q, Q2 and cr' are defined in eqs. (61), (62 ) and ( 59) respectively and the bar across the integral sign denotes the Cauchy principal value; compute M s pV4(s) cosads. (173) V V ' 0 0 J We therefore have four (five) integral equations to be solved, three (four) being of the first kind and one of the second, and five (six) derived quantities to be computed from their solutions: the numbers in parentheses refer to the unusual situationl where B is disjoint. Before attempting this task, there are certain features of the equations to be examined. 4.2 THE KERNELS AND THEIR SINGULARITIES The kernels K0 and K1 of the integral equations (163), (165), (166), (170) can be expressed in terms of complete elliptic integrals of the first and second kinds. From the definition of R given in eq. (28), we have R (+ 2 2+ (zz)2 )1/2 2 1/2 R (= p+p) + (z-z) ) (1 msin 2) / (174) where 4pp' -- (175) (p+p )2 + (z -z?)2 and ( - ). (176) 2

013630-9-T Hence -, 1 R 1 m 2 pp' )1/2 -1/2 (m sin ) (1 -m sin 0) (177) and when this is substituted into the definition ( 34) for K, we immediately obtain Ko m ppt )1/2 K(m) (178) where nr/2 K(m) = fO (2 - m in2 )1/2 (1-m sin 0) d (179) is the complete elliptic integral of the first kind (see, for example, Abramowitz and Stegun, 1964, p. 590). By a trivial manipulation, we also have cos J R implying cos ii -1/2 (mpp')1/2 R 2 m - 1 1 R 2pp' (180) {(i' - m (1-m sin2 ) 1/2 -(l-m sin 20)1 (181) and hence, from the definition (29) of K1, where E1 (mpp')1/2 (1 —) K(m)- E(m) 7T/2 E(m) = (1 -m sin2) 1/2 dO 56 (182) (183)

013630-9-T is the complete elliptic integral of the second kind (loc. cit. ). The above representations of K and K are exact. Since p, z, p', z' are all real with p, p' > 0, it can be verified that 0 < m < 1. Over this range E(m) is a finite slowly-varying function, having the values 7r/2 for m = 0 and unity for m = 1. A finite polynomial approximation sufficient for -8 computing E ( m) with an error of less than 2 x 10 is given in Section 17. 3. 36 of the above reference. Through the first three terms the precise expansion is (Jahnke and Emde, 1945): E(m) - 1- 4 m I + 0m + O(mI m r (184) 4 1 2 1 1 1 with m 1-m, (185) 2 2 i.e.(pp') +(z-z') 186) (P+p') +(z-z') and r= - In. (187) 2 m We observe that m = 0 if and only if p ' p, z '- z, that is, when the integration and observation points coincide. For an integration point in the immediate vicinity of the observation point, r2 )2 (188) where s is to a first order the arc length between the points. The elliptic integral K(m) also has the value Tr/2 for m = 0 but becomes logarithmically infinite as m -+ 1. A finite polynomial approximation sufficient 571 57

013630-9-T -8 to compute K(m) with an error of less than 2 x 10 is given in Section 17. 3. 34 of Abramowitz and Stegun (1964), and a precise expansion through the first three terms is (Jahnke and Emde, 1945 ): K(m) = + m + 2 ( 2189) 2 1m 4 O(ml ' m1 Because of the infinity of K(m) as m-+ 1 (ml- 0), K and K1 are also infinite in this limit, but their behavior in the vicinity of the singularity is easy to determine. Using (184) and (189 ) we have K 0T7 (F+m m r ) (190) and K 1 K1 72 + O(ml, 1, (191) showing that the singularity at p ' p, z ' z is an integrable one in each case. The contributions of the singular (or self ) cells to the integlrals in eqs. (163), (165), (166), and (170) are therefore finite and can be analytically approximated as follows. Consider for example the integral equation (165). If the self cell in the sampling procedure is centered on s = s (where p sp ) and is of arc n n length As, then s + As p'K0T2(s') ds' - 1 0T2( s ds' self s - s n 2 58

013630-9-T n T2(s n) 2 ( n) 1 s -n - As n 2 K ds 0 1 - As r8 n 2 21 AS /nI' s and hence self pt K T (s) ds' I T (s ) 22 n 16 Pn ( In + 1 As n As (192) It is desirable to retain the first correction, unity, to the logarithmic term to ensure the necessary accuracy when the sampling is relatively coarse and / or p is small. For the integral equations ( 166) and (170) the results differ (1) from the above only in having T3 ( s ) and T3 ( s ) respectively in place 3 n ( n of T ( s ); and for the integral equation (163): 2 n self p' K1 T1 (s')ds' ) - p T1 (s ) 1 n n + 1 A n +As s -- A s n 2 K dst 1 - As 2 / 8p T (s) In s \ ) -- As 59 2

013630-9-T giving I p' K1 T(S )ds ' T1 (s ) In n - '1 ) s 1 n As self (193) For the integral equation (172) the computation of the kernel is a more complicated task due partly to the presence of the functions 2 1 and 2. However, these also can be expressed in terms of complete elliptic integrals, and the resulting method of computation is much less time consuming than a direct numerical evaluation of the integral expressions for 1 and Q2 The definition of 1 is given in eq. (61), and using eqs. (177) and (181), the integraind can be written as I cos R3 3R (1-m 2 (1- n )(1-m sin 0) 2 -3/2 -21/2 -(1-msin 0) J 4 m -- from which we have 1 2 1 1 2m (1 -2 FT 2 -3/2 (1 - m sin 0) To d0-K(m). (194) To evaluate the remaining integral, differentiate the expression (179) for K(m) with respect to m to get 1 K'(m) - 2 7r/2 / sin2 0 2 3/2 0 (1 -m sin2 0) 60

1 2m 0 / 0 013630-9-T 2 1 1 -3/2 l(1-msin 0) (1-msin20) d 0 1 2m t/2 I (1-msin20)-3/2 dO 1 K K(m) 2m Hence 7T/2 ro -msin2 -3/2 (1-msin o) de = K(m) +2mK'(m) (195) and when this is substituted into eq. (194), the result is 1 i 1 pp 3/2 (1 Km K(m)1 1-"2' ) KT(m) - K(m) 24 (196) The procedure for n2 is similar. From eqs. (177) and (181), 2 2 cos = R3 1 2 2m m 2 -3/2 ppt i -(2-m) sin,- (2 - m)(1- msin2 ) implying -1/2 2 1/2 + ( 1 - m sin 0e) -3/2 (1-msin 0:) dO 1 2 2 m 3 /2 2 p t -2) r/2 I0 -(2-m)K(m)+ E(m), 61 / /

013630-9-T and when the expression (195) for the integral is substituted into this, we have / 3/2 2 {2m (l -) Kt(m)- 1 - 4 K(m)+E(m)). 2 m2 1PP' 2 (197) The finite polynomial approximations to E (m) and K (m) were mentioned earlier, and in particular, for the latter, 4 4 1 K(m)=(a +a +...+am + a m )+(b +blml +...+ m 4) In - 0 11 4 1 0 11 4 1 m 5 5 1 + O(m1, m1 In ) (198) 1 m1 where values for the coefficients a. and b., i 0,...., 4 are given in Section 1 17. 3. 34 of Abramowitz and Stegun (1964). Since d/dm r -d/dm1, it follows that 0 2 K (m) -+ (b al )+(b2 2a2 )ml +(b3 -3a )m +(b -4a) 2 3 1 -(b +2b m + 3b m +4b m ) In- + (m4 m n-) 1 2 1 1 4 1 m 1 1 n (1.99) which can be used to compute K' (m). We note the pole-like behavior of K t(m) when m 1 -m 0, and this is reflected in the non-integrable singularity of the kernel of eq. (172) at p'- p, z' - z. 4.3 THE BODY AND ITS VOLUME One of the many factors motivating the present study was the need to 62

013630-9-T compute the low frequency scattering behavior of missile-like targets. These are generally rotationally symmetric bodies ( or can be approximated as such to an accuracy which is adequate at low frequencies), and are often made up of several distinct parts, e. g. a cone mated to a cylinder which is terminated in a spherical cap. Although the complete profile of such a body is certainly not an analytic curve, each individual segment has a relatively simple equation whose form can be used to advantage in the nlUm!lelic:il process. It is therefore assumed that the profile is a finite piecewise smooth curve composed of straight line and circular arc segments. For definiteness, the number of segments is limited to 15 or less. At the end points of the profile where it intersects the z axis of rotation of the body, p = 0 (of course), and the nature of the program is such that segments which are perpendicular to the z axis can be handled, as can a 'disjoint' body having two separate parts provided each portion of the complete profile terminates on the axis. Every segment contributes to the total volume V0 which can be found by adding the individual contributions 6 V0. In certain cases, a volume contribution can be negative and subtract from the volume attributable to the other segments. Where this occurs, it must be noted as part of the input specification for the segment in question. In the following we list the input specifications of circular arc (Types 1 and 2) and linear (Type 3) segments, and give expressions for the corresponding volume conut i bu Ltions (assumed positive). The segments must be described sequentially starting at the intersection of the left hand segment with the axis, and the ordered sequence of segments defines the profile of the body. In some cases it may be desirable to regard a single linear or curved portion of the profile as two or more segments to permit a non-uniform 63

013630-9-T spacing of the sampling points over the whole. Type 1 Segment (Circular Arc, Concave Down) Specification: p1' 2 ( ) - p(zl) P2 3 P(Z2) (z(, p 2) ^\ \\ \ (ZO PO) 'v(z, 0(degrees), volume sense 0 < 0 < 1'.) If 0 is the angle subtended by the arc at the center of curvature, then the radius a is a -- --- p )+( 1 { 1 )2 + (z )2 1/2 2 sin - J 2 (200) Since we permit the specification of re-entrant circular arc segments we do not require z2 > z1. In order to obtain correct results for both standiard (z2 > z ) and re-entrant segments define the quantity 2 - 1i d - I2-Z l (201) 64

013630-9-T Then, the coordinates (z0, p ) of the center curvature are I I z0 ~2: P0 2 z z1 +z2 -d (P1 - p2) cot2 (202) P1 + P2 d (z - ) cot 0 2 The volume of rotation is given by z2 6V = v p (z)dz z1 and since the equation of the circular arc segment is -p0)2 2 (p-P +(z - Z ) 2 = a the incremental volume 6 V0 is 6V = 7 0 r2 2 (z- Z) + a 2 QO0 1 3 (U +u +u ) 2 1 2 1 + PO u2 (p2 P) - ul( P) + d a2 U2 = Z2 - 0 u1 = z1 - 0 (203) where 65

013630-9-T Type 2 Segment (Circular Arc, Concave Up) Specification: same as for Type 1 (Zo2, p) \ '~ (z2P 2) \2 (z1~ P) Eq. (200) gives the radius a of the type 2 segment, but the coordinates ( zO P ) of the center of curvature are now z 1 + z 2 + d(P lP2) cot 0 0 2 1 2 2 2 i +( +p + d(z.-z )cot20 -PO = 2 1Pi 2 2 - 1 2 (204) The incremental volu.me of the type 2 segment is 2 2 1 2 2 6V O = 7 ( z2 -Z1) P0 + aa (u2 u1 0 ( 2 12 + 212} (205) Note that only simple sign changes distinguish (202) from (204) and ( 203 ) from ( 205). Relationships that hold for both type 1 and type 2 segments may be derived by using a constant ~ defined as follows: 66

013630-9-T 1 type 1 segment 0 -. (206) 1 type 2 segment The center of curvature (z0, p ) and the incremental volume 6 V0 for circular arc segments of types 1 and 2 are then 1 r 0 o 2 + z2 + d (P - P2) cot P0 = 2 (P+2+d(z2-z1) cot20 6V0 = v (z2-) +a -( + (u +u22 )j ou2 (p O) - 2 2 02 -ul(pl-pO)+da20) where, as before, U2 = z2 -u1 z1- 0 (207) (208) Type 3 Segment (Linear) Specification: Z1 z; P1 P(Z1 ) (z2' P2) (z1' P1 P2 P( 2) 67

013630-9-T The equation of the segment is clearly P2 -PI and the volume contribution is (209) 6V = 2 7rT p dz i. e. V0 3 2 1 )(2 P2P1 +P2) (210) which is positive or negative according as z2 > z1, z2 < z1, respectively. 68

013630-9-T 5. NUMERICAL SOLUTION OF THE INTEGRAL EQUATIONS The numerical procedures involved in finding 33/V0, M11/V0 and where appropriate, 33/V0 are quite similar to those required for 11 /VO, and it is therefore sufficient to give full details only for P1 /V. 5.1 Pll/V COMPUTATION The primary task is the solution of the integral equation ( 163 ) for the function T ( s ) and this entails the determination of a sequence of values (i) T i 1, 2,...., N, approximating T (s) at the sampling points s = s. on the profile p " p( z). For this purpose the profile is divided into N cells Ci of arc length Asi and midpoints si corresponding to the coordinates (pi, z ) Within each cell we also define the points si and si+ where si- s - aO i (211) si+ si - Asi with the restriction 1 0 < <o. (i) By assuming that T (s) has the constant value T1 over the i th cell, 1 1 69

013630-9-T the integral on the left hand side of (163) can be evaluated as a linear com(i) bination of the T1 whose coefficients depend on the position (p, z) of the field point, leading to a linear system of N equations in N unknowns, viz. T1' ) p' Kds'+T( f pIKd8 +... +T1(N) p' Kds' 2rp C1 2 CN (212) ill, 2,...., N. Hence, the system to be solved is At = b 1 where t is a column vector with elements tli 1 A is a square matrix with elements (213) i 1, 2,...., No (214) f p'Klds' aij i, j "l 2,.... N, (215) Ci and b is a row vector with elements bj.= 2irp o b j = 1, 2, *.... No (216) Increasing the complexity of the quadrature technique used to evaluate the integrals f will generally improve the accuracy but will almost C. J 70

013630-9-T certainly increase the computational cost. What is therefore desired is the least expensive procedure capable of giving the required accuracy. The two simplest approaches are to integrate first and second order approximations to give (i j j): aj p K1(i, j) As (217) aij pW(i ) ( j Kl(iJI —)+ Pj+Kl(i, +) +Wo(ao) PjKi(ij)]Asj (218) respectively, where the subscripts j- and j+ correspond to the points s. and s of eqs. (211), and K (i, j) is the kernel defined in eq. (29) J J + 1 and evaluated at the points ( ) ( pp i,( z ). By requiring a <, we ensure that the sampling points s and sj+ do not coincide with the end points of the cell C., and thereby avoid any difficulty in the computations of 1 and 2 (see eqs. 196 and 197). When,U (219) eq. (218) reduces to the three-point Gaussian formula for which 4 5 w t w " (220) 9 ' 18 With this choice of w0 and wl, the advantages of eq. (218) vis-a-vis eq. (217) were now determined by computing Pll/V0 for a sphere using various values of N. Fig. 1 shows percent accuracy and C. P. U. time versus N for each integration scheme. It is apparent that for a given expenditure of C. P. U. time the Gaussian three-point technique is much more accurate than 71

5.0. 5.0 4.0 percent 30 error 2.0 1.0 0 4.0 3.0 C.P.U. time (sec.) 2.0 1.0 0 0 10 20 30 N 40 Fig. 1: Percent error and C. P.U. time of P1 /V0 calculation for a sphere: T denotes trapezoidal rule computation and G denotes three-point Gaussian. 72

013630-9-T the trapezoidal method, though the accuracies of both are severely degraded if N is too small (N < 5). Since the Gaussian scheme with N x 10 produces an accuracy of better than 99.8 percent for a sphere, there is no point in going to a more complicated procedure, and the computer program was therefore written using three-point Gaussian quadrature to determine the matrix elements a. In summary, the integral equation (163) is solved by conversion to the matrix system (213) in which ij [ Pj (i j-) +pj+K(iJ+) + pjK(i,) A i,j - l, 2,...., N; i,j (221) ai* [tn ( ) -1j^AS, i 1, 2,....,N. Having determined the sampled values T ) T s ) is 1 i 0 computed from eq. (164) by integration over each segment of the profile using a second order integration procedure (subroutine INTEG, described in the Appendix). 5.2 33 / V COMPUTATION The point sampling method of solution of the integral equations (165) and (166) requires us to find the sequences T ) T (s ) and T3) T ( s ) 2 2 1 3 3 i' i = 1, 2,...., N, from these equations. To determine the T, choose 2 CIO wO and w1 in accordance with (219) and (220) and thence solve the 73

013630-9-T matrix system A t = b where 0) t T 21i 2 1 1, 2,....,N (222) b " 2rz I i and ai, [4 (18 pK (i,) + P K iK(i+))+] 9 j As i, j- 1, 2,....,N; i j j, (223). ai [In F i + 1 As i 1, 2,...., N. (i) Asi / The Ti) are similarly determined by solving the matrix system At. b 3 3 where the elements ai are again given by (223), but T (i) 3i 3 i 1, 2,...., N. (224) b - 2. i The quantities C/, T and P33/V0! defined ineqs. (167), (168) and (169) respectively are computed using the same second order integration procedure employed in calculating 11/ V0 If the body profile consists of two discrete parts, it is also necessary to solve the integral equation (170). The corresponding matrix system is 74

013630-9-T almost identical to that in (224), and from the sampled values T (s i ) ( i andT (s ), 33/V (see eq. 171) is computedand, hence, P33/V0. 5.3 11/V0 COMPUTATION The basic approach is similar to the above in spite of the more complicated integral equation (172) that must now be solved. The matrix equation for the sampled values V4(s ) * V( is Avv b where 4 1 4 4 V4i 4 i " 1, 2,....,N (225) b " 2rp and a - j + P )j + i f(iJ) Asj aij 18 P f(i,-) + pj+f(i9+)+ f(i) i, ^ 1, 2,....,N; i ^ j, (226) 1 - p'f(s)ds' 1, 2,....N in which f(ij) [icosaj 2(i.J) + zj-zi)sin j-pjcos3 nl (i.j) iJ " 1, 2,.....N;. (227) I#j ls2 o i j. 75

013630-9-T We observe that the computation of each diagonal element of A requires the numerical evaluation of a Cauchy principal value (denoted by the bar across the integral sign in the above expression for ai ). As an approximation to this principal value, we remove from the cell C a slice defined by the interval 1 1 (s. -- As s + 3 As) where 3, 0 <,3 < 1, is the fractional 1 2 i1 2 ~ i exclusion; /3 1 implies no exclusions, i.e. that the principal value is not taken. We now have i 2 i aii tx - J i 2 i 1 p' f(s')ds' - fp'f(s') ds' 8 + As i 2 i (228) and these integrals are Defig si2 il also computed using three-point Gaussian quadrature. - st - i (1+ 3) As " si2 - (0 1 3) As8 12 2 0 1 a = 8 + 1 (1,8) A a i5 si + 4 (1+/) Ai 814 i5 - ao(1- /) Asi i6 si5 +2 a0(1- 3) As^ (229) 76

013630-9-T we obtain 1 a r s (1 - ) 2 ii2 8 (Pilf(,i)+Pi3f(ii3)+Pi4f(i 'i4) ", < + P f (it 1I ) +9 i2 f '( 2, 5 + I ' i A (2 - +p6f(ii) +4 (p 2f(i i 2)+p15f(i i ))] As1 (230) Equations (225) through (230) completely describe a system of N linear equations in N unknowns V4i, i x 1, 2,...., N. Their solution and sub sequent integration of the V4i according to eq. ( 173) yield Mll/ Experiments were performed to find an appropriate value for the fractional exclusion 3. As an example, for a sphereM 11/V0 1.5) with N a 20, the data in Table 1 were computed. If we exclude the fortuitous (?) Table 1 11/Vo 0 1.0 0.1 0.01 0.001 1.480 1.516 1.503 1.501 percent error -1.33 1.07 0.18 0.08 error zero occurring for, somewhere in the range 0.1 < 1 < 1, these data indicate that the choice 3B 0. 001 is sufficient to keep the error less than 0.1 percent. 77

013(;30-9-T 5.4 SAMPLING RATE Increasing the number N of points at which the surface is sam mipled will generally increase the accuracy of computation, but since the number 2 of matrix elements increases as N and the cost of a linear system 3 solution increases roughly as N, this improvement is obtained at the expense of an increase in computation cost. Unfortunately, there is no rule for specifying the minimum value of N sufficient for a given accuracy, and the information which follows is based only on our experience in using the program. The results in Fig. 1 and Table 1 show that for a sphere 11/V and 11 /V are accurately determined with N as small as 20, and this is also 0 true of 33 /V. On the other hand, if the body has a discontinuity in 0 '. dp/dz lying off the axis, it appears necessary to increase N to 50 or p more to nlintrlail] the same accuracy (error < 0. 5 percent) in the 11/V0 and 33 /V0 computations. This is illustrated by the results in Table 2 for Table 2 N N1 N2 11/Vo N1 2 11 7 4 2.752 20 15 5 2.801 40 30 10 2.872 70 50 20 2.888 Note: N is the number of sampling points on the generator of the cone (linear segment) and N2 is the number on the (half) base (circular arc segment). 78

013630-9-T a rounded cone with half angle 15. The small but not negligible (0. 58 percent) change in 11 /V0 as N increases from 40 to 70 suggests that such large values of N may be essential for bodies such as this for which T1 ( s ), T ( s ) and T3 ( s ) have infinities at one or more points on the profile. For the same rounded cone, the results for 11 /V0 are given in Table 3. Since an increase in N from 17 to 35 produces only an insignifiTable 3 N N1 N2 11 /V0 17 10 7 1. 680 35 25 10 1.678 cant change in M11/V0, the choice N 20 is now adequate. Observe M that the surface field V4(s) associated with 11 /V does not become infinite at a discontinuity in dp / d z, and this is undoubtedly the reason why in many cases a small value of N now produces the same accuracy as does a much larger value in the P11/V and 33 / V computations. No attempt has been made to exploit this finding in the general program. When tleating bodies composed of several segments, a strategy which has proved successful is to divide all segments into cells of approximately equal length. This serves to fix the allocation of any given number N of sampling points among the various segments. Tests so far performed have not conclusively shown the advantages of dividing a single segment into two or more smaller segments so as to effect a non-uniform sampling. It is, however, believed that such a sub-division may, for a given N, improve 79

013630-9-T accuracy in the 11 /V0 and 33/ V computations for bodies like the rounded cone having infinities in the surface field quantities. 80

013630-9-T 6. CONCLUDING REMARKS We have here considered the low frequency scattering of electromagnetic and acoustic waves by axially symmetric bodies. By concentrating on certain qjLulalltities. such as the normalised components of the induced electric and magnetic dipole moments, we have shown how it is possible to arrive at rather elegant expressions for the far zone scattered field in terms of quantities which are functions only of the geometry of the body. Each such quantity is expressible as a weighted integral of an elementary potential function which can be found by solving an integral equation. A computer program has been written to solve these equations by the moment method and to calculate the dipole moments, the electrostatic capacity, and a tl'u thle parameter y related to the capacity. Any body can be treated whose profile is made up of straight line and circular arc segments and it is even possible to have two distinct bodies with or without an electrical connection between them. Although no serious attempt has been made to optimise the program, only a few seconds are required to compute all of the above quantities to an accuracy of better than one half percent. We have already used the program to compute the scattering from a variety of shapes, and it may be helpful to list some of the results obtained so far. Data for a rounded cone consisting of the intersection of a cone of half angle 0 with a sphere centered on the apex are given in Table 4. L/w is the length-to-width ratio of the body. For 0 < 90 0, the values of P11/V and P33/VO are quite similar to those previously computed by Senior (1971) using a mode matching method, but since M11/V showed significant discrepancies, this quantity was determined for a variety of 0. Detailed checkilng has confirmed that the present data are accurate 81

013630-9-T 0 (deg. ) l/w 3 7.5 15 30 45 60 80 90 93 99.2 108 120 140 151.71 180 9. 554 3.837 1.932 1. 000 0.7071 0. 5774 0.532:1 0.5000 0.5262 0.5799 0.6545 0.7500 0. 8830 0. 9402 1. 0 Table V0 0.002870 0.01792 0.07137 0.2806 0.6134 1. 047 1.731 2.094 2.204o 2.429 r 2.742 3.142 3.699 3.938 4.189 4: Rounded Cones Pll/V P33/V 2.865 3.664 8.147 3.494 4.520 1.931 M 11/V 1. 884 1.813 1.678 1.484 1.366 1.312 1. 334 1. 373 1.386 1. 416 1.458 1.507 1. 547 1.540 1.5 5. 406 6. 386 7.428 7.303 7.123 6. 889 6. 586 6.441 6.283 4.428 4.368 4.261 4.071 3.789 3. 370 3. 187 3.0 2. 184 2.242 2.372 2.553 2. 769 3. 006 3. 042 3.0 82

013630-9-T to three significant figures. When 0 = 900 the body is a hemisphere for which precise values of 11 /V and 33 / V are available: 11/V0 4.430...., 33/V0 2.189.... (Schiffer and Szego, 1949, p. 152). The corresponding values in Table 4 are within 0. 2 percent of these. For 0 > 900 the cone is a re-entrant one, i.e. a sphere with a conical region removed, and when 0 1800 the body is a sphere for which exact data are also known. Results for ogives and symmetrical lenses whose arcs subtend an angle 0 at their centers of curvature are shown in Table 5. The transitional shape is a sphere for which 0 x 180. To illustrate the computations when two bodies are present, Table 6 gives data for two identical spheres separated by a distance e d where d is the sphere diameter. When the two spheres are touching (c = 0) itisknownthat 33/Pll P 8/3 (SchifferandSzego, 1949, p. 154); the ratio deduced from Table 6 is 2. 678, which is within 0.4 percent of the exact value. As e increases, ll / V, M11 / V and P33/V 0 0 rapidly approach the values appropriate to a single sphere in isolation. P33 /V0, on the other hand, is proportional to the axial component of the induced electric dipole moment for two spheres which are electrically connected by an infinitesimal wire, and with increasing c this increases indefinitely, as expected (Kleinman and Senior, 1972). The same is true of C/( Vx The parameter y has also been included in Table 6, and 0 since its exact value can be shown to be -(1 + e/2 ), the accuracy of computation can be judged. 83

013630-9-T Table 5: Ogives and Lenses Shape 0(deg. ) l/w 22.8 10.02 36 6.314 ogive 56 4.011 88 2.475 132 1.540 150 1.303 sphere 180 1 64.4 0.6297 43.6 0.4000 lens 28 0.2493 17.2 0.1512 11.4 0.09981 v Pll /V 0 0 0.004146 2.089 0.01058 2.100 0.02650 2.189 0.07148 2.363 0.1966 2.647 0.2847 2.775 0.5236 3 1.212 3.779 2.586 5.182 6.448 7.649 17.30 11.88 39.55 17.49 P33/V0 49.88 24.15 12.57 6.778 4.136 3.595 3 2.161 1.674 1.390 1.225 1. 144 1 /Vo C/(e w) 1.943 6.128 1.910 5.696 1. 843 5.451 1.739 5.424 1.611 5.696 1. 564 5.880 1.5 6.283 1. 368 6. 570 1.252 7.375 8.758 1. 098 10.759 1. 061 12.982 84

0136;30-9-T Table 6: Two Spheres P ~ 11/V O _33/V D 0 0 0.005 0.01 0.02 0.05 0. 1 0.5 1.0 5.0 2. 702 2.706 2. 709 2. 715 2.732 2.759 2. 891 2.950 2.994 7.237 7.280 7.319 7. 402 7. 655 8. 026 12.02 18.19 120.6 7.237 5. 021 4.800 4.557 4.210 3.922 3.299 3. 142 3. 045 3.046 11/V( 1. 605 1.586 1. 607 1. 592 1.579 1. 528 1.511 1.500 1.501 c E~T-)M y 6. 153 6. 151 6. 150 6. 148 6. 141 6.130 6. 015 5.822 4.384 3.470 -1. 0004 -1. 0029 -1. 0054 -1. 0100 -1. 0250 -1. 0503 -1. 2502 -1. 5001 -3. 5001 -6. 0000 10.0 3.002 383.1 85

013630-9-T REFERENCES Abramowitz, M. and I.A. Stegun (Eds. ) (1964), Handbook of Mathematical Functions, NBS Appl. Math. Series No. 55. Forsythe, G. and C.B. Moler (1967), Computer Solution of Linear Algebraic Systems, Prentice Hall, Englewood Cliffs. Jahnke, E. and F. Emde (1945), Tables of Functions, Dover, New York. Karp, S. N. (1956), Relation of Electric and Magnetic Dipole Moments (abstract only), Quart. Status Report No. 1, Contract AF 19(604)1717, Institute of Mathematical Sciences, New York University. Keller, J. B., R. E. Kleinman and T. B. A. Senior (1972), Dipole Moments in Rayleigh Scattering, J. Inst. Math. and Applics. (to be published). Kleinman, R.E. (1965), The Rayleigh Region, Proc. IEEE 53, No. 8, 848-856. Kleinman, R.E. and T. B.A. Senior (1972), Rayleigh Scattering Cross Sections (to be published). Mautz, J. R. and R. F. Harrington (1970), Computation of Rotationally Symmetric Laplacian Potentials, Proc. IEE (London) 117, 850- 852. Payne, L. E. (1956), New Isoperimetric Inequalities for Eigenvalues and Other Physical Quantities, Comm. Pure Appl. Math. 9, 531-542. Schiffer, M. and G. Szego (1949), Virtual Mass and Polarization, Trans. Amer. Math. Soc. 67, 130-205. Senior, T. B. A. (1971), Low Frequency Scattering by a Finite Cone, Appl. Sci. Res. 23, 459-474. Van Bladel, J. (1968), On Low Frequency Scattering by Hard and Soft Bodies, J. Acoust. Soo. Amer. 44, 1069-1073. 86

013630-9-T APPENDIX: THE COMPUTER PROGRAM The program computes 1/V0, C/P, y, 33/V 11 /V and, where appropriate, 33/V, and consists of a main program and six subroutines. A. 1 DATA SET A data set is made up of one control card and a number of segment specification cards, one for each segment (or sub-segment) of the profile. The segment specifications conform to the convention stated in Section 4.3. Control Card Columns Description 1 3-4 The number (1 or 2) of bodies. Two digit integer (right justified): the number of segments on the first body (the body to the left). When there is only one body, use these columns. Same as columns 3-4, but for body to the right. A printing key: 1: print T from 33/VO com3 0 putation. 0 or blank: do not print T3. 6-7 9 87

013630-9-T 11 A computation key (0, blank or 1 ) 1: suppresses computation of P11/V C/, P33/V0 13 A computation key (0, blank or 1) 1: suppresses computation of 11/V0. 21-30 A real number: the fractional exclusion J3. If these columns are blank,, defaults to 0.001. Segment Specification Card Columns Description 1 - 2 Two digit integer (right justified): the number of sampling points or cells on the segment. 4 Segment type key: 1: circular arc, concave down 2: circular arc, concave up 3: linear. 6 Volume sense: + or blank: additive volume -: subtractive volume. 11-20, 21-30 Two real numbers: respectively, the end coordinates z and z2 of the segment. 1 2 88

013630-9-T 31-40, 41-50 Two real numbers: respectively, the end coordinates P1 and p2 of the segment. 51-60 A real number: for circular arcs, the included angle in degrees. There are the following restrictions: (i) the total number of segments must not exceed 15, and (ii) the total number of cells over all segments must not exceed 80. The profile is specified in the direction of increasing profile-length, beginning at its left-hand intersection with the z-axis and ending at its right-hand intersection with the z-axis. Re-entrant segments are permitted, allowing z1> z2. 89

013630-9-T A. 2 MAIN PROGRAM The main program reads and prints data and supervises all computations. A rough flow chart showing the interaction of the subroutines is given below DATA For each segment 1) read specification card Solve linear systems 90

013630-9-.T RFAL APlI (80,,80),AP33(8n,80),AMI 1(80,80),X(80),8(80), 6 Z7F P (2?),HOFP (2),ST3 (A),T3 (80),MllI TI (80) INr\TF(,Fk NUMPTVS (15 ), PL I V/'I+'/,tRL/' I/,INOX(2) COniMOVN RHo(80,t9 ),(8099),ARC,(80)t,C(8O9),S(80,9)/SUL/IPS(80) DA TA MI NiWIWOPI vP I,,woWI/' -' t,A283I853,3141593,p.44444449.27:7777771 37 kFAD)(5,34,FNI)=999 )N800DNSLN~S2, IPkljNTKFYP11,PKFYM11,FR 3~4 FnPMAT(11IX,2(12,lX),3(11,1X),6XFl2.7) WR I TFE( A,t4) NFROONS1I 4 FORMA T('I BEG I NN I NG OF I)ATA S ET:'/'0 ',5X t 1ODI E S I, X,''2/ I F (NROO.LE. 0.0k. NB(1D.GT.s 2) GO TO 990 NSEGS=NS1+NS2 IF(NBOn.EQ. 2) WRTTE(6,I)NS? I FORMA T(' I,5X IBODY #2?',95X 1',2) IF( NBOf). EQ. 1) GO) TO 10 IF(NSI.LF. 0.0k. NS?.LF. 0) GO TO 990 10 TF(NSFGS.LE. 0 nOR. NSEGS oGT. 15) GO0 TO 990 WRITE(69352)IPRINTKEYPII,vKEYMI1 352 FORMAT(' I,5X,'PRINT KEY',3Xv'='I?1/' 'v5X,1COMP KEY P''S=',V'2/ 6' ',S1X,'COMP KEY M''S=',12) IF( ER.LE. 0..OR. FR GT. 1.) FR=.001 IE(KEYM11.EO. 0) WRITE(6,5)FR 5 FORMAT(' ',w5Xr'EXCLIJSI0N = IvF7v4) IF(KEYP11.NE. 0.AND. KEYM11.NE* 0) GO TO 990 M=0 NC 1 =0 v0=0.0 00 11 I=1,NSEGS RFAD(5vl?) NIJMPTS(I),ITYP,TSIG7NZEPRHOEPTHETA 12 FOnRmAT( 12,LX,Ii,1XAl,,4X',5FlO.7) TF(NIJMPTS(I).LE, 0.OR* ITYP.LEv 0.OR. ITYP.Gi., 3) GO TO 990 IF(ISIGN.E0, BL) ISjIGN=PL(JS WRITE (6,1l3)1,NUJMPTS(1I),1TYP, ISIG-Nr7-EP,tRHOEP 13 FORMA(I(OSEGMENT #',I2,':'/ I 15Xv'CELLS',7X,'=',12/' ',5X, A 'TYPE KEY',4Xj'=',I2/' ',5X,'VOLIJME SENSE= ',A4/' ',t5Xv 8 'z-cnFlRDINATE END POINTS =(' #F12.-7v',',Fl2.7,')*/' 1,5X, 5 ORHO-COOROINATE END P0INTS=(',F12.7v1','F12.7,')') IF(ITYP.NE. 3) WRITE(6o,14) THETA 14 FORMAT(' ',5X,'THETA (DEG) =',FLO.5) THEFTA=P I`THETA/l180. M=M+NUMPTS( I) IF(I -LE. NSIL) NCl=NCl+NLJMPTS(I) IF(M.GT. 80) GO TO 990 CALL DA4TA( ITYPNM,ZEPRHOEPTHETAFRPVOLINC) IF( ITYP.NE. 3 *AND. ISIGN.EQ. MIN) VOLINC=-VOLINC 11 VO=VO+VOLINC WRI TE(6A,52)VO 52 FQRMAT('OCOMPUTED RESULTS:'/* ',5X,'VOLUME',p6X,'=',FlO.5) 91

013 63 0-9-T Dnn2N~, I NDX( (1)=Nl A N =AR C (NI) TN=RH0( N,8) F)fl 3 L=NvM IF(L.FO. N) GO TO 8? AL= ARC (L) TL=RHn(, 8) IF(KFYM1l-l)llOl09q109 1.10 AMlil(NL )=0.0 AMi 1 (L,PN) =0.0 1.09 IF(KFYPII-l1)111,1129112 I A PII( N,9L )= 0.O APLL1(L,N)=0.0 AP33( NvL)=Oo0 AP33( LgN)=0.0 112 INDX(2)=L 92

013630-9-T Dfl 103 J=1,3,2 P6=,l + 6 fo 104 LL=l,2 I =3-LL I1=INnX( LL) I?= IN X ( I) TI2=RHn (I2,lP6) CALL SETUtP(KEYPll,KEYMll 1,1,2,JP6,AP11l1,AP133,AMI 11,1) IF (KEYP11-1)105,106,106 105 APll( I1, I2)=APll( I1,I2)+API11ll*TI 2 AP33( I 1,12)=AP33(I1,1 2)+AP133*TI2 106 IF(KEYM1l-1) 107,104,104 107 A 11( I 1, 2)=AMll( I1,I 2)-AMI11*T12 104 CnNTINUE 103 CIN TINUF CALL SE TIJP(KEYP11,KEYM11,N,L,8,AP11,AP133,AMIll1 ) IF( KEYP11-1) 108,209,209 108 H=WO*API 11 AP1 (N,L)=AL* {WI*AP1 (N,L )+i)*TL ) AP11 (L,N)=AN*(Wl*AP11(L,N)+U*TN) U=W0O*API 33 AP33(NL)=AL*(W1-AP33(N,L )+U*TL ) AP33(LN)=AN*(Wl*AP33(L,N)+UJ*TN ) 209 IF(KEYM1l-1)210,3,3 210 AMll(N,L)=AL*(Wl*AM11 (N,L)-AMI111*WO*TL) CALL SETU)P( lOLN,8tAPIllAP133AMI11,O0) AM1l (L,N)=AN*(W1*AM11 (LN)-AMI1l*WO*TN) GO TO 3 82 IF ( KEYP11-1)83,84,84 83 tJ=ALOG( 16.TN/AN) APl1 (N,N)=(IU-1.0)*AN AP33(N,N )= (IJ+1.0) *AN 84 IF(KEYM11-1)85,3,3 85 IF(FR.EQ. 1.0) GO TO 3 Dn s86 I=1,6 86 CALL SETUtP( 1,O,N,N,tI t),Ut)ST3 ( I ),1) lU=.5*( 1.O-FR) *AN AM1 (N,N)=PI-lIJ*(WO*(RHO(Nt2)*ST3(2)+RHO(N,5)*ST3(5) ) 3 +W1*(RHO(N, 1 )*ST3(1 ) +RHO(N,3 )*ST3(3)+RHO(N,4)*ST3(4)+ 4 RHO(N,6)*ST3(6))) 3 CONTINUE 2 CON TI NUE 93

01363 0-9-T MI1 20 1=1,M 20 F ( I )TW IPIP-,RH(1( 1,v8) IF(KFYP11-1)21,24,24 21 CALL DFcnmp (APllM) CALL snLVE(AP11,XF~,M) nn 22 I=lM 22 X( I)=RH0( I,8)**2*X( I) CALL JNTFG(XvNSEGS9NUJMPTStPll) P1ll=PlI *PI /VO 24 IF(KEiYMl1-1 )25,28,28 25 CALL DFcnmp(AM11I,M)_ C ALL Ls L-vF.(AMli,XI~,M) MI 26 T=1,M 26 X( I )=RHn( i,8)*c( i,)*x( i) CALL INTEG(XNSEGSPNIJMPTSMll) 28 I F (KEYPi11-1)32,45,45 32 nfl 29 I=lM 29 B(I)=TWOPI*Z(II,8) CALL DECOMP (AP33 9M) CALL S0LVE(AP33,XBM) nfl 134 I=1,M X( I)=Z(1,v8 )"RHO( 1,8 )X( I) 134 B(I)=Twnpi CALL INTEG( XNSEGSN(JMPTSP33) CALL S0LVE( AP33,T3tBM),DOl 35 1=1,vM X( I)=RHO( I 8)*'T3( I) 35 B( I)=Z( 1,8)'X( I) CALL INTEG(XNSEGSNUMPTStCAP) CAP= TWOIP I*CAP CALL INTEG(FB,NSEGSINUMPTSGAM) GAM=-TWOP I*GAM/CAP P33=( TW0PI*P33-CAP*GAM*GAM)/V0 IF(NBnD - 1) 39939,54 54 DO 36 I=1,PM IF( I-NCLI)305,305,306 305 Pi(1)=TWOPI GO TO 36 306 B( I)=000 36 CONTINUE_ _ _ _ _ _ 94

013630-9-T C ALL S~OLVE( AP33,TI,B,9m) Dn 307 1=1 M -307 X( I)=RHl( I,8)"T1 (I) CALL lNTFG( XNSFG7SNUMPTSTL) CIALL TNTE(( X rNS1,N(iMPTSvTN) DO 308 T=lm 308 X (I)=Z (1,8)*X (I) CAL L I NTEG ( X,9NS GS, NtJMPTS,1J) nEL TAP=-( Twnlpi/VO)*(J+CGAM*TL )**2/(TN-TWOPI *TL*TL/CAP) U=P33+DELTAP 39 WRITE(~,v40)CAP,GCAMP1llP33 40 FORMAT(' 11,SX,'CAPACITANCE =1,F1O.5/' ',5X,'GAMMA',7X,'=~9F10.5/ 8 ' ',5X,'P11/V',7X,'='vFl~o5/' ',5Xv'P33/V',7X,'=',FIo.5) I F (NBOD EO. 2)WRITE(69309 )DELTAP9tJ 309 F~fnRMAT(t ',5X,'DELT P33/V',2X,'=lFIO.5/' ',5X, IDISJNT P33/V=I, 6 F10.5 ) 45 I F( KFYfA1L-1 )4?,337,337 42 WRITE(6,43)Mll 43 FORMAT($ g,5X, 'Ml/V3,7X, 6=IFl~O5) 337 IF(IPRINT.EQ. 1 *AND, KEYPll.EO. 0) WRITE(6,44) 6 ( Z( I1,8),RH0( 1,8),PT3( I ),pI=1,M) 44 F-ORMAT( '0' 5X,'Z',1OX, 'RHO',1l2X, 'T3'/(' ',3(F12.6,y2X))) GO TO 37 99 0 WR ITE (6,r991) 9 91 FnlRMAT( 10"-** ERROR IN DATA' 999 CALL SYSTEM ENn 95

013630-9-T A. 3 SUBROUTINE DATA (IN, MX, MY, ZEP, RHOEP, THETA, B, VOL) This subroutine is called once for each segment of the profile. From the input specification for the segment, DATA computes the (z, p) coordinates of the necessary sampling points on the profile, the quantities cos a and sin a at these points and the incremental volume of the segment. Arguments: IN Type key for segment. MX Total number of cells in segments to the left. MY MX + (number of cells in this segment). ZEP z- coordinate end points of segment: ZEP (1)" z ZEP(2) z2 RHOEP p- coordinate end points of segment: RHOEP(1) = RHOEP(2)" 2 THETA Angle (in radians) subtended by a circular arc at its center. B Fractional exclusion,. VOL Incremental volume of segment. Comments: Stored in COMMON are the arrays RHO(80, 9), Z(80, 9), ARC (80), C (80, 9) and S (80, 9) which contain the numbers computed by DATA. For the Ith cell, the subscripts (I, J) correspond to the points si of (229) when 1 < J < 6. For J - 7, 8, 9, the subscripts (I, J) refer to the points i- i and i+ respectively of (211). 96

01363 0-9-T SUB~ROUTINE DATA(I1NMXMYZEPRHOEPTHETABVOL) DIMFNSinN ZEP(2),RH0EP(2) COMMON RHO(80,9),Z(8099),ARC(80),C(80,9bS(80,9) DATA STFP/.3872988/ MXPl =MX+l EN=FLOA T( MY-MX) IF(B.NE. 1) SUI3STP.*5*(1.0-B)*STEP IF( IN-?2)19293 cc=-i'o GO0 TO 1 0 2 cC=1.0 10 S T?2= SIN (T H FT A /?2.O A=ZEP(2 )-7EP( 1) RAD=o.5*SORT ((RHOEP (1)-RHOEP (2)) **2+A*A )/ST2 DD)=A/ABS (A) T=CC~nDD*COS (THE TAt?2*) 1ST2 ZCNT=0. 5*~(ZEP (1)+ZEP( 2)+T*( RHOEP (1)-RHOEP (2))) RHOCNT=0.5*,(RHOEP(1)+RHOEP(2)+T*A) tJ2=ZEP(2 )-ZCNT U..l=ZEP(l1)-ZCNT VOL= 3. 141 593*ABS (A* (RHOCNT**2+RAD*RAD- ( U2**2+Ul*U.2+Ul**2) /3.0) 3 -CC*RHOCNT*(U)2*(RHOEP(2)-RHOCNT)-Ul*(RHOEP(l)-RH(3CNT) +RAD*RAD 3 *DD"-THETA)) BETA=CC*DD*THETA/EN THETI=ATAN2(RHOEP(l1)-RHOCNTZEP(1)-ZCNT) U=ABS (BETA*RAD) B3=S TEP*kRETA DO 902 I=MXP19MY PHI=THETI+( I-MX-~.5)*BETA IF(B.EQ. 1.0) GO TO 1905 DO 1902 J=192 4NG=PHI+.5*(J-l.5)*BETA*(1.O+B) DO 1903 L=1,3 PSI=ANG,+(L-2 )*SUJBSTP*BETA M=L+3*(J-1) C(1I M)=-CC*SIN(PSI) 5(1,M)=CC*COS(PSI) Z ( I vM ) =ZCNT+RAD*CC*S ( I vM) 1903 RHO( 1,M)=RHOCNT-CC*RAD*C(1,M) 1902 CONTINUE 1905 DO 903 J=7,9 ANG=PHI+( J-8 )*fB3 C( I v J) =-CC*S I N(ANG) S( I J ) =CC*COS ( ANG) Z( IJ)=ZCNT+RAD*CC*S(1IJ) 903 RHO(IvJ)=RHOCNT-'CC*RAD*C(1vJ) 902 ARC(I)=U RETURN 97

013630-9-T 3 L)X=( LtY(2 )-ZFP(1) )/EN DY= (RHnEP (2) -RHOEP(1)) /EN U =SQRT(DX*DX+DY*DY) SI =DY/U) CT =DX/U nO 917 I=MXPlMY PHI=FLOAT( I-MX)-.5 IF( F3.EO 1.0) GO TO 1800 DO 1802 J=192 ANG=PHI+.5*(J-l.5)*(1*0+B) 00 1803 L=193 PSI=ANG+(L-2 )*SUBSTP Z (1,M)=ZEP( 1)+PSI*DX RHO( IM)=RHOEP( 1)+PSI*DY S (I,M ) =5 1803 C(lM)=CI 1802 CONTINUE 1800 00 913 J=7,9 ANG=PHI +( J-8)*STEP Z (I,J)=ZEP( 1)+ANG*DX RHO(I1,J)=RHOEP(1)+ANG*DY C( IJ)=CI 913 S(IPJ)=SI 917 ARC(I)=U VOL=1.047198*(ZEP(2)-ZEP(l))*(RH0EP(1)**2+RH0EP(1h*RH0EP(2)+ 8 RHOEP(2)**2) RETURN E N D_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 98

013630-9-T A. 4 SUBROUTINE INTEG (V, NSEG, NUMPTS, SUM) INTEG numerically integrates quadratic interpolating polynomials approximating the data on each segment of the profile. When the profile is composed of several segments, no interpolation is performed across segment boundaries. Hence, the integration is accurate even for disconnected segments, e. g. the circular arcs of two spheres. Arguments: V Real vector of function values, ordered as the cells. NSEG Total number of segments in the profile. NUMPTS Integer array containing in NUMPTS (I) the number of cells on the Ith segment: I a 1, NSEG. SUM Integral of V across the profile. Comments: Stored in COMMON are the ar lengths ARC (I), I 1..., N required to compute the integral. 99

01 3630-9-T S IiFAR0IJ IT I N I NTFG ( Vp NS FG, N(IMP TS,SUJM (AIMMON RHO(BA0,9 ) Ia ( 80t9 ) A RC,(8O),C(80,9),tS(80,9 DTlMFNS iO~N V( 80),NUIMPTS( 15) SU)M=O0. 0 J AC, C= 1 nnl 3000 I=lNSFG T= ARC, ( J ACC ) L =NIMP TS ( I ) N=I.-JACC-l SU)M=S1IM+T*(0.625*(V(JACC)+V(N))h.125*(V(JACC+1)+V(N-1))) IF(L/2 *NF. (L+1)/2) GO TO 3001 SI,)M= SUM+T* ( 0.66 666 67*V( N-I )-0. 08333333*V( N-2 )+0.4 166667*V(N)) 3001 I-ml=N-l J Ln0= JA CC +1 no 3002 J=J)LnOL M I,? 3~002 SUM=SUM+0.3333333*T*(V(J-1)+4.0*V(J)+V(J+1)) 3000 JACC=JACC+L RETURN END 100

013630-9-T A. 5 SUBROUTINES DECOMP (A. N) AND SOLVE (A, X. B. N) Used together, DECOMP and SOLVE solve the linear system AX " B. DECOMP performs a L-U decomposition of the N x N matrix A and SOLVE performs back-substitution. These routines are adapted from Forsythe and Moler ( 1967, pp. 68 - 69 ). SIUBROUITINE DECOMP(IJLN) DIMENSION UL(80,80) COMMON /SOL/IPS(80) nn 5 I=1,N 5 IPS(I)=I NM1=N-1 DO 16 K=1,NM1 BIF=O.0 DO 11 I=K,N IP= I PS (I) IF(ABS(JL(IP,K)).LE. BIF) GO TO 11 BIF=ABS( UL(IPK )) IDXPIV=I 11 CONTINUE IF(IDXPIV.EQ. K) GO TO 15 J=IPS(K) IPS(K)=IPS(IDXPIV) IPS( IDXPIV)=J 15 KP=IPS(K) PIVOT=()L (KP,K) KP1=K+1 DO 16 I=KPI,N IP=IPS( I) EM=-)L ( IP, K)/P IVOT lJL( IP,K)=EM DO 16 J=KP1,N L ( IP,J )=UL ( IP,J)+EM*UL (KP,J) 16 CONTINUE RETURN END 101

01 3630-'9 —T SIIBRn0JT T NEsnLVE (IL, X,R,N) DTME-NSifN UtL(80,80),R(80),X (80) CfIMMON /SOL/IPS(80) NPI=N+l IP= I P S(I1 X( 1) )::( IP) 00 2 T=2,N I P= IPS (I) I ML= I-I SIJM=0. DC) 1 j=i,im1 1 SU)M=SU)M+UL(IPJ)*X(J) 2 X(T)=(ITP)-S(JM I P= IPS (N) X(N)=X(N)/U)L(IP,tN) DCn 4 IBACK=?,N I=NPI-IB~ACK I P= IPS (I) 1 P1=I+1 S(UM= 0.0 00 3 J=IPLN 3 SLJM=SUM+U)L(IPJ)*X(J) 4 X( I )(X(lI)-SUM)/tJL( IPI) RETURN END 102

013630-9-T A. 6 SUBROUTINE ELLI (Ml, K, E, KPR, KEY) This computes the elliptic integrals K (m) and E (m) and the derivative '(in) from their power series approximations (see Section 4.2). Arguments: Ml Real, the quantity (1 ain). K Real, K(m). E Real, E(m). KPR Real, K'(m). KEY Integer: 0 Compute K, E and KPR; 1 Compute KI E but omit KPR. SUBRRoTINE ELLI(ML,KtEKPRKEV) REAL M KtKPR T=-ALOG C MI) K=l.3 294+.5*T+M1*(9.666344E-2+.l249859*T+Ml*(3.59OO92E-2 5 +-6,880249E-2*cT+Ml*(3.742564E —2+3.328355E-2*T+Ml*(1.451196E.-2 5 +4.41787E-3*T)))) F=1O+Ml*(.4432514+,2499837*T+Ml*(6.26O601E-2+9.200l8E-2*T+Ml*( A 4,757384E-2+4.O69698E-2*T+Ml*(l.736506E-?+5.264496E-3*T)))) IF(KEY.EQ. 1) RETURN KPR=.5/Ml + 2.8'3225E-2 -.1249859*T + Ml*(-2.999362E-3-.137605*T 6 +Ml*(-7.899336E.-2 - 9.985066E-2*T * Ml*(-5.362998E-2 5 1767148E-2 * T ))) RE TURN END 103

013630-9-T A.7 SUBROUTINE SETUP (KEYPll, KEYMll, I, J, L, APIll, API33, AMIll. IJ) This is essential in computing the linear systems. Specifically SETUP, after calling ELLI, computes the quantities API 11 (K1 of eq. 182), API33 (K0 ofeq. 178), f21 (eq. 196) and f2 (eq. 197). The quantities Sfl and 2 are used to compute f(i, J) (AM11) of eq. (227). Arguments: KEYP11 0 when computing API11 and API33, else 1. KEYM 11 0 when computing AMI 11, else 1. I Subscript of observer (unprimed) cell. J Subscript of remote (primed) cell. L Index of the point within remote cell for which the kernels are to be computed (see DATA, Comments). API11, API33 Described above AMI11 IJ 0: use last value of Ml in kernel computations; 1: compute new Ml. 104

01363 0-9-'T StJFRRlhTINE SETIJP(KEYP~lKEYM~ll, JLAPIIIAP133,AMI11,IJ) COMMflN RHn(8o,9),Z (80,9),ARC(80),C(80,9),S(80,9) PEAL MtMItKKPR 70=7 (JvL )-~Z( 1,8) R=RHfl( i,8) RP=RHl(,J,L) IF(IJ.EQ. 0) GO0 TO 115 RRP=R*RP Al=RRP+RRP A2=R *R+R P*RP+ZOD*ZDf Ml=( A2-Al )/(A2+A1) CALL ELL I(MlKrEtK PRKEYM11I) AO=M/RRP Al=SO)RT( AO) A2=M+M 43=2.-M IF (KEYPll-l)l13,114, 114 113 APT111=Al*(A3*K-E-E) /M 1 14 A I=AI*AO IF (KEYM11-1) 115,116,116 115 AO=C(JL) nmi —Al*(.25*K-A3*KPR) flM2=Al*(E-A3*((A3+M)*K -A2*A3*KPR))/(M*M) AMI111R*AO*0M2+(ZD*S(JL)-RP*AO)*OM1 116 RETURN END 105

UNCLASSIFTED(I., irt. DOCUMENT CONTROL DATA.- R & D (Sr %,l I 'ct NC ' If -0...' Of fIl I fI h If. I'.r I tint 0. b I, III Ft i -iI,i It I mI hv 'f., for FI win VIt vnI/nevri,l I r eot) rI~ rI v c N,.I I 0I 1. IN A IIN C A C I IV I I IC r, ri' iI C- mithfi..) T7 0 4L 11 RT SEC U 141 1Y C LAS SI FCAT ION heUn1ivers1 —:ity of Milg 1'i Rii~tionl L;1l)oU-,tolry UNCLASSIFIED 2216 Spce Reserch Blr.,North C.i inpuz;Szb Ann Arbor i ihi, 18105 --- 3. H LPO RI T I IL THE NUMERICAL SOLUTION OF LOW FREQUENCY SCATTERING PROBLEMS 4.0 DE SC RI P TIVE N 0T E S (T7pe of report flndlrnclh vI ~ve datom) Scientific Inte rim S. AU THOR(S) (Forat name, mlddlu' uiitinl, last #,m)IMThoniwis B.A. Senior David J. Ahig-ren F EbORTuAry 197i7*. TOTAL NO. OF PAGES 7b. NO. OF RCFS $8P. CONTRACT OR GRANT NO. 9o. ORIGINATOR'S RLPORT NUMBE4RrR51 F19628-68-C-0071 013630-9-T.b. PROJECT NO. Project, Task, Work Unit Nos. Scientific Report No. 12 5635-02-01 9h. O THCR RCPOR T NO(S) (Any other,smomcra thot nuty hi' &Asl4neu d. DoD Element 61102F ACL7-06 DoD Subelement 681305 10 IST RIBDUTION STATEMENT A - Approved for public L'elease; distribution unlimited. II. SUPPLEMENTARY NOTES 12. SPONSORING MILITARY ACTIVITY TECH THERAir Force Cambridge Research Laboratories (LZ: TECHOTHERL. G. Hanscom Field II -.3 ABSTRACT Bedford, INa'ssa'Chusetts 01730 I'm. jfto I "^- I The low frequency sca-tterinhg of electromagnetic and acoustic waves by roln-tiol:-Illy symmetric bodies is considered. By concentrating on certain quantities such as the normalised component of the induced electric and maignetic dipole moments., it. is shown how the first one or two terms in the far zone scattered fields can be expressed in terms of quantities which are frunctions only of the geometry of the body. Each of these is the weighted integral of an elementary potential function which can be found by solving an integral equation. A computer program has been written to solve the aippr-opriaite equations by the moment mnethod., and for C.1cUIctil~Itg the dipole moments, the " I ci Giv~I capacity, and a further quantity related to the capacity. The program is described and related data are presented. %W r M-Nown. I rN M F O R m d a -v o% - - - mmum L)LJ I NOV r 1 4/7.i UNCLASSIFIED ___

See IIrit I C1.u I t Low requency Scattering. Electromag-netic Acoustic IRotntional Symmetry Computer Prog-ram. Numerical Data I I I I i I I i I T TCTIAT L, ELEL___