~EE]r-GAP T rans. RL-5 013630 NAPS 01891 7-6-72 RL-553 = RL-553 COMPUTER PROGRAM: RAYLEIGH SCATTERING T.B.A. Senior and D. J. Ahlgren The Radiation Laboratory Department of Electrical and Computer Engineering The University of Michigan Ann Arbor, Michigan 1. INTRODUCTION This program solves integral equations to compute the three independent polarizability tensor elements Pll, P33 and Mll, together with the capacitance C/e and a quantity y related to this for rotationally symmetric, perfectly conducting shapes composed of either one or two bodies. The tensor elements are as given in Keller et al (1972), but see also Senior and Ahlgren (1972) where the electromagnetic capacity C and the quantity y are defined; E is the permittivity of the surrounding medium. The numerical procedures are discussed and the various subroutines presented along with a flow chart of the program. 1.1 BODY DESCRIPTION AND ITS VOLUME The body B is required to be rotationally-symmetric about the z-axis of the cylindrical coordinate system. The profile of the body in the p, z plane is given by the function p = p (z), and is constructed of segments which may be linear arcs or circular arcs of either concavity. Distance along the profile is measured by the variable s while the tangent at any point forms the angle a with the z-axis; -7r/2 <a < 7r/2. The volume of rotation of a linear arc (type 3) segment whose endpoints in the (z, p) plane are (zl, p ), (z2, p2) is z =V+ zTrp (z)dz 3 i2 1 P2 + * (1

For a circular arc which subtends an angle 0 at its center ( z0, p ) and is of radius a, define a quantity ~ where j= ( (2) 1 segment concave up (type 2). Then, 6V0 = ( Z l) + a - (u1U2+ uu2 -0 2 1 3 + 2 2 (3) P p 2(P2 -Po ul(PlP )+da20 0 2 2 0 1 1 0 where z2-z1 d = 1-_ (4) and u2 = z2- z u1= z1- z0 1.2 THE INTEGRAL EQUATIONS Given the profile p p (z) and the volume V0, quantities m and m1 are defined as * 4pp' m = m(z, p, z', p') = 4ppt2 (5) (p+ ')2 + (z - z')2 and m = 1 - m (6) where (z, p) and (z', p') are distinct points on the profile p(z), the unprimed variables denoting an observation point and the primed variables a remote point. Note that 0 < m, m < 1. Denote by K(m E(m) respectively, the complete elliptic integrals of the1 -Denote by K(m.), E (m) respectively, the complete elliptic integrals of the 2

dm below: i. To compute 11/VO, solve s p' K1T1(s')ds'= 2rp (7) 0 where K1 K (z, p, z, pt 2 {(1 - ) K(m) -E(m) '/2 ~2 (mpp') (8) and compute S p p pl p2T (S) ds (9) ii. To compute 33/VO, C/, and y, solve f p'K0T2 (s)ds'= 2z, (10) and p'K0T3 (s')ds' = 27r (11) n(11) where K 0 = K(m) (12) sC/T(s)ds C/e= 27r| p T3(s) ds, (13) Jo 3

55: 2= - zpT3(s)ds (14) C 3 and P33 2 C2 33 2 zpT(s)ds- e V (15) V0o. (15) iii. If and only if B consists of two separate (electrically unconnected) closed parts B1 and B2, solve s 1 2 ir for points on B p'K T (1)(s) ds' (16) 0 for points on B2 compute 2 V-2 (z + y) pT (1) (ds 6P33 s (17) V (1) 2 ()(17) JpT~) (s)ds- 2 < pT3(s)ds where the symbol (1) below the integral sign means that the integration is carried out over the profile of B1 alone. Then the tensor element P33 for the disjoint shape is P33 P33 6P33 -- - - 3 (18) 0 VV V0 iv. Solve 4

,J U p cos a' + f(z'- z) sinc' - p' cos a' = 7 TV4(s)-2p i1 p'ds' (19) where = 1 ) K'(m) - K(m) }2 4 (20) m Ci/2 "2 P 2m (I ) 1 2 KT(m) 2 ) - ( - )K(m)+ E (m) (21) and the bar across the integral sign denotes the Cauchy principal value. Compute 011 v 0 s f pV4(s)cosads Vo 0 (22) 5

2. NUMERICAL SOLUTION OF THE INTEGRAL EQUATIONS The numerical procedures involved in finding P33/V, /V 0 11/V0 and where appropriate, 33/ V are quite similar to those required for P P 11 /V, and it is therefore sufficient to give full details only for 11 /V 2.1 Pll/V COMPUTATION 0 The primary task is the solution of the integral equation (7) for the function T (s) and this entails the determination of a sequence of values (i) 1 T1, 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 1 N cells C. of arc length As. and midpoints s. corresponding to the coordinates (Pi., z) Within each cell we also define the points s. and Si+ where S. = S. - a As i- i 0 i (23) s s - a As i + i i with the restriction 0 < < i (i) By assuming that T (s) has the constant value T1 over the ith cell, 1 6

the integral on the left hand side of (7) 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. T f pt Kds'+ T1 f ptK1 ds? +....+T N p' K ds 2 rp C1 C2 CN (24) i 1, 2,.d.., N. Hence, the system to be solved is At = b 1 where t is a column vector with elements t T (i) t = T ( li 1 A is a square matrix with elements (25) (26) ai. 1J - J p'Kld' C. J i j - 1, 2,...., N, (27) and b is a row vector with elements b.J 27rp, J J j = 1, 2,...., N. (28) Increasing the complexity of the quadrature technique used to evaluate the integrals I will generally improve the accuracy but will almost C. J 7

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 i j): a.. pK1 (i, j) As. (29) ij _ 1 -Kl(i j-)+ pj+ + Wo()}+) pjiKl(i, j AS (30) respectively, where the subscripts j - and j + correspond to the points s. and sj of eqs. (23), and K (i, j) is the kernel defined in eq. (8) and evaluated at the points (Pi. zi) (p., z.). By requiring a0 < 2 we ensure that the sampling points sj and s 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. 20 and 21). When -1 2 (31) 0 2 ( eq. (30) reduces to the three-point Gaussian formula for which 4 5 WO 9 w1 18 (32) With this choice of w0 and w1, the advantages of eq. (30) vis-a-vis eq. (29) were now determined by computing l/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 8

G 5.0 5.0 - 4. 0 percent 3.(0 error 2.0 1.0 0 4.0 C.P.U. time 3( (sec.) 2.0 1.0 0 0 10 20 N 30 40 Fig. 1: ercent errorandC.P.U. time of P11/V calculation for a sphere: T denotes trapezoidal rule computation and G denotes three-point Gaussian. 9

553 the trapezoidal method, though the accuracies of both are severely degraded if N is too small (N < 5). Since the Gaussian scheme with N = 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. ij In summary, the integral equationi(7) is solved by conversion to the matrix system (25) in which r, 5 r j.)+ K (i, +).+ 4 aijL 18 P li - j li 9 PjK + asj i, j 1, 2,...., N; i j (33) 16p. a [n i -l\As ' i 1, 2,....,:N. ai i / nj ii As (i) P Having determined the sampled values T T (s), 11 is computed from eq. (9) by integration over each segment of the profile using a second order integration procedure (subroutine INTEG, described in Section 3.4). 2.2 33 / V0 COMPUTATION The point sampling method of solution of the integral equations (10) (i) T (si) and (i) T and (11) requires us to find the sequences T2 T( s) and 3 T s i = 1, 2,...., N, from these equations. To determine the T, choose a 0 w0 and wl in accordance with (31) and (32) and thence solve the matrix 10

system At b where 2 t T (i) 2 i 2 i - 1, 2,...., N (34) b 2r z i i and i, j 1, 2,...., N; i 7 j, (35) 3 3 aii = I [n (6p i 3 + /si ' i -1, 2,...., N. where the elements aij are again given by (35), but = T(i) 3i 3 i - 1, 2,....,N.(36) b 2 i The quantities C/e, and P33/V! defined in eqs. (13), (14) and (15) respectively are computed using the same second order integration procedure employed in calculating 11/V If the body profile consists of two discrete parts, it is also necessary to solve the integral equation (16). The corresponding matrix system is 11

553 almost identical to that in (36), and from the sampled values T3 /V. 33/V0 0'" and T ( i), 33/V0 (seeeq. 17) is computed and, hence, 2.3 11/V COMPUTATION The basic approach is similar to the above in spite of the more complicated integral equation (19) that must now be solved. The matrix equation for the sampledcvalues V (s.) V4( is Av 4 b where (i) 4 i 4 ii N i s 1, 2,....,N (37) b - 27rp i i and I5 a.- 18 ij 1 8 p f(i, j-) + pj+f(i,j+) + 4 f(i,j) As. i, 1, 2,.... (38) s +- As. 1 2 7T - T i 1 S - 2 A S Op' f(s') dst i 1, 2,...., N in which f(i, cos (i, j) +osj 2 ) + j zj- zi) sinc j- pjcos aj (i, j)..., N; i j.(39) i, j - 1, 2,. 12

55 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 i). As an approximation to thisprincipal value, we remove from the cell C. a slice defined by the interval 1 1 ( - /3 Asi s. + - /3 As ) where 3, 0 < / < 1, is the fractional i i i 2 i exclusion; /13 1 implies no exclusions, i. e. that the principal value is not taken. We now have S. - A. S.+ - A. i 2 i i 2 i a - p' f(s')ds' - p'f(s) ds' nJ1 - s -- As + - As i 2 i i 2 i (40) and these integrals are also computed using three-point Gaussian quadrature. Defining 1 Si2 - -- (1+/) As. 12 1 4 1 1 si3 i2 2 1O i 1 (41) s,5 S + S (1+/) A s 1 si4 Si5 2 0 3) Ai 1 Si si5 + 2 % (1- 1) A S si6 i5 2 ai 13

55U we obtain aii 2 (1- ) f(i,i)+p.3f(i,i3)+p f(i i) 2i1i Pil ' 1 i 3 '3 i4 J 4 +Pi6+ f(i, i2)+ 5 f(i, i) As (42) Equations (37) through (42) completely describe a system of N linear equations in N unknowns Vi, i = 1, 2,...., N. Their solution and subsequent integration of the V4i according to eq. (22) yield Mll/V 4i 11/0 14

55t 3. THE COMPUTER PROGRAM The program computes 11/V CP P7 33 /V A1/0, C., ly, 3V M /V and, where appropriate, 33 /V, and consists of a main 0 0' program and six subroutines. 3.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. Control Card Columns Description 1 The number (1 or 2) of bodies. 3-4 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. 6-7 Same as columns 3 - 4, but for body to the right. 9 A printing key: P 1: print T from 33/V computation. 0 or blank: do not print T. 3 15

11 A computation key (0, blank or 1) 1: suppresses computation of 11 / V0 C/, P 33 /V 13 A computation key (0, blank or 1) 1: suppresses computation of 11 /V. 21- 30 A real number: the fractional exclusion 13. If these columns are blank, j3 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 1 and z2 of the segment:. 16

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 > z. 17

553 3.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 I>DA DATA -- < -- I 18

553 F Al '~l1 (8 n0 Ap~8,P)Av (,8)X(8)-1 71 ~-P (2 R HFIEP (2 S T3 (T 3(0,~,U( 'I N\TFG-F-P NIJMPTS( 15),PUJo/I I+/ k8L I ' N 2 CPviON Ho( 80, ) 7 (80,9 )A, AkC 80 ) H (0,9) S (80) ~L/ J TS( 0 PA T A T1,1 V~I PIP,PI /'-',4 8153 * 141-59-3, *44/t4-4 4 * 2 7 77 117/ 37 FAP(5R4, P=99)N PPN-S I -vf\i S2,vI P kIN TKEFY P I K F Y1 P1,FkI 34 F F1R NA T( 111 (I 23.X I, -t( I,X )6, X I F 12I 7 kIP. I T F( 6,t4)N8OC,NIS1, 4 F 'P M A T i~ EG I NJNIViG OF 1 ATA S FT:*' / '0l' 5X, '4 CIDILS6,X, ' 2 '-' -,X, 'S FENlPNT S: I / 'I,X 'I OY 19;,Xt, '= '1) I F (NIP R.l0L F' 0.0 Rk. R II H.6 T 2 ) (0 10 000(4 NI S F (; = N 1~ I + N\ 2, FOR M A T( SX, ' ROY ti' 5 X, 1 12 I F ( N! RF2) 0 F * 1) (f T r 10 IF(NS11.LE-0F0. NS2.LE. 0) 60 TO- (90( 1 T E( P1S FG.L. 0 PRP. INISEFS.61. 15)G 60 10 99q0 352 FRMATTF( ~ ' ISX 'P PINT tK FY' P,3XKF 0 ''12/ S,3OvP yfU~:I, 5 'RN ',5X 'CX PE R M I'= K F I? C M Y: IF( FE.L F..0. F R.61. i. )FR=.001. TF(KEYWV 11.. 0) WPITF(6,51)FR, I F (KEFYPI 1.*E 0 *A V0. KFY NHI *NF. 0) GO TO0 9 9 \I 0=0,. P 12 F AIT( 2, I X, I II X,9AI,4 X',q5F1 07) I FU\ J VP TS( I).L. 0 P,0 T ITY P -LF. 0 0k, IITY P.61.T 3)6 10 99 0 IF (1 SiTG N EQ. RL T) IJ G NIPi_ US5 1I FO PMIA, '(I OSL-;GPAENJT H' I J I2 I / ', 5Xv, 'GEL LS ',7X, '= ' 1 /' ' -, X, 4 TY P K FY',X,=1/''5 X V /OL I I AE SEFN SE F iL\4/,X, 8 ' 7CP f II PJ IN AT' F FND F POIN lT S =( F,1.7, ' F127, 5 i x, I F(ITY~.MEF 3) WiRI T 6,1 4) THETA 14 FOMA T ( X T 'HEFT[A (PE) F = I, FO5 I F(KFY M, _I F FFPP,1 ER=1. THEF T A =P 1 T HF A/ 10. viM \1 rP T S ( I IF: I.I.NS1 N Cl=.N C1~I NitJ MP T S (I IJF(.61 80)( TV 99 0 C',ALL[ P)A-'A( (I TYP,N,If, 7 FP, PHiDER,THEfTA FR, VO-ltN IF ( T 'TYP.N\]F 3 IANP 5IIG M FC-0 vIM). VO0L I k\C=-V\IIL I C 11 \0,= VO~'/OF) INyC, tJQI TF (64 52 ) V 52 __ P, M A T(C ' O P't UTFF)P, RSU(JLT S: / 5,X, 'VL iJM I', 4>: 6 = ', Fl 5) 19

553 Fil 2. NIl, TN 1P, H fMj T~ 1_ (P, i \j )l =0. i. 09 T F (K FYP I 1-1)1 1, I, I1.,I A P 1 1 (L N')=0 *0 AP3 NI 1L)=. A P3~3(IN)=0. 1 2 If1 NPX(2)=L 20

5 5 - T 104 X=, IIp p II. C AL-L ~F 11 JP KFYPIi IKFYMli J-?,P~AIi A13AJ1~ 10 ~ I( I, 1 2)=API I(11912 )~APi3. 1TI I 7 A`133I( Ii, 1 2 )=IPI3( I11,1I )APMI33T I I04 CON TI NI H.iF CALL I S FTIJ P(K FY PlI IKFYM~Ii,,I AP I 1, IA PI 3 3 A Iv 4)II TIF KEFY P 1 1. 0 8,209, 209 Ai 1 NL )A L ~WI4A P II (NL ) +i JT L A P11 ( L N, N) N~ (Id ~AP ii(_N)+VTN N0 A T 33 AP313 (3 I IN) AN I (N AP33 ( L N) + Iw — Tr 20(9 TIF(K FY MI- 2 P)213 p3 20 AI11(NL )A(1AA flq(NL)-AMl1l 'L) AN! 1 ( L,Ni =ANO I ( I0 1NM 1 ( LN,N)AJ11W 0T CF Tn T F (KEFY P II-1)2.3,8 4,t8 4 3:A LnCZ 1 6 T N/AN! 84 IFKYl- 85~ 15 F(F R *F ". 10) (I1 T U 3 O 26 =16 1:*5(1*0- F R ~A \ Ai1(N\)P I -U(W0JA (RH (N ~T 3 (? +R H(1N5 ~T3( 5 3- +W~ I ( HJ,1 'T 3 (1 + RH O (N, 3 T 3 (3) RH (N, 4) S T 3 (4) 4 SH- I,4 ST 3(6) 2 CF-1INT INI 21

5 5 3 21 CAl I- FV l I fP( P1j 1,r CAL I ALV(A i,Xk, F) 22 2 2 XI R PHO (I 8) * X(I CALL lT F x NUE~ MfvPTS,P 1)I p I I p I p I /,- V 0 2 4 lF(KFYMT1i-1 ),28,? 5 C A LL F C n M P ( -41 )_ CALL SOL-VF( A f~I 1 X Pf I-) 2 6 J =I?,M CALL R TNFG( I, F8,XJMPTI II11= P I 8 11/ v 2 8 I F (VK F YP I 1-1 ) 32,45,p4 5 32 Ffl 2 9 1 =1?,A -29 I(J = Ti.OR nPI7 ( I,8 CA! L -) FC ~(IVPP 33,[VI CAI STLV ( AP3 XI, X, I LZ (I N, 8) N RH (I'SN M 8 ) SX (ID F) f 35 1=1,I X( I )RHO( I,8VT 3( I) 35 8~ (1 I=7 (1 I 8 )- X (I) CA LL I N TEC x G, X- N! SF S ', t JM PTS C A P CA P= Twi f-1P {CA P CA4L L I N TFGC (8,NS,NoivipTrs,c.; A v G A T = - UIw P I c; A 8 C A P P 33 TWOJf PI P 3 3- CAP (GA V G AWfA I F(NRnF) - 1) 39,?39,';4 5 4 DO (36 1 =1,8M IF0 - ( I-T,,.lCi) 305 3I,0 306 PI=. 36 C O NITINHE__ _ _ _ _ _ _ _ _ 22

553 LV. T I,T Xn 0- I T 1, T NX I)= TI G, fX(N I) NitINlPT TN (T LI NJ TLO X, NS F( G S NUMPTS,1i FD FL TA P= (TW O1P I 1/ V ) u+AM T I. ) T~ (\ I-T W P T ~T TL /CIP =03 3+ O)F L TA P 39 V'PITE( 6 4 0) CAP,?G A vPi,3 3 4 0 FFPAA(' 5X,'CAPACITANCE =',FIO.5/' I,5Xp 'DAMMA' 1,7Xt v=,F 1O5/ 8,X'P/V 7X,?=,F1O5' 57'P3V 7X, ' JF(NfD *F I )PTE(6,309)DFLTAP IJ 309C F R -!'IAt\T(' X ',5),'DFT P 3 3/\/V X ''F I 0.5/,X LiS J T P 32/V6 FlO5 45 IF(KFYM1iA-1)4~,r337-t337 42 WRI TF(6,43)MI1I 33 F ( I PPV;',T IE *AN O. KEY PlII.EO C) RWPIT F((44) 6 ( 7 (I I 8),tRHO( 18),T 3 ( I )= I= -,NI 44 F OPMA T 'C,5X '7'I X OX, 'HU',1X T3/ ',3F 12.6,2 X) CD T O 37 990 T TjTE( 6,9 91) 99 F fPTA! A T( n F PRHP I N D)AT A F L SY DFM 23

553 3.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 )= ZEP (2) = z. RHOEP p- coordinate end points of segment: RHOEP (1) = P1, RHOEP (2) = p THETA Angle (in radians ) subtended by a circular arc at its center. B Fractional exclusion, f. 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 s.. of (41) when 1 < J < 6. For J = 7, 8, 9, the subscripts (I, J) refer to the points s., s. and s respectively of (23). 1- 1 i+E 24

5 5 1 SI,)RR OU TINE DA TA( IN,0 X9rY, Z F.PRHCOFEPPTHFTAI, hVOL) DIMFNSiCnp ZEP(2),wRHDFP(2) COMMON RHO(80,9),Z(80,9),ARC(80)) C(80,99),S(20,9) DA TA S TF F'38 72988/ MXP 1.=MX+l FN=FLn)ATBMY-MX) IF(3.NE. 1) SIJEBSTP=.5; I10-B "'STEP IF(lIN-2)i,-2,p3 cc=-1I.0 GO TO 10 2CC=1I.0 IV) S T2 =S I Ni(T H F T4A/2.0) A=ZEP(2 )-ZEP( 1) RAD)=O0*5"SQRT ((RHOEP (1.)-RHnEP (2))*-"2+A"-A )/ST2 DD)=A/ABS (A) T=CC`0D*C0S ( THETA/2,*) /ST2 ZCINT=0.5'-IZEP 1 )+ZEP( 2)+T-"(RHOFP(lI)-RHC)EP (2))) RHP, CN T=0. 5,-"( R HCE P( 1 ) +RHCf) EP ( 2)+T4*A) tJ2=ZEP (2 ) —ZCNT U1l=ZFP( I)-7CNT V L = 3. 14+1 5 93 'A BS (A'( R HO0C NT *2+ RA D R A D- ( t2 *2 +U 1 U2+UI 2) /3.0) 3 ICCRH)CT 91i2 IV,2 H)CT-1(RCE(1-HiNT RDRAD 3 '4DO"TETA) KEF TA = CC4- )D-",T H T-A/E N I-HI' T1=ATAN2(RHC)FP(l)-RHO)CNTZEP(l) —ZCNT) UJ=ABS ( KFTA*RAD) K'3 = S TE P -:- E TA DO 902 l=MXPINlY PHI=THFTl+(I-MX —.5)`BFTA IF(R EO. 1.0) GC) TO 1905 DC) 1902 J=l,2 [ANG=PHI1+. 5" (,J-1,*5),BETA-* (1.O+B) DC) 1903 L=1,3 PSI=ANG+(L —2)' SUJKST P4BFTA m = L +3 —( J- I C(1,vM) =-CC —',SIN( PSI) s(1 m) =c~C:*cC)s( PsJ Z I M =7.CNT +R AD* C CS ( I M) 1 903 RPHO (IMl)=RHOC)NT-CC-R AD,-C(IM) 19 02 CC)NTI NUEF 19l)0 5 DC) 903 J=7,P9 AN G=PHI-I+ ( J -8 )', 3 C(,=-CC"S I N (ANG) S(I tJ )=CC*-,CDS (ANG) Z IJ)=ZCNT+RAD-CC-',S (I,J) 9 03 RH(,) (IJ ) =RHC)CNT —CC RAD"4C(IJ 9 02 ARC( I)=U PRE TUR N 25

5 5 3 3~ JX~T7EP 2-7 E (i ) EN D)Y P (H FI FP (2 R RHO0E P (1) / EN P =S.kT (D O'DX + DWDY)0 S I DY/ U C I = DX/ t I WI 917 I=MXP1,PMY PH I =F:LOIAT( IL-MX ) -,5 IF( B FEO. 1.0) GO IC) 1-800 DO 1802 J=1,2 r(1 1803 L=1,3 M=+3~J P S i= A11G,+(L-2 )S UBRST P Z (1,M)=ZEP( 1)+PSI"'DX P\HO (I v t-' ) =RHC1EP (1I) +PS I "DY S (I,M ) =SI1 1803 C(IM~)=CI 18 02 CON TI NU E 1,800 OF) 91:3 3=7,t9 AN(;=PHI+(J3-8 )*STEP 7 (I,p3 )=7EP (1) +ANG*DX RHn(i ~J )=RH0EOP (I) +ANG,'DY C (I, 3)=C I 913 S(IJ)=SI 917 ARC (I)=U VOL=1,0471(98,(ZEP(2)-ZEP(l))-"(RHOEP(l)""2+RHOEP(1)-'RHOEP(2)+ 8 R H [)EP (2 )r42) R ETIJR\IJ E N D E_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 26

3.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 - 1, NSEG. SUM Integral of V across the profile. Comments: Stored in COMMON are the arc lengths ARC (I), I 1,...., N required to compute the integral. 27

5bt IFA~l T HO I NT,9 ) 7 F G( ). I '\J (T 80 )~, ' 0 &(i 'ARC, Ho 0 00 300) I =1NS F G T APRC, ( I I%[A ) I ~=1 46JC -1 I F(-/ 2.\F. (1-+1)12) (n TO( 3 00 1 S I i StJ+vIT~ 0.66 (- 6 617\ (NI R0 3 3 3 333V (N-? ) + 0.4 1o6-6~ iofn 3,o n n=L,LA,? PF TI J ) 28

3.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 ). RJPPHJ N F 0 F C(NA ( P I-) fl I IF I rN L COnmmn /snO/IPS,(20) F)U 5 I=1,N\I JPS(II)=I NM I =NlDO l K=1,NM1 I F=0.0 On i I =K N I P= I PS (I) IF(ABS(U)L(IPqK)).LF, RIF) O TOll1 RIF=4RS(UJL(IPK)) l)XPI V= I 11 C, N TI N(UF IF(IDXPIV.FQ. K) c;o TO 15 J= I PS(K) I P S(K)=IPS( I X,IV IPS( IDXPTV)=J) 15 KP=IPS(K) P 1 \/ 0T= J L (KP,K) KPI=K+1. on 16 I=KPlN IP=IPS( I) FM=-L(t IPK/P IvOT iJL( IP, K)=EM DO 16 J=KP1,N UL ~(I P = )=UL(I P, I ) +FMW(L ( KPp J 16N TCO NTINU RE TURN END 29

553 S-iWpnItT I fNES)LVF (1- i1Xr (,4) f')IMFNSIflN tJL —(80,80),FR(R),X(8O) Nj P1I= N +I I P= I PS( 1) X (1) =k( I P) Dfn 2 I=2,?N I P= I PS (I) SI A 1 = I -1 SUJM=SUM+UL (I PJ ),-,X (J) 2 X(I)=R(IP)-SIUM I P = I P S (N ) X (N)=X (N)/UL ( IPN) Dfn 4 IPiAC-K=?,wN I =NP1- IR4ACK I P= I PS (I) fIf 3 J=IP1,pN 4 X (I )= (X (I)-SUJM) /IJLlkITPp ) RE TURN F Nn 30

553 3.6 SUBROUTINE ELLI (Ml, K, E, KPR, KEY) This computes the elliptic integrals K (m) and E (m) and the derivative K' (in) from their power series approximations (see Abramowitz and Stegun, 1964, p. 590). Arguments: Ml1 Real, the quantity (1 -i). K Real, K(m). E Real, E(m). KPR Real, K'(m). KEY Integer: 0 Compute K, E and KPR; 1 Compute K, E but omit KPR. Y FH Rf1UTITT'f t:F FL- ( MI, KgF K PR % K FY) RFAL MI,V PR T=-A4 LfG ( mlI K =l.4 38 6 29; 4-. 5 T + M (9 56 6 344 E - 2 + 1. 2 4 9 8159T+M( + 3M.3 9 0 9 2 E - 2 5 + 6 8l0 P 4 E' -? ' T4.Ml1-:-( 3,74 564E-2+3.328'355E-2*T Ml.451196E —2 S +4.41797'F-3"-T)) -=l0. 0+Ml: (.443251 14+. 2499 837T+Ml* ( 6.26061E-2+9. 200 01E-2'T+MI ( F 4.757384 F- 2+4.0 69698F-2,T+Mi' (l(.736506E-?2+5.264496E —3-.T)))) JF(KE-Y F). I) RFTURN K P F=.5/M1 + 2.83??5F-2 -.l249859 -T + MI (-2,999362E-3 —.1376O5"i 6 +MI*(-7.899336F-2 - 9.985066F-2:-:T + M1I(-5.362998E-2 - 5 1.767148 F-? T ) RF TURN F ND 31

553 A. 7 SUBROUTINE SETUP (KEYP1l, KEYMll, I, J, L, APIll, API33, AMI11, IJ) This is essential in computing the linear systems. Specifically SETUP, after calling ELLI, computes the quantities API 11 (K of eq. (8), API33 (K0 of eq. 12), n (eq. 20) and ~2 (eq. 21). The quantities Ql and Q2 are used to compute f (i, j) (AMIll) of eq. (39). Arguments: KEYP11 0 when computing API11 and API33, else 1. KEYM 1 0 when computing AMI11, 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 AMI 11 IJ 0: use last value of M 1 in kernel computations; 1: compute new Ml. 32

55 SHRR171tJTJINIF SFTU P(KEYP 1l tKF Y M II, IJ,JL, AP II,?\PI133,9Af III, 1.J) CflMMON RHFI(W p9),7(FO9),APC,(8n),C(8o9),s(80,9) PEA,~L M V, I K PK P 7 D = 7- (3 L ) -Z 7 1 I =PRHO ( I:r ~8 I _IJ ~FO 0) 0 Tn ll5 PPR P = R~ P P Al R Pi P + R R PI A4= - R P P-~ R PI + 7 D 04,7 0 M1= ( 2 —AI)/I( A2+41I M=l1. -MI C A LL F L.L I ( MIwK vEtK PPRK FY MlII A\ 0= M / P.F P A 1 = SQ P, T ( A. 0),A? = 4A +\M A 3= 2.-m I F (KFY PII -I) 113,vI114 pI114 113 API`Ill=A1*(A3`K-E-F) /M A PI 3 3: AlI "K I114 A I=AlI:4A0 A 3=5Aa ll'i AO=C(,JL) n MI=.A1*I.5' K- A 3 *KP R nM2=Ai* (F-A'(3M-, - A24'A3*KPR6))/(M,'M) Am Ill1:R,,,A, 0M 2 + ( 7 D0*5,(3J ~L R P4*AO'4F!M 1 1 1 6 P FTI RN~ F N D 33

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. Keller, J.B., It.E. Kleinman and T.B.A. Senior (1972), Dipole Moments in Rayleigh Scattering, J. Inst. Math. and Applics., 9, 14-22. Senior, T. B. A. and D. J. Ahlgren (1972), The Numerical Solution of Low Frequency Scattering Problems, The University of Michigan Radiation Laboratory Report No. 013630-9-T. 34