RL-706 = RL-706 RAYLEIGH SCATTERING BY DIELECTRIC BODIES* Purpose: To compute the tensor elements specifying the low frequency (Rayleigh) scattering of a plane electromagnetic wave by a homogeneous dielectric body with axial symmetry. Language: Fortran IV Authors: Thomas B.A. Senior and Thomas M. Willis III, Department of Electrical and Computer Engineering, University of Michigan, Ann Arbor, MI 48109 Avail abi 1 i ty: ASIS-NAPS document NAPSDescription: When a plane electromagnetic wave is incident on a homogeneous dielectric body immersed in free space, the leading (Rayleigh) term in the low frequency expansion of the far zone scattered field can be expressed in terms of the electric and magnetic polarizability tensors P and M. These are functions of the relative permittivity cr and relative permeability }r, respectively, and of the geometry of the body, but are independent of the excitation. If the body is symmetric about the x axis of a Cartesian coordinate system 3 x,x,x, the tensors diagonalize and P = P, M = M. As shown in 1 2 3 22 11 22 11 [1], P and M are particular cases of a general polarizability tensor X(T) with P = X(r ), M = -X(p ). The tensor elements are weighted integrals of potential functions which are obtained by solution of the corresponding integral equation by the moment method. The program computes X (T) and X (T), corresponding to the transverse and axial 11 33 components of the dipole moments respectively, for any complex T, * This work was supported by the U.S. Army Armament Research and Development Command, Chemical Systems Laboratory, under Contract DAAK11-81-K-0004.

-2 - and is a modified version of an analogous program [2.3] for perfect conductivity (cr = c, r = 0). Although it is difficult to state the highest frequency at which the Rayleigh scattered field is applicable, a reasonable criterion is that the free space and interior wavelengths should both exceed the maximum linear dimensions of the body by a factor of 10 or more. The body is described by its profile in the z(=x ),p plane of a 3 cylindrical coordinate system. The profile is assumed to be composed of up to 15 straight line or circular arc segments. The former are described by their end points in the z,p plane, while the latter are described by their end points, the angles subtended at their centers of curvature and their concavities. Each segment is uniformly subdivided into cells and the total number N of cells across the entire profile must not exceed 102. For many shapes however, a much smaller value of N is adequate. Concave (re-entrant) shapes can be treated, as well as multiply connected ones, e.g., a torus. Bodies having two or more disjoint portions can also be considered provided their material parameters are the same and they have a common axis of symmetry. An input data set consists of one control card and a number of segment specification cards, one for each segment of the profile. The control card lists: SEGMENTS: number (<15) of segments EXCLUSION: a decimal fraction involved in the computation of a self cell; defaults to 0.001 TAU: complex material parameter PRINT; prints (PRINT = 1) or suppresses (PRINT = 0) the two potentials obtained by solution of the integral equation; defaults to 0

-3 - A segment specification card lists: CELLS; number of sampling TYPE KEY: type of segment: 1 arc, concave down; VOLUME SENSE: + on blank: additi -: subtractive volL cells on the segment circular arc, concave up; 2 circular 3 straight line ve volume; Z-COORD. ENDPOINTS: RHO-COORD. sel f-explanatory ENDPOINTS: self-explanatory THETA(DEG.): angle subtended by a circular arc at its center of curvature As an illustration, consider a sphere of radius 0.5 and T = 2+i with N = 20. The corresponding profile is a semicircle in the z,p plane of radius 0.5 centered (say) at (0.5,0). The print out is as follows:..... 111:ft..f:c::c::...L i...I I I. I II A C11 *II,. "....... ) ). I:.::1. -I I I

-4 - The known (exact) results are X /V = X /V = 0.88235... + 11 33 iO.52941.... In the sample run the imaginary part of X /V was 11 computed with least accuracy (0.2 percent). The c.p.u. time for the Amdahl 470/V8 computer at the University of Michigan was 0.14 s. Increasing N to 30 reduced the errors by about a factor 2, but at the expense of a factor 2 increase in the c.p.u. time. References [1] T.B.A. Senior, "Low frequency scattering by a dielectric body," Radio Sci., Vol. 11, pp. 477-482, 1976. [2] T.B.A. Senior and D. J. Ahlgren, "Rayleigh scattering," IEEE Antennas Propag. (Comp. Prog. Descr.), Vol. 21, p. 134, 1973. [3] T.B.A. Senior and D. J. Ahlgren, "The numerical solution of low frequency scattering problems," The University of Michigan Radiation Laboratory Report No. 031630-9-T, February 1972.

RL 706 COMPUTER PROGRAM: DIELCOM T.B.A. Senior and T. M. Willis III Radiation Laboratory Department of Electrical and Computer Engineering The University of Michigan Ann Arbor, Michigan 48109 1. Introduction The program solves two integral equations and then computes the two independent elements X and X of the general polarizability tensor X for 11 33 a homogeneous dielectric body whose material parameter T,representing either the relative permittivity or relative permeability,may be complex. The body is assumed symmetric about the x (=z) axis of a Cartesian coordinate 3 system x,x,x. The theoretical background and the definition of the I 2 3 the tensor X are given in [1], and the program itself is a modified version of one previously developed [2] for a perfectly conducting body. The numerical procedures will be discussed and the various subroutines presented along with a flow chart of the program. 2. Formulation 2.1 Body Description It is convenient to introduce a cylindrical polar coordinate system p,4,z whose z axis is the axis of symmetry of the body. The profile of the body in the p,z plane is defined as p = p(z), and is constructed of segments which may be linear arcs or circular arcs of either concavity.

-2 - Distance along the profile is measured by the variable s, and the tangent at any point forms the angle a with the z axis: -f < a < -. The volume of rotation of a linear arc (Type 3) segment whose endpoints in the z,p plane are (z,p ), (z,p ) is 1 2 2 z 2 6V = Tp2(z) dz 1 = - (Z - Z )(p2 + p p + p2) 2 2 1 2 1 2 1 (1) For a circular arc which subtends an angle 0 at its center (z,pO) and is of radius a, let ~-1 I if segment is concave down (Type 1) if segment is concave up (Type 2) Then V = ( - ) + a - 1 (u2 + u + u2) 2 IP ~ 1 1 2 2 + p u (p - p )- u (p - p )- ~da2o 0 2 2 0 1 I 0o (2) where Z - ZI d 2 1 d = - 2 1 and u = z - z 2 2 0 u = Z - Z 1 1 0

-3 - 2.2 Integral Equations Given the profile p = p(z) and the volume V, quantities m and m are defined as 1 m = m(z,p,z',p) = 4pp (p - p,)2 + (z - z')2 (3) and m = 1 - m, 1 (4) where (z,p) and(z',p') are distinct points on the profile. The unprimed variables denote the observation point and the primed variables the point of integration. Note that 0 < m,m < 1. If K(m) and E(m) are the complete elliptic integrals of the first and second kinds respectively, let 3/(m + 2mK'( Qo =( I f(m) + 2m K' (m) (5) 3/2 =1 ( (6) - K'(m (1) - Km) (6) 3/2 2 2 = p p(2,) {2m ( K'(m) - (1 - K(m) + E(m (7) where K'(m) = (d/dm)K(m). For a body whose complex material parameter is T, the integral equations to be solved are (see [1,3])

-4 - 1+ Tr 1 - TW (s) = 2p + w (s') lp cos a' + [(z' - z)sin a' - p' cos a'] p' ds' (8) and -1 W (s) = 2z + -rlW (s') cos a'Q + [(z' - z)sin a' 1 -T 3 7j 1 - p' cos ] p' ds' (9) oj where the integration is over the entire profile of the body and the bar across the integral sign implies a Cauchy principal value. In terms of the solutions, 1 V JP'W (s') cos a' ds' (10) V3 - 2 p'W (s') sin a' ds'. (11) where V is the volume of the body. 3. Numerical Solution of the Integral Equations The numerical procedures involved in computing X /V are 33 similar to those required for X /V, and it is therefore sufficient 11 to give the details only for the latter.

-5 - 3.1 X11/V Computation The primary task is the solution of the integral equation (8) for the function W (s) and this entails the determination of a sequence of values i) i = 1,2,...,N, approximating Wl(s) at the sampling points s = Si 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 (zi,p.). Within each cell we also define the points s. and s. where 1- 1+ si = si - 0A5. (12) S = + asj i+ i 0 1 with the restriction 0 < a < 1/2. O - By assuming that W (s) has the constant value W(j) over the jth cell, the integral on the right-hand side of (8) can be evaluated as a linear combination of the W) whose coefficients depend on the 1 position (Z,p) of the field point, leading to a linear system of N equations in N unknowns, viz. N - (j p. cos a' 2 + [(z' - z.)sin a' - p' cos a'] p' ds' j= C. 2rp, - T (13) 2api- 1 - T 1 for i = 1,2,...,N. Hence, the system to be solved is

-6 - Aw = b 1 where w is a column vector with elements 1 j = 1 A is a square matrix with elements (14) a.. a oi =- f f(s ')p' ds' C. J (i f j) (15) a. = - f f(s Cj where f(s') = pi cos a' 2 1 2 1+ T ')p' ds' + I T (i j) - -T + J(z' - zi)sin a' - p' cos a']Q 1 (16) and b is a row vector with elements bj = 2TPj j = 1,2,...,N The integrations were performed using three point Gaussian quadrature for which the parameter ao in (12) is o = 1/2 /3/5. For i f j, we then have a j [ 1 i pj f(i,j-)+ pj+f(i+ +f(i,j) As (17) for i,j = 1,2,...,N, in which

-7 - f(ij) = pi cos aQ (i,j) + (zj - z )sin a. - p. cos a. Q (i j) (18) where the subscripts j- and j+ correspond to the points s. and s. of (12). The diagonal elements a of the matrix A are (12). The diagonal elements aii of the matrix A are si +(1/2)AS i 1 +T _ aii = - T si- (1/2)Asi p' f(s' )ds' and require the numerical determination of a Cauchy principal value. As an approximation to this principal value, we remove from the cell C. a slice defined by the interval (s. - (1/2)BAs., s. + (1/2)3As.) where A, 0 < g < 1, is the fractional exclusion; B = 0 implies no exclusion, i.e., that the principal value is not taken. We now have a Tr1 + T ai i = 1 ' si —(V2)Asi si- (1/2)Asi p'f(s')ds' - s.+(l/2)As. p "f(s') ds' (19) and these integrals are also computed using three-point Gaussian quadrature. Defining

-8 - i2 S. 11 S. 13 i5 sj6 1 = s - - (1 + B)As S2 - (1 - B)as. Si2 + ao(l - 3)As 1 Si= S+ + si + = S - - o(1 - ) aS 1= + ( - ): si5 + %o(1 - )ASi I (20) -I we obtain a..(1 - - T B) + I 1 - T 2 + Pi6f(ii) + 4Pi - pi, 18 P f(ii ) + P 3f(i i3 ) + Pif(ii ) -8L fii 1 i3 '3 i4 '14 f(i i ) + p 5f(i i) As i i)] (21) Equations (17) through (21) completely describe a system of N linear equations in N unknowns W(i), i = 1,2,...,N. Having determined the sampled values Wi) = W (si), X /V is computed from (10) by 11 11 integration over each segment of the profile using a second order integration procedure (subroutine INTEG).

-9 - 3.2 X33/V Computation The procedure for solving (8) and then computing X /V from (11) is 33 about identical to that described above. The matrix equation for the sampled values W (s) W(i) is Aw b where:3 1 3 3 w3 WJ ) j 1,2,...,N, (22) b. = 2rz. 3 3 and for i,j = 1,2,...,N with i f j, [ f 4 ] aij - ~ Pj g(i,j-) + pj+g(i j+) + g(ij) ASi (23) with g(i,j) = P cos a (i,j) + (z - z )sin a. - pi cos a. (i j) 1 j i j jo '3 (24) The same approximation as before is used to compute the diagonal elements a i of the matrix A, and thus a.. -- lp 1 (1 i5 p g(Lii) + Pi gi + pg(ii ) 1 2 18 1 3 i 4 I 3 4 + Pi g(i,i )} + g(i,i ) + g(i,i ) AS(25) 6 2 5 for i = 1,2,...,N.

-10 - Equations (22) through (25) completely describe a system of N linear equations in N unknowns W( ) j = 1,2,...,N. Once the sampled 3 values Wj) = W (Sj) are found, X /V is computed from (11) by 3 3 33 integration over each segment of the profile using a second order integration routine as before. 4. Computer Program The program computes X /V and X /V, and consists of a main program and six subroutines. 4.1 Data Set A data set is made up of a control card and a number of segment specification cards, one for each segment (or sub-segment) of the profile. Control Card Columns 1-5 6-17 18-41 42-46 Description Integer number of segments in the body (must be <15). The fraction exclusion B. If these columns are blank, e defaults to 0.001. Complex material parameter T. Integer print key; 1: print W and W potential functions 1 3 0: do not print W and W (default is 0). 1 3

-11 - Segment Specification Card Col umns 1-5 6-10 11-13 14-25, 26-37 38-49, 50-61 62-73 Description Integer (right justified): the number of sampling points or cells on the segment. Segment type key: (right justified) 1: circular arc, concave down 2: circular arc, concave up 3: linear. Volume sense: (single character, right justified) + or blank: additive volume -: subtractive volume. Two real numbers: respectively, the end coordinates z and z of the segment. 1 2 Two real numbers: respectively, the end coordinates p and p of the segment. 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 102. The profile is specified in the direction of increasing profile-length, beginning either at its left-hand intersection with the z-axis and ending at its right-hand intersection with the z-axis, or beginning and ending at its own left-most point. Re-entrant segments are permitted, allowing

-12 -z > z, and the "body" may consist of two or more disjoint parts. 1 2 For maximum accuracy in the X33 computation the body should be approximately centered about the origin.

-13 - 4.2 Main Program The main program reads, prints and supervises all computations. A rough flow chart showing the interaction of the subroutines is given below. START Read control card L For each segment (1) read specification card (2) compute sample END OF -- STOP DATA DATA -^I Construct linear systems DECOMP SOLVE - Print Results - Solve linear systems \ I — N --- SETUP ELLI INTEG I Integrate weighted sequences

-14 - 4.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. Comments: Stored in COMMON are the arrays RHO(102,9), Z(102,9), ARC(102), C(102,9) and S(102,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 si, si and si+ respectively of (12). 4.4 Subroutine INTEG (V, NSEG, NUMPTS, SUM) INTEG numerically integrates quadratic interpolating polynomials approximating the complex valued 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. Comments: Stored in COMMON are the lengths ARC (I), I = 1,...,N required to compute the integral.

-15 - 4.5 Subroutines CDLUD (N,ADIM,A,IV) and CDBS (N,TDIM,T,IV,B,X) Used together, CDLUD and CDBS solve the linear system AX = B. CDLUD performs a L-U decomposition of the N x N matrix A and CDBS performs back-substitution. These routines are adapted from Forsythe and Moler [3]. 4.6 Subroutine ELLI (Ml, K, E, KPR) This computes the elliptic integrals K(m) and E(m) and the deriviative K'(m) from their power series approximations [4]. 4.7 Subroutine SETUP (I, J, L, AXIll, AXI33) This is essential in computing the linear systems. Specifically SETUP, after calling ELLI, computes the quantities S (Eq. 5), s (Eq. 6), and Q (Eq. 7). These are then used to compute f(i,j) (AXIII) of 2 Eq. (18) and g(i,j) (AXI33) of Eq. (24). A listing of the entire program follows:

-16 - C C C DIELCOM DEC. 1981 C C C C********************************************** ********************C C******** *********************** * ****************** ******************C C C C The PROGRAM DIELCOM computes the elements of the C C polarizability tensor for a rotationally symmetric C C homogeneous dielectric body. See "Radio Science", vol. 11 C C No. 5, pages 477-482, May 1976: Low-frequency scattering C C by a dielectric body, by T.B.A. Senior C C C C For a discussion of a similar program (MINK) see The C C University of Michigan's Radiation Laboratory report, C C "The Numerical Solution of Low Frequency Scattering C C Problems", by Thomas B. A. Senior and David J. Ahlgren; C C report no. 013630-9-T. C C C C(************************ ****************************************C C C C LOGICAL DEVICE ASSIGNMENT C C C C 5: FORMATTED INPUT PARAMETERS (SEE INPUT FORMAT) C C 6: FORMATTED OUTPUT OF TENSOR ELEMENTS C C C c ***************************************************************** *C c C C INPUT FORMAT C C C C CARD#1 FORMAT(I5,3E12.6,I5) C C NSEGS NUMBER OF SEGMENTS IN THE BODY C C (NSEGS MUST BE <=15) C C FR FRACTIONAL EXCLUSION (BETA) USED WITH CELLS C C TAU COMPLEX MATERIAL PARAMETER C C PRINTW INTEGER PRINT KEY; PRINTW=1 PRINTS THE W C C POTENTIALS, PRINTW=0 DOES NOT (DEFAULT IS 0) C C C C CARD#2 TO CARD#1+NS1+NS2 FORMAT(2I5,2X,A1,5E12.6) C C NUMPTS NUMBER OF CELLS THIS SEGMENT IS TO BE DIVIDED C C UP INTO (THE SUM OF ALL NUMPTS MUST BE <=102) C C ITYP TYPE OF SEGMENT; ITYP=1 CONCAVE DOWN SEGMENT, C C ITYP=2 CONCAVE UP SEGMENT, ITYP=3 LINEAR C C ISIGN CHARACTER CORRESPONDING TO THE VOLUME SENSE OF C C THIS SEGMENT; ISIGN='+' POSITIVE VOLUME, C C ISIGN='-' NEGATIVE VOLUME (DEFAULT IS '+') C C ZEP(1) Z CO-ORDINATE OF THE 1ST END POINT OF SEGMENT C C ZEP(2) Z CO-ORDINATE OF THE 2ND END POINT OF SEGMENT C C RHOEP(1) RHO CO-ORDINATE OF THE 1ST END POINT C C RHOEP(2) RHO CO-ORDINATE OF THE 2ND END POINT C C THETA ANGLE (IN DEGREES) SUBTENDED BY SEGMENT C C C C**********************!*******************************************C COMPLEX AX11(102,102),AX33(102,102),X(102),TAU,POM,X11,X33 DIMENSION B(102),ZEP(2),RHOEP(2),ST1(6),ST3(6) INTEGER NUMPTS(15),PLUS/'+'/,BL/' '/,INDX(2),PRINTW COMMON RHO(102,9),Z(102,9),ARC(102),C(102,9),S(102,9) COMMON /SOL/ IPS(102)

-17 - DATA MIN,TWOPI,PI,WO,W1/'-',6.283185,3.141593,.4444444,.2777777/ C CCC READ, CHECK, AND WRITE INPUT VALUES C 37 READ(5,34,END=999) NSEGS,FR,TAU,PRINTW WRITE(6,4) NSEGS IF(NSEGS.LE. 0.OR. NSEGS.GT. 15) GO TO 990 IF( FR.LE. 0..OR. FR.GT. 1.) FR=.001 WRITE(6,5)FR WRITE(6,150) TAU M=0 V0=0.0 IC CCC READ AND WRITE INPUT PARAMETERS FOR EACH SEGMENT C DO 11 I=1,NSEGS READ(5,12) NUMPTS(I),ITYP,ISIGN,ZEP,RHOEP,THETA IF(ISIGN.EQ. BL) ISIGN=PLUS WRITE(6,13) I,NUMPTS(I),ITYP,ISIGN,ZEP(1),RHOEP(1),ZEP(2), & RHOEP(2) IF(NUMPTS(I).LT.0.OR. ITYP.LE.0.OR. ITYP.GT.3) GO TO 990 IF(ITYP.NE. 3) WRITE(6,14) THETA IF (NUMPTS(I).EQ.0) GOTO 11 N=M THETA=PI*THETA/180.0 M=M+NUMPTS(I) IF(M.GT. 102) GO TO 990 CALL DATA(ITYP,N,M,ZEP,RHOEP,THETA,FR,VOLINC) IF(ISIGN.EQ. MIN) VOLINC=-VOLINC VO=VO+VOLINC 11 CONTINUE C CCC SET UP LINEAR SYSTEMS OF EQUATIONS TO BE SOLVED C DO 2 N=1,M INDX(1)=N AN=ARC(N) TN=RHO(N,8) C DO 3 L=N,M IF(L.EQ. N) GO TO 82 AL=ARC(L) TL=RHO(L,8) AX11(N,L)=CMPLX(0.,0.) AX11(L,N)=CMPLX(0.,0.) AX33(N,L)=CMPLX(0.,0.) AX33(L,N)=CMPLX(0.,0.) INDX(2)=L C DO 103 J=1,3,2 JP6=J+6 C DO 104 LL=1,2 I=3-LL I1=INDX(LL) I2=INDX(I) TI2=RHO(I2,JP6) CALL SETUP(I1,I2,JP6,AXI11,AXI33) AX11(I1,I2)AXX11(I1,I2)-AXI11*TI2 AX33(I1,I2)=AX33(I1,I2)-AXI33*TI2

-18 -104 CONTINUE IC 103 CONTINUE C CALL SETUP(NL,8,AXI11,AXI33) AX11(NL)=AL*(W1*AX11(N,L)-AXI11*WO*TL) AX33(NL)=AL*(Wi*AX33 (NL)-AXI33*WO*TL) CALL SETUP(LN,8,AXI11,AXI33) AX11(LN)=AN*(W1*AX11(LN)-AXI11*WO*TN) AX33(L,N)=AN*(Wi*AX33(L,N)-AXI33*WO*TN) GO TO 3 82 IF(FR.EQ. 1.0) GO TO 3 C DO 86 I=1,6 86 CALL SETUP(N,N,I,ST1(I),ST3(I)) C U=.5*(1.0-FR)*AN POM=PI*(CMPLX(1.,0.)+TAU)/(CMPLX(1.,0.)-TAU) AX11(N,N))=POM-U* (WO* (RHO(N2)*ST1(2)+RHO(N5)*ST1(5) & +W1*(RHO(N,1)*ST1 (1 )+RHO(N,3)*ST1(3)+RbO(N,4)*ST1(4)+ & RHO(N,6)*ST1(6))) AX33(N,N)=POM-U*(WO*(RHO(N,2)*ST3(2)+RHO(N,5)*ST3(5)) & +W1*(RHo(N 1)*ST3(1 )+RHO(N,3)*ST3(3)+RHo(N,4)*ST3(4)+ & RHO(N,6)*ST3(6))) 3 CONTINUE 2 CONTINUE CCC SOLVE MATRIX EQUATIONS AND INTEGRALS TO GET Xli AND X33 C C(CC COMPUTE Xli DO 20 I=1,M 20 B(I)=-TWOPI*RHO(I,8) CCALL DECOMP(AX11,M) CALL SOLVE(AX11,X,B,M) IF (PRINTW.NE.1) GOTO 30 WRITE(6,44) (I,Z(I,8),RHO(I,8),X(I), I=1,M) 30 CONTINUE C DO 26 I=1,M 26 X(I)=RHO(I,8)*C(I,8)*X(I) IC CALL INTEG(X,NSEGS,NUMPTS,xll) Xl 1 =PI /VO*~X1 1 IC CCC COMPUTE X33 C DO 24 I=1,M 24 B(I)=-TWOPI*z(I,8) C CALL DECOMP(AX33,M) CALL SOLVE(AX33,X,BM) C IF (PRINTW.NE.1) GOTO 35 WRITE(6,46) (I,Z(I,8),RHO(I,8),X(I), 1=1,M) 35 CONTINUE

-19 - DO 22 I=1,M 22 x(I )=-2.*RHO(I,8)*X(I )*S(I,8) IC CALL INTEG(X,NSEGS,NUMPTS,X33) X33=PI/VO*X33 IC CCC WRITE OUTPUT IC WRITE(6,52)V0 WRITE(6,200) Xli WRITE(6,210) X33 GO TO 37 990 WRITE(6,991) 999 STOP C (.CcC FORMATS C 34 FORMAT(15,3E12.6,I5) 4 FORMAT('1***** PROGRAM DIELCOM *** & /'0',5X,'SEGMENTS',4X,'=',I2) 1 FORMAT(' ',5X,'BODY #2',5X,'=',12) 5 FORMAT(' ',,5X,''EXCLUSION =',F7.4) 150 FORMAT(' ',5X,'TAU ='E12.6j,' +J',El2.6) 12 FORMAT(2I5,2X,,A1,5E12.6) 13 FORMAT('OSEGMENT #'r12,';'/ ' ',,5X,,'CELLS',9X,.'= '1I2/' '.,5X, & 'TYPE KEY',6Xr rv'='I2/' ',5X,'VOLUME SENSE =,?A4/' ',5X#, & '1ST END POINT = ('F1O.5,', ',rFlO.5,')'/' I',5X, & '2ND END POINT = (', F1O.5,', ',r1F1O.5,'v)') 14 FORMAT(' ',5X,,'THETA (DEG) ='0,F10.5) 52 FORMAT('OCOMPUTED RESULTS;'!' ',5X,'VOLUME = 'F10.5) 200 FORMAT(' ',5X,'Xil/V = ',F1O.5,' +J',F1O.5) 210 FORMAT(' ',5X,'vX33/V = ',FiO.5,'I +J',FiO.5) 44 FORMAT( '0',6X, 'I',lix,,'z', 15X, 'RHO',21X, 'WIl(I) '/ & (' ',5X,I3,5X,2(F12.6,5X),E12.6,' +J',E12.6)) 46 FORMAT( '0',6X, 'I''lix,' z', 15X, 'RHO',21X, 'W3(I)'I/ & (' ',5X,I3,5X,2(F12.6,5X),E12.6,' +J',E12.6)) 991 FORMAT('0*** ERROR IN DATA') END

-20 - C C C C SUBROUTINE DATA(IN, MX, MY, ZEP,RHOEP,THETA,B,VOL) iC C C C C SUBROUTINE DATA is called by the MAIN PROGRAM for each C C segment profile. It computes the z and rho co-ordinates, C C sin(alpha), cos(alpha), and the arc length for each cell. C C Arguments: C C IN Type key for segment C C MX Total number of cells in segments to left C C MY MX + (number of cells in this segment) C C ZEP z - coordinate end points of segment: C C ZEP(1)=zl, ZEP(2)=z2 C C RHOEP Rho coordinate end points of segment: C C RHOEP(1)=rho1, RHOEP(2)=rho2 C C THETA Angle (in radians) subtended by circular C C arc at its center C C B Fractional exclusion (beta) C C VOL Incremental volume of segment C C C DIMENSION ZEP(2),RHOEP(2) COMMON RHO(102,9),Z(102,9),ARC(102),C(102,9),S(102,9) DATA STEP/.3872988/ C MXP1=MX+ 1 EN=FLOAT(MY-MX) IF(B.NE. 1) SUBSTP=.5*(1.0-B)*STEP IF(IN-2)1,2,3 C CCC TYPE KEY 'IN' = 1 (CONCAVE DOWN SEGMENT) C 1 CC=-1.0 GO TO 10 C CCC TYPE KEY 'IN' = 2 (CONCAVE UP SEGMENT) C 2 CC=1.0 C CCC COMPUTE THE RADIUS OF CURVATURE C 10 ST2=SIN(THETA/2.0) A=ZEP(2)-ZEP( 1) RAD=0.5*SQRT((RHOEP(1)RHHOEP(2))**2+A*A)/ST2 C CCC COMPUTE THE CENTER OF CURVATURE CO-ORDINATES C DD=A/ABS(A) T=CC*DD*COS (THETA/2.)/ST2 ZCNT=0.5*(ZEP)EP+ZEP(2)+T*(RHOEP( 1 )-RHOEP(2))) RHOCNT=0.5*(RHOEP( 1 )+RHOEP(2)+T*A) C CCC COMPUTE THE INCREMENTAL VOLUME C U2=ZEP(2)-ZCNT U1=ZEP( 1 )-ZCNT VOL=3.141593*ABS(A*(RHOCNT**2+RAD*RAD-(U2**2+U1*U2+U1**2)/3.0)

-21 - & -CC*RHOCNT*(U2*(RHOEP(2)-RHOCNT)-U1*(RHOEP( 1 )-RHOCNT) +RAD*RAD & *DD*THETA)) C CCC COMPUTE COS(ALPHA), SIN(ALPHA), Z, AND RHO AT THE BOUNDARIES CCC OF EACH CELL, AND THE ARC LENGTH OF EACH CELL IN THE SEGMENT IC BETA=CC*DD*THETA/EN THET1=ATAN2(RHOEP(1)-RHOCNT,ZEP(1)-ZCNT) U=ABS(BETA*RAD) B3=STEP*BETA C DO 902 I=MXP1,MY PHI=THET1+(I-MX-.5)*BETA IF(B.EQ. 1.0) GO TO 1905 DO 1902 J=1,2 ANG=PHI+.5*(J-1.5)*BETA*(1.0+B) C DO 1903 L=1,3 PSI=ANG+(L-2)*SUBSTP*BETA M=L+3*(J-1) C(I,M)=-CC*SIN(PSI) S(I,M)=CC*COS(PSI) Z(I,M)=ZCNT+RAD*CC*S(I,M) 1903 RHO(I,M)=RHOCNT-CC*RAD*C(I,M) C 1902 CONTINUE C 1905 DO 903 J=7,9 ANG=PHI+(J-8)*B3 C(I,J)=-CC*SIN(ANG) S(I,J)=CC*COS(ANG) Z(I,J)=ZCNT+RAD*CC*S(I,J) 903 RHO(I,J)=RHOCNT-CC*RAD*C(I,J) C 902 ARC(I)=U C RETURN C CCC TYPE KEY 'IN' = 3 (LINEAR SEGMENT) C CCC COMPUTE COS(ALPHA), SIN(ALPHA), Z, AND RHO AT THE BOUNDARIES CCC OF EACH CELL, AND THE ARC LENGTH OF EACH CELL IN THE SEGMENT C 3 DX=(ZEP(2)-ZEP(1))/EN DY=(RHOEP(2)-RHOEP( 1) )/EN U =SQRT(DX*DX+DY*DY) SI=DY/U CI=DX/U C DO 917 I=MXP1,MY PHI=FLOAT(I-MX)-.5 IF( B.EQ. 1.0) GO TO 1800 C DO 1802 J=1,2 ANG=PHI+.5*(J-1.5)*(1.0+B) C DO 1803 L=1,3 M=L+3*(J-1) PSI=ANG+(L-2)*SUBSTP

-22 - Z (I,M)=ZEP( 1 )+PSI*DX RHO(I,M)=RHOEP( 1)+PSI*DY S(I,m)=Si 1803 C(IM)=CI C 1802 CONTINUE C 1800 DO 913 J=7,9 ANG=PHI +(J-8) *STEP Z(I,J)=ZEP( 1)+ANG*DX RHO(I,J)=RHOEP( 1)+ANG*DY CUij)=CI 913 S(IJ)=SI IC 917 ARC(I)=U IC.CCC COMPUTE THE INCREMENTAL VOLUME VOL=1.047198*(ZEP(2)-ZEP(1) )*(RHOEP( 1)**2+RHOEP( 1)*RHOEP(2)+ & RHOEP(2)**2) RETURN END

-23 - C C C C SUBROUTINE INTEG(V,NSEG,NUMPTS,SUM) C C C C C SUBROUTINE INTEG is called by the MAIN PROGRAM to C C numerically integrate quadratic interpolating polynomials C C approximating the data on each segment of the profile. C C Arguments: C C V Vector of function values, ordered as cells C C NSEG Total number of segments in the profile C C NUMPTS Vector containing in NUMPTS(I) the number of C C cells on the Ith segment: I = 1, NSEG C C SUM Integral of V across the profile C C C COMMON RHO(102,9),Z(102,9),ARC(102),C(102,9),S(102,9) COMPLEX V(102),SUM DIMENSION NUMPTS(15) SUM=0.0 JACC=1 DO 3000 I=1,NSEG T=ARC(JACC) L=NUMPTS(I) N=L+JACC-1 SUM=SUM+T*(0.625*(V(JACC)+V(N))-.125*(V(JACC+1)+V(N-1))) IF(L/2.NE. (L+1)/2) GO TO 3001 SUM=SUM+T*(0.6666667*V(N-1)-0.08333333*V(N-2)+0.4166667*V(N)) 3001 LM1=N-1 JLO=JACC+ 1 DO 3002 J=JLO,LM1,2 3002 SUM=SUM+0.3333333*T*(V(J-1)+4.0*V(J)+V(J+1)) C 3000 JACC=JACC+L C RETURN END

-24 - C C C******************************************************************C C C SUBROUTINE DECOMP(UL,N) IC C c ******************************************************** *********C C C C SUBROUTINE DECOMP is called by the MAIN PROGRAM to perform C C a L - U decomposition on the N x N matrix UL. C 4C C COMPLEX UL(102,102),PIVOT,EM,CDUM COMMON /SOL/IPS(102) CABS1(CDUM)=REAL(CDUM)+AIMAG(CDUM) C DO 5 I=1,N 5 IPS(I)=I C NM1=N-1 DO 16 K=1,NM1 BIF=0.0 DO 11 I=K,N IP=IPS(I) IF(CABS1(UL(IP,K)).LE. BIF) GO TO 11 BIF=CABS1(UL(IP,K)) IDXPIV=I 11 CONTINUE C IF(IDXPIV.EQ. K) GO TO 15 J=IPS(K) IPS(K)=IPS(IDXPIV) IPS (IDXPIV)=J 15 KP=IPS(K) PIVOT=UL(KP,K) KP1=K+1 C DO 16 I=KP1,N IP=IPS(I) EM=-UL(IP,K)/PIVOT UL(IP,K)=-EM C DO 16 J=KP1,N UL(IP,J)=UL(IP,J)+EM*UL(KP,J) 16 CONTINUE C RETURN END

-25 - C C C C SUBROUTINE SOLVE(UL, X, B, N) C C c C C SUBROUTINE SOLVE is called by the MAIN PROGRAM to solve C C the linear system UL * X = B by back substitution after C C the UL matrix has been decomposed by the SUBROUTINE DECOMP. C C C COMPLEX UL(102,102),X(102),SUM DIMENSION B(102) COMMON /SOL/I PS( 102) NP1=N+ 1 IP=IPS( 1 ) X(1)=B(IP) DO 2 I=2,N IP=IPS(I) IM1=I-1 SUM=CMPLX(0.,0.) C DO 1 J=1,IM1 1 SUM=SUM+UL(IP,J)*X(J) C 2 X(I)=B(IP)-SUM C IP=IPS(N) X(N)=X(N)/UL(IP,N) C DO 4 IBACK=2,N I=NP1-IBACK IP=IPS(I) IP1=I+1 SUM=CMPLX(0.,0.) C DO 3 J=IP1,N 3 SUM=SUM+UL(IP,J)*X(J) C 4 X(I )=(X(I)-SUM)/UL(IP,I) C RETURN END

-26 - C( C C C SUBROUTINE ELLI(M1,K,E,KPR) C( C C( C C SUBROUTINE ELLI is called by the SUBROUTINE SETUP to C C compute the elliptic integrals K(m), E(m), and derivative C C K'(m) from their power series approximations. (See C C sec. 4.2 of "The Numerical Solution of Low Frequency C C Scattering Problems".) C C C REAL M1,K,KPR T=-ALOG(M1) K=1.386294+.5*T+M1*(9.666344E-2+.1249859*T+M1*(3.590092E-2 & +6.880249E-2*T+M1*(3.742564E-2+3.328355E-2*T+M1*(1.451196E-2 & +4.41787E-3*T)))) E=1.0+M1*(.4432514+.2499837*T+M1*(6.260601E-2+9.20018E-2*T+M1*( & 4.757348E-2+4.069698E-2*T+M1*(1.736506E-2+5.264496E-3*T)))) KPR=.5/M1 + 2.83225E-2 -.1249859*T + M1*(-2.999362E-3-.137605*T & +M1*(-7.899336E-2 - 9.985066E-2*T + M1*(-5.362998E-2 - & 1.767148E-2 * T ))) RETURN END

-27 - C C C C SUBROUTINE SETUP(I,J,L,AXI11,AXI33) C C c*********************** ***************************************** *C C C C SUBROUTINE SETUP is called by the MAIN PROGRAM to compute C C the quantities AXI11 (f(i,j) of eq. 227 in "The Numerical C C Solution of Low Frequency Scattering Problems") and AXI33 C C (a similar quantity for x33). C C Arguments: C C I Subscript of observer (unprimed) cell C C J Subscript of remote (primed) cell C C L Index of the point within remote cell for C C which the kernels are to be computed C C (See SUBROUTINE DATA) C C AXI 11 See above C C AXI33 See above C C C COMMON RHO(102,9),Z(102,9),ARC(102),C(102,9),S(102,9) REAL M,M1,K,KPR REAL*8 DA1,DA2,DM1,DZD,DRP,DR,DRRP DZD=DBLE(Z(J,L) )-DBLE(Z(I, 8)) ZD=SNGL(DZD) R=RHO(I,8) DR=DBLE(R) RP=RHO(J,L) DRP=DBLE(RP) DRRP=DR*DRP RRP=SNGL(DRRP) DA 1=DRRP+DRRP DA2 =DR*DR+DRP*DRP+DZD*DZD DM1=(DA2-DA1)/(DA2+DA1) M1=SNGL(DM1) IF (ABS(M1).LT.1.E-50) WRITE(6,100) M1,Z(J,L),Z(I,8),R,RP,DA1,DA2 100 FORMAT(7E10.3) M=1. -M1 CALL ELLI(M1,K,E,KPR) AO=M/RRP A1=SQRT(AO) A2=M+M A3=2.-M A1=A1*AO A3=.5*A3 AO=C(J,L) OMO=.25*A1 *(K+2.*M*KPR) OM1=-A1*(.25*K-A3*KPR) OM2=A1*(E-A3*((A3+M)*K - A2*A3*KPR))/(M*M) AXI11=R*AO*OM2+(ZD*S(J,L)-RP*AO)*OM1 AXI33=R*AO*OM1+(ZD*S(J,L)-RP*AO)*OMO RETURN END C CCC MODIFICATIONS OF MINK TO DIELCOM BY TOM WILLIS C

-28 - References [1] T.B.A. Senior, (1976), Low frequency scattering by a dielectric body, Radio Sci., Vol. 11 (5), pp. 447-482. [2] T.B.A. Senior and D. J. Ahlgren (1972), The Numerical Solution of Low Frequency Scattering Problems, The University of Michigan Radiation Laboratory Report No. 013630-9-T. [3] G. Forsythe and C. B. Moler, (1967), Computer Solution of Linear Algebraic Systems, Prentice Hall, Englewood Cliffs; pp. 68-69. [4] M. Abramowitz and I. A. Stegun (eds.), (1964), Handbook of Mathematical Functions, NBS Appl. Math. Series No. 55; p. 590.