I,, z - - I - -t,. " TECHNICAL REPOR I*.... 4 -~ I _ I I I _ _ I-.: Dielectric-Covered, Cavity-Backed Aper tures, -. Pisti B. Katehi Department of Electrical Engineering anc Computer Science RL-,842 *': *:.... ' 5 ';,)',,; ' ~ f.:' r: r ', p '.... 4~f -G,t,. ~, - ~ - " '*;;? -'j ' "-. 11~ -^-'^^c^. ~.,,r.. -?'r..,,.;,1:~:; 5.~:. * 4 *~;- ^ *-t.-. >o ti.1 —*<- $iX"*-,...,....., % s-,,,, 4 ~,,-,.,.,, -,..,.,, -, -s r", -D, 1 e e; sU; 4v~~ j 6 t --: '- ';- *-; K v * r, *- r s, < t^ *; -,? * - ' ~-^ s s -- 3k i "'~~~ -. - * t "- 4./ >- ^^. -. ~;l,: i -, b ~' _~ t t /1 * * ~~-~~is: ~~-,S \,;~~~~~,", -.sr, r ~+ 7 - *:^~ ''I {'. ( t -;,..,, r,,... J * * { The University of Michigan Ann Arbor, Ml Lli' rr~a:1. rr~:~ `~ r cr. i. ~'a''64-3i '~ 7' J, `F: t~ --?~; -r RL-842 = RL-842 June 1987

DIELECTRIC-COVERED, CAVITY-BACKED APERTURES 1. General formulation of the problem An important special case of a cavity-backed aperture is that where the cavity is a finite rectangular waveguide. The aperture exists at one end of the the cylindrical waveguide and is covered by a dielectric slab. The other end of the waveguide is completely covered by a conductor. The cavity can then be viewed as a short-circuited waveguide and waveguide theory can be applied. Figure 1 represents a typical problem, where the excitation is by an incident plane wave from the half-space region. The equivalence principle is used to divide the problem into two equivalent problems as shown in figure 2. In region a, the field is produced i -i by the incident field E, H and the equivalent magnetic current -- A -* M =n xE (1) over the aperture region, with the aperture covered by an electric conductor. In region c, the field is produced by the equivalent.-+ magnetic current -M over the aperture region, with the aperture covered by an electric conductor. The fact that the equivalent current in region c is the negative of that in region b ensures that the tangential component of the electric field is continuous across the aperture. The remaining boundary condition to be applied is continuity of the tangential component of the magnetic field across the aperture. -.b The tangential component Htof the magnetic field in region b over the aperture is the sum of that due to the impressed sources

%/M I.J. -P 0 0 I1 -0 0 0 \ v oo o g:\@<<<\\ > \ I-. to a ur HI - - ri I 0 0 a a C- 0 tq 0 x 1-, l< 00. N

x ^*> A n region a (a) Equivalent problem b I H I I region c (b) Equivalent problem c Figure 2

-2 - -*i - -Ob Ht plus that due to the equivalent source M, denoted by Ht (M). Therefore -b - i -b -~ Ht = Ht + H (M) (2) i -b - Both H and Ht(M) are computed with a conductor covering the -> aperture. In region c the equivalent source -M is the only source. Hence, the tangential component of magnetic field in region c over the aperture is -4C -4C - c Ct t Ht = Ht(-M) = -Ht(M) (3) -ob -~c Since Ht should be equal to Ht, equations (2) and (3) give -b - -c -- — i ( 4) Ht(M) + H (M) =-H or ^ f-b - -c -~ A ^ i ('5) n x H (M) + H (M) = -n H Equation (5) is the basic operator equation for determining the equivalent magnetic current M. -ob -Oc The magnetic fields H and H are given in terms of the magnetic current M as shown in the following two integral equations. H.= | - k I + V V F(r/r).M(r')ds' (6) slot cOE)Ld and H = i f. k: I + V V G(r/r') M(r')ds' (7) slot c( c

-3 - In equation (6), Ed is the dielectric constant of the dielectric = —e — e substrate, (Ld is the magnetic permeability and F(r/r') is the dyadic Green's function for the equivalent problem of Fig. 2a. Similarly, e is the dielectric constant of the material in the cavity, Kc is = — -4 -the magnetic permeability and G(r/r') is the dyadic Green's function for the equivalent problem b. In view of (6) and (7), equation (5) gives the following coupled integral equations fxx (r/r') M (r') + Kxy (r/r')M (r')ds' = -H (z=-h) (13) xy y j slot and [K /)M ' + K(r/r' )M )]ds' = -H (z=-h) (14) slot where Kir = j ) (F + G) + VV + ) 3 kd k 2c and i,j = x,y. For a thin slot which is open at the center of the cavity wall as shown in Figure 3, the transverse component of the electric field is much smaller than the longitudinal E<<Ey. This results in one directional magnetic current and therefore a single integral equation: Kx r/r' Mx r' ds' = -H (z=-h) (16) slot

- 4 - In this case, K r/r ) takes a simpler form as shown below K = (F + G) + XX G + ~ F k c xx kd ( 17) with [1 ] F - xx 2 22tkd ( XIr-rf I) 0 ucosh(uz) - Ed cu~sinh (uz) EdU cosh(uh) + usinh(uh) U (18) F - zX 2 2n k co4 1 (xr7 i) oErd U cosh(uh) + uo si~nh(uh) sinh [u(h+z)] ( u cosh(uh) + u0 sinh(uh) and (Appendix A) (19) G r/r') = 4cOe ab n41l _= E E 1 n m k sin kc) - sin ~ Lao nic sin XI + X a 0 1 - ji (2 0) - 0 Cos rra y I - y b 0 Cos [k (C + z + h)].

-5 - 2. Method of Moments To obtain the unknown magnetic current Mx(r' in equation (16), we apply the method of moments. The unknown current is assumed to be in the form: Mxr') = V f (x)g(y') (211) n where w/2 2 g 1 (22) g(y') = W. J -2 -w/2 / 2y! and {f (x'), n=1, 2..., N} is a set of expansion functions overlapping in nature. These expansion functions are given by sink xl-x) I- X < xI < X sin (kl) f(x') = sin[k(x x I)] (23) I — -X < x < X sin(kl)n n+l 0 elsewhere In view of equation (21), equation (16) takes the form L w/2 Vn dx dy' y K(r/r') fn(x')g(yl) = -H + AHR (24) n o -w/2 In order to minimize the error AHX, we impose the following condition

-6 - <AH, > = 0 m = 1,2,...,N (25) or Xm+l | m(x) AHi dx = O at y=O (26) X Xm-1 The testing functions m (x), in equation (26), have been chosen equal to the expansion functions (Galerkin's method). The above condition imposed on equation (24) reduce the integral equation into a set of linear equations which can determine the coefficents V and the magnetic current M according to equation (21). Once M is known, the fields and field-related parameters may be computed by standard methods. The above solution can be put into a matrix notation as follows: [n]= []' [ m (2 6) where [In] is the excitation vector given by I = - f (2 7) The admittance matrix [Y] can be split into two parts: [Ymn = aperture admittance matrix for equivalent problem b b and [Ymn] = aperture admittance matrix for equivalent problem c. C

-7 - In this manner, equation (26) takes the form: -l n (mn + [mn] m (28) [v] =, ([I+ [YI)i i N (28) with Ymn <fJ' Kx > f g dx'dy' (29) slot Ymn = JJ Km > g dx'dy (30) slot and b 1 _ K = 1 + - F + xx Xx 2 xx zx (31) f 1a2 \ 2 a2 K + - F + F X 2 xx 2 aXa zx (32) k' Bx k dxdz It is important to note that computation of [Y]b involves only regions a, b and computation of [Y]c involves only region c. Therefore we have divided the problem into two parts each of which can be formulated independently. The elements of the two admittance matrices are given by: w/2 b 2 1 1 Ymn = f.(. dy' in n ( dx) /2 T 222 1 - { lx ix k d dx dx sin [k (lx-x' sin kd (l-x)]. ( -' ) 0 0

-8 - I x lx ) + k Jdx Jdx' sin -kd1l-x)] S( 7'-'i) 3I 0 0 (3 3) z=z I = —h y=O and 2 ~c =~64 k d mnn T1 ab k, n m n=l m=0 mw i 0 2b Cos (iY 2 2 k - c a - 2 2 2 kd - a L J Cos( kc) [2 (d + I 2 1x n7c sin - k d - 2 a L j [in (x + x). sian (x + (341) lIn equation (33), the functions D( i-7' - I) and S(I I-r I have the following expressions D r- r'-I)I j axd 2irkd 00 ucosh(uh) + Erduosinh (uh) erduocosh(uh) + usinh(uh) 4 Jo Xp, dX (3 5) S 9-'r II edkd,u cosh (uh) + usinh (uh)

-9 - U,cosh(uh) + e — sinh(uh) + * 5(x 2 ( rd-1) } ucosh(uh) + u sinh(uh) 0I (36) I -1)- 2cos (kdlx) 5(x )] with P = {(x+x I) + P2 = ({[-x+x ') P3 = (X —X') (x -xn ) ]2+ + (Xm-n) 2 ] ' 2 1/2 1/2 ' 2 + y2} 1/2 + y} (37) 2 - (x- x)] m n' = (-x-x ) + (-x )]2+ y }

-10 - 3. Excitation of the Slot As it was mentioned in section 1, the slot is excited by an incident plane wave. This wave can be either parallel or vertically polarized with respect to the plane of incidence. The x-component of the incident magnetic field for both polarizations will be in the form: i " -jkd X sinOt (38) H =H e x 0 with i 2E cos[kd (z+h) cost] sini H (39) C m cos90 - cos(kdhCSt)+ sin (kd hcoset) ~d cos9, se d 1 l 1 2 Eo cos [kd(z+h) cosOt] cosoi H (4 C) rId Tlo coset cos (kdhcosOt) + j sin (kdhcosOt) In view of equations (38)-(40), the elements of the excitation vector [I] are given by the following expression lH cos (kdlXs t esinOxt) cs kl) -jkin I = 2 H e (41) m o

-11 -4. Electric Field Distribution on the Slot-Numerical Results Based on the formulation shown in previous sections, the electric field distribution on cavity-backed slot is shown in figures 1-5 for various slot, substrate and cavity characteristics.

Fiaures —1-3 Er = 2.35 h = O.lko Ca Yc= O.25Xo z= O.124X0 C = 4. 0 rc ga = 0.0 a = 0.-0

Figure 1 Electric Field Distribution on a Cavity-Backed Slot (- = 0o, 0 = 0o).J w 0 u3 O 0,5 E -.- Amplitude -e- Phase Lu 0 (.0 0 (A uL 0 (0 Cu 0.000 0.100 0.200 0.300 Length/Wavelength

Electric Field Distribution on a Cavity-Backed Slot (O= O0, 0 = 00) Figure 2 0.15 0.10 1.4 Lmu 0 (A 0 E.4 0 (A co 0.05 -El- Amplitude -*- Phase 0.00.;j;0.000 0.200 0.400 0.600 Length/Wavelength

Electric Field Distribution on a Cavity-Backed Slot () = 90~, 0 = 60~) Figure 3 *.A LJ 0 0 0 E - Amplitude -- Phase 180 90 LUJ 0 (n LU 0 0 cu -90 -180 0.0000 * — I I 0.000 0.100 0.200 0.300 Length/Wavelength

Figures 4. 5 e = 2.35 r h = O. lXo xc-0. 4X0 Yc= O.25Xo E = 4.0 rc ~L = 1.0 6 = 0. 0 6g = 0.0 Incident Wave = 120 e=o00

Amplitudes for various Zc 's 0.25 Figure 4 -a- A(Zc=.1) -4- A (Zc=.186) - A (Zc=.2714) - A (Zc=.3143) U- A(Zc=.4) 9 -o -0 I 3 E ( 0 5 10 15 20 Length/Wavelength

Phases for various Zc 's 180 90 a-.... Figure 5 -6- P (Zc=.1) -*- P (Zc=. 186) -is- P (Zc=.2714) -*- P (Zc=.3143) I-U P (Zc=.4) w 0 'L 0 4 --- a 2 6 E l~ am A ~ A A - - -A0 - * A O *.... 1 0- 0 0 _on.1 --- —-- m - - - - - - 0- a - - - - - - - 0. 0 O 0- - - i 0 - -4&-",M ID 0 low 0 0 0 -0 i * -0 F-M ---4 4* —m0- MP-," — MF —W - * —4* - 1 RA A) DU 0 1 0 20 Length/Wavelength

E = 2.35 r h = 0 k Yc= O. 25X0 zc= O.125Xo Kor = Oo 6E = 0.0 6~ = 0.0 Incident Wave = 120 e = 00

Amplitudes for various Erc 's, Figure 6 0.3 -0 a ~-2- Amp(2. 0) L4- Amp(4.57) -i-Amp(7.143) 4-Amp(9.714) EE 0.1 -0.000 0.100 0.200030 Length/Wavelength

Phases for various Erc 's Figure 7 180 90 0 -90.0 0 a -C * * *0* * * * * ***..** **- * * -0- Phase(2.0) 4*- Phase(4.57) -0- Phase(7.143) 40- Phase(9.714) - B --- - -— mw —4-" -1 An I.1I. 0.000 0.100 0.200 0.300 Length/Wavelength

Amplitudes for various Erc's Figure 8 0.04 w 0.02 -— '' --- —-------------- Amp(14.86) \ - Amp(17.43) -V- Amp(20.00) E 0.01 0.00 0.000 0.100 0.200 0.300 Length/Wavelength

Phase Eslot / E;' 0 0 0 0 0 0o Io I I 0 0r 0i 00 y ~":3' <f' E! -----— 1 ---- g 0 P P190 -e_ -~9 00 (:a 0- C O c)) 0mm 0~ I'0'-1, I::~:r;3 I~ m(0 o.~co 0o _~

[1] P.B. Katehi, "The Green's Function for a slot cut on the Ground of a Dielectric Substrate," Rad. Lab. Report.

Appendix A Green's Function for a Magnetic Current on the Wall of a Rectangular Cavity Let us begin by considering an x-directed Hertzian magnetic dipole (HMD) at the position (x',y',z') as shown in figure (A.1) With an assumed ejOt time dependency, Maxwell's equations take the form, V XH= jOD -- -- -4 V XE = -jOB - M (1) V -D =0 V B = 0 and D = ~E c (2) B = cH with e~= ~ (1-jtan6 ) (3) c c c = C c(l-jtan ). (4) We can now introduce an electric vector potential function F such that D = Vx F. (5) When we substitute equation (5) into (la), the following relation is derived Vx[ - jCO]. (6)

-A. 2 - Equation (6) is equivalent to -4 -- H = jCOF + V7. (7) From equations (7) and (lb), we can get the following result -4 - - -4 -- -2 -VXVXF=V (V F)-V2F= kCF - jCECVm - EM. (8) Since equation (5) alone cannot specify F, we can select its divergence so that - - V-F = -jEO (9) m. Equations (8) and (9) then give (2 + k) F = ~M. (10) The electric and magnetic fields, in terms of F, are given as follows - 1 - E = V xF (11.) C 2 H = ( 02 - V V-. L} CJ (12) with 2 2 k = C For a magnetic current M with more than one component, the electric potential can be consider as a dyadic function of the form AA AA AA F = xx F + yy F + zz F. (11) xx yy zz For the case of Hertzian dipoles F is being considered as the Green's function of the problem and satisfies the following differential

equation + k G = I (12) where I is the unit dyadic. For the problem of figure A.1, G =G =0 yy zz and E =0 x 1 aGxx Y (13) c 1 xx E = - z ay C Also, the magnetic field is given by 2 H = J k + G sx x C 2 c H = - a G (14) Y i) ayax xx H.= G z 0~ L azax xx c c The boundary conditions for this problem are: I E = 0 z = 0 x = 0,a Y I E = 0 z = 0,a y = 0,b z (15) E = 0 z = c x = 0,a II = 0 x = 0,a y = 0,b z

The solution to equation (12) which satisfies the above boundary conditions is in the form: region I I c 4 nx' in y' G E - Z sin os xx ab n m a b n=l m=0 (nnx mny ~ sin cos - a (nt b c (16) cos[kc-z)] cos (k z) k sin(kzc) z z region II FcI nKx ' mny' G f ~ ~ e E sin i cos ' xx ab n a b n= 0 m=0 ( nKxx mIy\ ~ sin cos. (a c b (17) cos( kz') cos[ k(c-z) z z J k sin(kc) with 2 2 k2 k2 _ (18) Using the above expressions and for the case of a magnetic dipole Mx at z' = 0, the x-component of the magnetic field in the cavity

takes the form HX = 2II 2y G3X (xfyiz;x',y') MX(x',y )dxcldy'. (19) "'i L k axj

Print file "SLOT. FTN" >.................... Slot.ftn................ This program evaluates the elements of the admittance matrix in the the following problem: - "Scattering by two slots on infinite dielectric substrate with cavities at the back." This program is good for any substrate thickness h, er and any dimensions of the cavities. SUBROUTINE SLOT IMPLICIT REAL*8 (A-H,O-Z) REAL*4 CONST,GSK,GS1S2K COMPLEX ZS,ZS1S2,CI,BMATR,CUR COMMON/CTAIL/S1 (4,205,7),D1(4,205,7),D2(4,205,7), *T1(4,205,7),T2(4,205,7),T3 (4,205,7) T4 (4205 COMMON/COMP/ZS(50),ZSlS2(350),NS,NS1S2 COMMON/OUT/GS(50),GS1S2(350) COMMON/MAT/PLI,AI,TI,V(3),IY COMMON/PUT/SSJO(250,7),SAJO(250,7),YSIN,YCOS COMMON/ADON/DIST(250,7,10),RCOE(20,250,7,10),AX,SERS(5),SERA(5), *DARG(10, 4), S(10,2),WREAL,NSER, NMAX(7) COMMON/DAT/ER,H,T,DLX,XC,YC,ZC,ERC,RMC,XO,YO,A,TPI,TPI2,PI,W,E1, *E2,EER,AK0,AK,AKK,FA,OFFSET(7),WDELTA,OFFLIM,ERROR,THI,FI, *NS1,NS2,NSS2,NOFF COMMON/DATT/COAL(20),POINT(20),CN(51),BM(51),POLTM(20), *POLTE(20),AM(41), DM(41),POLES (40),VXXM(20),VZXM(20),VZXE (20), *BPOINT (10), BCOAL(10), MPOINT,NPOINT,NKO, MA,NTM,NTE,NKOK,IFIRST COMMON/COEF/RX, XX, RZ, XZ, FRX, FRZ, F1X, F1Z COMMON/IOFF/INS, INS1S2 COMMON/B01/BJ, BJ1 CALL DATA SLOT Subroutine POLES evaluates the poles of the Green's function and orders them according to their magnitude IFIRST=0:dominant mode is a TM wave 1:dominant mode is a TE wave 2:only one TM wave CALL SPOLES PRINT *,'SPOLES has run' H=H*DSQRT(ER) CI=(0.00,1.00) NS=NS1 IF (NS1.LT.NS2) NS=NS2 Page 1

Print file "SLOT.F27JUF Pa ge 2 C C C C C MS=NS IF (NOFF.EQ.1) GO TO 50 NSlS2=NS2+NSS2-1 MS1S2=NS 152 IF (NSlS2.GT.200) NS1S2=200 50 IF(NMAX(INS).LE.(NS+2)) NMAX(INS)=NS+2 IF (NOFF.EQ.1) GO TO 51 IF(NMAX(INS1S2).LE. (NS1S2+2)) NMAX(INSlS2)=NSlS2+2 51 ADL=AKK*DLX YS IN=D SIN (AD L) YCOS=DCOS (ADL) For the normalization of the current along the y axis CVON=W*PI/2.DO Computation of lamda-integration limits between 0 and A CALL LIMIT Evaluation of the Green's function at different points in the interval [0,A]. The Bessel function has been excluded CALL GREEN Evaluation of the tail contribution (from a to infinity) CALL TAIL C"ONST=(1l.DO/CVON) *DSQRT (ER) /(480.DO* (PI**3) *YSIN*YSIN) WRITE (6,9) 9 FORMATf ( --- —--------------------- Zs --- —------------------- f WRITE(6,10) MS 10 FORMAT(11X,14) DO 11 K=1,MS ZS (K) =ZS (K) *CONST GSK=REAL (GS (K) )*CONST WRITE (6,30) K,ZS(K),GSK 30 FORMAT(lXI4,2XO'ZS=O',El4.7,2XEl4.7,2X, * ~'GSK=',E14.7) ZS (K) = (ZS (K) +GSK) *CI WRITE (6,12) ZS (K) 12 FORMAT(68XE14.7,1XE14.7) 11 CONTINUE IF (NOFF.EQ.1) GO TO 52 WRITE (6,15) 15 FORMAT(' -------------- S12 --- —----------— ZSS WRITE (6,10) MS1S2 DO 16 K=1,MSlS2 ZSlS2 (K) =ZSlS2 (K) *CONST GS1S2K=REAL (GS 1S2 (K) ) *CONST WRITE (6, 32) K, ZS1S2 (K),fGS1S2K 32 FORMAT(lXI4,2Xf,ZS1S2=',E14.7,,2X,,El4.7,,2X, *f'GSlS2K-', E14.7) WRITE (6, 12) ZS1S2 (K) ZSlS2 (K) =(ZS1S2 (K) +GSlS2K) *CI 16 CONTINUE 52 CONTINUE 100 0 CONTINUE RETURN END

Print file "SLOT.FTN" Page 3 C.................................................................... C...................................................... C This subroutine evaluates the limits of integration in C the interval [0,A]. C Specifically: C 1) It divides the interval [0,k0] to 10 equal C subsections and then apply fixed-point Gaussian C Quadrature C 2) It divides the interval [k0,k] into so many C subsections as the number of poles and in C such a way that each subsection includes one C pole only away from the ends of the subsection C 3) It divides the interval [k,A] into 20 equal C subsections and then apply fixed-point Gaussian C Quadrature C..................................................................... SUBROUTINE LIMIT IMPLICIT REAL*8 (A-H,O-Z) EXTERNAL WSPE,WTPE,WSPM COMMON/DAT/ER,H,T,DLX,XC,YC,ZC,ERC,RMC,XO,Y0,A,TPTPI,TPI2,PI,W,El, *E2,EER,AKO,AK,AKK,FA, OFFSET (7),WDELTA,OFFLIM,ERROR,THI,FI, *NS1,NS2,NSS2,NOFF COMMON/DATT/COAL(20),POINT(20),CN(51),BM(51),POLTM(20), *POLTE(20),AM(41),DM(41),POLES(40),VXXM(20),VZXM(20),VZXE(20), *BPOINT(10),BCOAL(10),MPOINT,NPOINT,NKO,MA,NTM,NTE,NKOK,IFIRST ------------------------------------------------— + Step 1: Evaluation of vector CN I it gives the end points of the I intervals considered in (0,kO) I ----------------------------------------------— + DELTA=AKO/FLOAT(NKO) CN(1)=0.DO DO 1 I=1,NKO CN(I+1)=DELTA*FLOAT(I) 1 CONTINUE ------------------------------------------------— + Step 2: Evaluation of vector BM it gives the end points of the intervals considered in (k,A) ---------------------------------------------— + DELTA= (A/DSQRT (EER)-AK)/FLOAT (MA) BM(1)=AK DO 2 I=1,MA BM(I+1)=DELTA*FLOAT(I)+AK 2 CONTINUE ------------------------------------------------— + Step 3: Evaluation of the vectors AM,DM "AM" gives the end points around the TM poles "DM" gives the end points around the TE poles IFIRST= 2 only one TM pole 1 TEO<TMO 0 TMO<TEO ----------------------------------------------— + AM(1)=AKO DM(1)=AKO NMAX=NTE+NTM- 1 IF (IFIRST.EQ.2) GO TO 3 DO 4 I=1,NMAX AM(I+1)=(POLES (I+1)+POLES(I))/2.DO DM(I+1)=AM(I+1)

Print file "DATA. SLOT" Page 1 C C ---- Dielectric constant --- C 2.35 C C ---- Substrate Thickness --- C 0.1 C C ---- Conductor Thickness --- C 0.0001 C C ---- Dimensions of the Cavity ---- C 0.4 0.25 0.124 C C ---- Dielectric Permittivity and Permeability of the Medium in the Cavity ---- C 4.00 1.00 0.00 0.00 C ---- Slot width ---- 0.1 ---- Position of the Center of the Leading Edge of the Slot ---- 0.072796 0.075 Subsection Length 0.022 ---- Lower Limit of the Tail Contribution:100.0 Number of Points on Each Slot ---- 23 0 Number of Offsets ---- 1 Offset Between the Two Slots ---- 0.0 ---- Longitudinal Displacement Between the Two Slots ---- 0 ---- Theta angle 0.0 ---- Fi angle ----

Print file "DATA. SLOT" Page 2 0.0

Print file "MAIN SLOT.FTN"Pg Page I C....................................... o MAINSLOT.FTN o This program evaluates the reflection coefficient of C cavity-backed slot as a function of the characteristics C of the cavity and the incident wave C............................................................... IMPLICIT REAL*8 (A-HO-Z) %INCLJDE '/SYS/INS/CAL. INS.FTN' INTEGER*4 CPU SECONDS INTEGER*2 TIMEDATEREC(6) COMPLEX ROREFLRRATIOCURBMATRZSZS1S2 CHARACTER*10 PARl, PAR2 LOGICAL THFLAG1, FIFLAGi, XCFLAG1, YCFLAG1, ZCFLAG1, ERCFLAG1 LOGICAL THFLAG2, FIFLAG2, XCFLAG2, YCFLAG2, ZCFLAG2, ERCFLAG2 COMMON/CTAIL/Sl (4,205,7), DJ (4,205,7),D2 (4,205,7), *T1 (4,205, 7),T2 (4,205,7),T3 (4,205, 7),T4 (4,205, 7) COMMON/COMP/ZS (50), ZSlS2 (350),NSNS1S2 COMMON/OUT/GS(50),GS1S2 (350) COMMON/MAT/PLI, AI,,TI, V(3), IY COMMON/PUT/SSJO (250,7),SAJO (250,7),YSINYCOS COMMON/ADON/DIST (250,7,l0), RCOE(20,250, 7, 10), AX, SERS (5), SERA(5), *D'ARG(l0,4),S(10Q,2), WREALNSERNMAX(7) COMMON/DAT/ER, H, TDLX, XC, YC, ZC, ERC, RMC, XO, YOfA, TPI, TPI2, PI,W, El, *E2rEERfA0,AK,(K, FAOFFSET(7),WDELTAOFFLIMERRORTHIFIf *NSlNS2,NSS2,NOFF COMMON/DATT/COAL(20),POINT(20),CN(51),BM(51),POLTM(20), *POLTE(20),AM(41),DM(41),POLES (40),VXXM(20),VZXM(20),VZXE(20), *BPOINT(10),BCOAL(10),MPOINTNPOINTNKOMANTMNTENKOKIFIRST COMMON/COEF/RX, XX, RZ, XZ, FRX, FRZ, F1X, F1Z COMMON! IOFF/INS, INS1S2 COMMON/MAN/BMATR(260,260), IA(260), fB (260) COMMON/PAT/CUR (260) COMMON/BO 1/BJO, BJl COMMON/LOSS/TLOS E,TLOS M OPEN (UNIT=0lFILE='OUT.PLOT',STATUS='OLDf) OPEN (UNIT=05,FILE='DATA.SLOT',STATUS='OLD') OPEN (UNIT=06,FILE='OUT.SLOT' STATUS=fOLD') CALL PROC1_$GET CPUT(CLOCK) CPUSECONDS=CAL_$CLOCKTOSEC(CLOCK) THFLAG1=.FALSE. FIFLAG1=.FALSE. XCFLAG1=.FALSE. YCFLAG1=. FALSE. ZCFLAG1=.FALSE. ERCFLAGl=.FALSE.

Print file "sMAIN SLOT. FTN" Page 2 C C THFLAG2=.FALSE. FIFLAG2-.FALSE. XCFLAG2-.FALSE. YCFLAG2=.FALSE. ZCFLAG2=.FALSE. ERCFLAG2=.FALSE. PRINT *, 'You can plot the reflection coeff. as a function of:' PRINT *, 'THETA, FI,XC, YC, ZC,ERC' PRINT *, 'The reflection coefficient will vary as a function of:' PRINT *, 'First parameter:' READ (*,1) PAR1 PRINT *, 'Second parameter:' READ (*,1) PAR2 1 FORMAT(A10) PRINT *, 'What is the maximum value READ (*,2) PAR1MX 2 FORMAT(E14.7) PRINT *, 'What is the minimum value READ (*,2) PAR1MN PRINT *, 'How many points should be PRINT *, '(number of points >=1)' READ (*,23) NPAR1 PRINT *, 'What is the maximum value READ (*,2) PAR2MX PRINT *, 'What is the minimum value READ (*,2) PAR2MN PRINT *, 'How many points should be PRINT *, '(number of points >=1)' READ (*,23) NPAR2 23 FORMAT (I5) of ', PAR1 of ',PAR1 considered' of ', PAR2 of ',PAR2 considered' Subroutine DATA reads the given values for variables different CALL DATA PRINT *,'Subroutine DATA run' Subroutine slot evaluates the parts of the elements of the impedance matrix which are coming from the dielectric substrate WREAL=W W=W*(1.DO+2.DO*WDELTA/W) CALL SLOT PRINT *,'Subroutine SLOT run' DELTA1=0.DO DELTA2=0.DO IF (NPAR1.GT.1) DELTA1= (PAR1MX-PAR1MN) /FLOAT (NPAR1-1) IF (NPAR2. GT. 1) DELTA2= (PAR2MX-PAR2MN)/FLOAT (NPAR2-1) NTHI=-1 NFI=-i NXC=-1 NYC=-1 NZCM-1 NERC=- 1 IF ((PARl.EQ.'THETA').OR.(PAR1.EQ.'theta')) THFLAG1=.TRUE. IF ((PAR2.EQ. ' THETA').OR.(PAR2.EQ. ' theta' ) ) THFLAG2. TRUE. IF ((PAR1.EQ.'FI').OR.(PAR1.EQ.'fi')) FIFLAG1-.TRUE.

Print file "MAIN SLOT.FT"' Pg3 Pa ge 3 IF IF IF IF IF IF IF IF IF ((PAR2. EQ. 'FI' ). OR. (PAR2.EQ.' fi')) ((PARl. EQ. XC').OR. (PAR1.EQ. xc')) ((PAR2. EQ. XC').OR. (PAR2.EQ. xc')) ((PAR1.EQ.'YC').OR.(PAR1.EQ.'yc')) ((PAR2. EQ. ' YC'). OR. (PAR2.EQ. yc' ) ((PAR1.EQ. 'ZC').OR.(PAR1.EQ. 'zc')) ((PAR2. EQ. 'ZC' ). OR. (PAR2. EQ. ' zc' )) ((PARl. EQ. 'ERC' ). OR. (PARl.EQ.' erc' ((PAR2.EQ.'ERC').OR.(PAR2.EQ. 'ERC' FIFLAG2=. TRUE. XCFLAG1=.TRUE. XCFLAG2=. TRUE. YCFLAG1=. TRUE. YCFLAG2=. TRUE. ZCFLAG1J=. TRUE. ZCFLAG2=. TRUE. )) ERCFLAG1=. TRUE. )) ERCFLAG2=.TRUE. IF (THFLAG1) THEN THIMX=PAR1MX*PI/180.DO THIMN=PAR1MN*PI/180.DO DELTHI=DELTA1*PI/180.DO NTH I=NPAR1 IDTHI~1 END IF IF (THFLAG2) THEN THIMX=PAR2MX*PI/180.DO THIMN=PAR2MN*PI/ 18. DO DELTHI=DELTA2*PI/180.DO NTHI=NPAR2 IDTHI=2 END IF IF (FIFLAG1) THEN FIMX=PAR1MX*PI/180.DO FIMN=PAR1MN*PI/180.DO DELFI=DELTA1*PI/180.DO NF I=NPAR1 IDFI=1 END IF IF (FIFLAG2) THEN FIMX=PAR2MX*PI/ 18. DO FIMN=PAR2MN*P I/ 18. DO DELFI=DELTA2*PI/ 180 DO NF I=NPAR2 IDFI=2 END IF IF (XCFLAG1) THEN XCMX=PAR1MX XCMN=PAR1MN DELXC=DELTAJ NXC=NPAR1 ID XC=1 END IF IF (XCFLAG2) THEN XCMX-PAR2MX XCMN=PAR2MN DELXC=DELTA2 NXC=NPAR2 IDXC=2 END IF IF (YCFLAG1) THEN YCMX=PAR1MX YCMN-PAR1MN DELYC=DELTA1 NYC=NPAR1 IDYC=1 END IF IF (YCFLAG2) THEN YCMX-PAR2MX YCMN=PAR2MN DELYC=DELTA2 NYCNUPAR2 ID YC-2 END IF

Print file IIMAIX SLOT.FTN"a Page 4 IF (ZCFLAG1) THEN ZCMX=PAR1MX ZCMN=PARlMN DELZC=DELTA1 NZC=NPAR1 IDZC=l END IF IF (ZCFLAG2) THEN ZCMX=PAR2MX ZCMN=PAR2MN DELZC=DELTA2 NZC=NPAR2 IDZC=2 END IF IF (ERCFLAG1) THEN ERCMX=PAR1MX ERCMN=PAR1MN DELERC=DELTA1 NERC=NPAR1 IDERC=1 END IF IF (ERCFLAG2) THEN ERCMX=PAR2MX ERCMN=PAR2MN DELERC=DELTA2 NERC=NPAR2 IDERC=2 END IF C WRITE (1,100) PAR1,PAR2 100 FORMAT (5X, Al O, 5X, Al 0, 5X,f'Rel. Refl. Coef f icient-Real',5X, *IRel.Refl.Coefficient-Imaginary'//) C SLENG=(NS1+1) *DLX IF (NXC.EQ.-1) THEN XCMN=XC NXC=1 END IF DO 3 IXC=1,NXC XC=XCMN+FLOAT (IXC-1) *DELXC XO=(XC-SLENG)/2.DO IF(IDXC.EQ.1) WPAR1=XC IF(IDXC.EQ.2) WPAR2=XC IF (NYC.EQ.-1) THEN YCMN=YC NYC=1 END IF DO 4 IYC=1,NYC YC=YCMN+FLOAT (IYC-1) *DELYC YO=YC/2.DO IF(IDYC.EQ.1) WPAR1=YC IF(IDYC.EQ.2) WPAR2=YC IF (NZC.EQ.-1) THEN ZCMN=ZC NZC=l END IF DO 5 IZC=1,NZC ZC=ZCMN+FLOAT(IZC-1) *DELZC IF(IDZC.EQ.1) WPAR1=ZC IF(ID ZC.EQ.2) WPAR2=ZC IF (NERC.EQ.-1) THEN ERCMN=ERC NERC=l END IF DO 6 IERC=1,NERC ERC=ERCMN+FLOAT(IERC-1) *DELERC

Print file "MAIN SLOT. FTN" IF(ID ERC.EQ.1) WPAR1=ERC IF(ID ERC.EQ.2) WPAR2=ERC CALL CAVINV(NOR) PRINT*, 'Subroutine CAV INV run' IF (NTHI.EQ.-1) THEN THIMN=THI NTHI=1 END IF DO 7 ITHI=1,NTHI THI=THIMN+FLOAT (ITHI-1) *DELTHI IF(ID THI.EQ.1) WPAR1=THI*180.DO/PI IF(ID THI.EQ.2) WPAR2=THI*180.DO/PI IF (NFI.EQ.-l) THEN FINM=FI NFI=1 END IF DO 8 IFI=1,NFI FI=FIMN+FLOAT (IFI-1) *DELFI IF(ID FI.EQ.1) WPAR1=FI*180.DO/PI IF(ID FI.EQ.2) WPAR2=FI*180.DO/PI CALL MULTEXC(NOR,RO) PRINT *,'Subroutine MULTEXC run' CALL REFLEC(WREAL, R, REFL) PRINT *,'Subroutine REFLEC run' R RATIO=REFL/RO WRITE (1,10) WPAR1,WPAR2,R RATIO 10 FORMAT(5X,E14.7,5X,E14.7,5X, (E14.7, * 1X,E14.7),5X,(E14.7,1X,E14.7)) 8 CONTINUE 7 CONTINUE 6 CONTINUE 5 CONTINUE 4 CONTINUE 3 CONTINUE WRITE (6,9000) CPUSECONDS 9000 FORMAT(1X,'CPU TIME=',I10,'SECS') STOP END The name of this subroutine is DATA and gives all the data used by the main program and the other subroutines. SUBROUTINE DATA IMPLICIT REAL*8 (A-H,O-Z) COMMON/DAT/ER, H,T, DLX, XC,YC, ZC, ERC,RMC,XO, YO, A, TPI, TPI2, PI, W, El, *E2,EER, AK0, AKAKK, FA,OFFSET(7),WDELTA,OFFLIM,ERROR,THI,FI, *NS1, NS2,NSS2,NOFF COMMON/IOFF/INS, INS1S2 COMMON/LOSS/TLOS E,TLOS M PI=3. 141592653589D0 Page 5 ---- Dielectric constant ---

Print file "MAIN SLOT.FTN" READ (5,1) ER 1 FORMAT (///6X,D16.9) EER=ER AER=DSQRT(ER) WRITE (6,2) ER 2 FORMAT(10X,'Dielectric Constant of the Substrate'/10X,E14.7/) C C ---- Substrate Thickness --- C READ (5,1) H WRITE (6,3) H 3 FORMAT(10X,'Substrate Thickness'/10X,E14.7/) ---- Conductor Thickness --- READ (5,1) T WRITE (6,4) T 4 FORMAT(10X,'Conductor Thickness'/10X,E14.7/) T=T*AER -- Dimensions of the Cavity ---- READ (5,1) XC READ (5,10) YC READ (5,10) ZC 10 FORMAT(6X,D14.7) WRITE (6,5) XC,YC,ZC 5 FORMAT(10X,'Dimensions of the Cavity'/10X,'XC=',E14.7/ *10X,'YC=',E14.7/10X,'ZC=',E14.7/) XC=XC*AER CYC=YC*AER ZC=ZC*AER --- Relative Permittivity and Permeability of the Cavity ---- READ (5,1) ERC READ (5,10) RMC READ (5,10) TLOS E READ (5,10) TLOS M WRITE (6,6) ERC,RMC,TLOS E,TLOS M 6 FORMAT(10X,'Dielectric of the Cavity'/10X,'ERC=',E14.7,2X,'RMC=' *,E14.7/10X,'TLOS E='E14.7,2X,'TLOS M=',E14.7) -- Slot width ---- READ (5,1) W WRITE (6,7) W 7 FORMAT(10X,'Slot Width'/10X,E14.7/) W=W*AER ---- Position of the Center of the Slot ---- READ (5,1) XO READ (5,10) YO WRITE (6,8) XO,YO 8 FORMAT(10X,'Position of the Center of the Slot'/lOX,'XO=', *E14.7/10X,'YO=',E14.7/) XOXO*AER YO=YO*AER ---- Subsection Length --- READ (5,1) DLX WRITE (6,9) DLX 9 FORMAT(10X,'Subsection Length'/10X,E14.7/) Page 6

Print file "MAIN SLOT.FTN" C Page 7 C C C C C. --- Lower Limit of the Tail Contribution ---- READ (5,1) A WRITE (6,11) A 11 FORMAT(10X,'Lower Limit of Tail Contribution'/10X,E14.7/) ---- Number of Points on Each Slot ---- READ (5,20) NS1 READ (5,21) NS2 20 FORMAT(///6X,14) 21 FORMAT(6X,14) WRITE (6,12) NS1,NS2 12 FORMAT (10X,'Number of Points on Each Slot'/10X,'NS1=',I4/ *10X, 'NS2=', I4) ---- Number of Offsets ---- READ (5,20) NOFF WRITE (6,15) NOFF 15 FORMAT(10X,'Number of Offsets'/10X,I4/) --- Offset Between the Two Slots --- READ (5,1) OFFSET(1) WRITE (6,13) OFFSET(1) 13 FORMAT (10X,'Offset Between the Two Slots'/10X,E14.7/) OFFSET (1)=OFFSET (1) *AER DO 16 I=2,NOFF READ (5,10) OFFSET(I) WRITE (6,17) OFFSET(I) 17 FORMAT(10X,E14.7/) OFFSET(I)=OFFSET(I)*AER 16 CONTINUE Order of Offsets ---- INS=1 INS1S2=2 ----- Longitudinal Displacement Between the Two Slots ---- READ (5,20) NSS2 WRITE (6,14) NSS2 14 FORMAT(10X,'Longitudinal Displacement Between the Two Slots'/ *lOX,I4/) -Theta angle ---- READ (5,1) THI WRITE (6,23) THI 23 FORMAT(10X,'Theta Angle'/10X,E14.7/) ---- Fi Angle ---- READ (5,1)FI WRITE (6,22) FI 22 FORMAT(10X,'Fi Angle'/10X,E14.7/) RETURN END

Print f-le "SZLOT.TPTN" Page 8 COMMON IOFF/INS, INS1S2 C CI=(O.O,1.O) NCON=O X=AUP -ALOW Y=AUPJ+ALOW AK02=AKO*AKO AK2=AK*AK ER11.DO-ER NPOL=N IF (IFIRST.EQ.2) NPOL=1 IF (IAD.GT.2) GO TO 1 ALI=O05D0*(TI*X+Y) GCONX=AI*X*O 5DO FCONX=GCONX GCONZ=GCONX*ER1 FCONZ=FCONX AIMA1. DO CALL GREI(ALI,O.DO.O.DO, IADO.DO) GO TO 10 1:F (IAD.NE.3) GO TO 2 ALI-0 5D0 O*(TI*X+Y) XTM-POLTM (NPOL) TMTM-(2.DO*XTM-Y) /X GCONX=AI/ (TI-TMTM) GCONZ=GCONX*ER1 FCONX=GCONX FCONZ=FCONX AIMA=0.DO CALL GREI(ALfIXTM,0.DO, IAD, O.DO) GO TO 10 2 IF (IAD.NE.4) GO TO 3 ALI=POLTM (N) TM= (2. DO *ALIY) /X GCONX=-AI/(TI-TM) GCONZ=GCONX*ER1 FCONX=.DO FCONZ=0.DO AIMA=0.DO RX=VXXM (N) RZ=VZXM (N) GO TO 10 3 IF (IFIRST.EQ.2) GO TO 5 IF (IAD.NE.5) GO TO 4 ALI=0.5D* O(TI*X+Y) XTE=POLTE(N) TMTE=(2.DO*XTE-Y) /X GCONX=AI*X*0 5DO GCONZ=AI*ER1/ (TI-TMTE) FCONX=GCONX FCONZ=FCONX AIMA=0 DO CALL GREI (ALl, 0.DO, XTE, IAD, TMTE) GO TO 10 4 IF (IAD.NE.6) GO TO 5 NCON=6 ALI=POLTE (N) TM= (2. DO *ALIY) /X GCONX=O.DO GCONZ=-AI*ER1/ (TI-TM) FCONX=O.DO FCONZ=O.DO AIMA=0.DO RZ=VZXE (N) GO TO 10 5 IF (IAD.NE.7) GO TO 6

Print file "SLOT.FTN" Page 9 ALI=0.5D0*(TI*X+Y) GCONX-AI*X*0.5D0 GCONZ=GCONX*ER1 FCONX=GCONX FCONZ=FCONX AIMA=0.DO CALL GREI(ALI, 0.D0,0.D0,IAD, 0.D0) GO TO 10 6 NCON=8 ALI=POLTM(N) TM=(2.D0*ALI-Y)/X FCONX=0.DO FCONZ=0.DO AIMA=0.DO RX=VXXM(N) RZ=VZXM(N) GO TO 28 10 CONTINUE GXXR=GCONX*RX-FCONX*FRX GXXX=AIMA*GCONX*XX GZXR=GCONZ *RZ-FCONZ *FRZ GZXX=AIMA*GCONZ *XZ 27 CONTINUE VARX-AK2*GZXR VARZ-AKK*(GXXR-GZXR) GXXR-VARX GZXR=VARZ VARX=AK2 *GZXX VARZ=AKK*(GXXX-GZXX) GXXX=VARX GZXX=VARZ PLI=ALI CALL ADONIS DO 13 K=1,NS S1=REAL(GXXR*SSJO(K,INS)+GZXR*SAJ0(K,INS)) S2=REAL(GXXX*SSJO(K,INS)+GZXX*SAJO(K, INS)) ZS (K)=ZS (K) +S1-CI*S2 13 CONTINUE DO 20 K=1,NSlS2 Sl=REAL(GXXR*SSJ0(K,INSlS2)+GZXR*SAJO(K,INSlS2)) S2=REAL (GXXX*SSJO (K, INSlS2) +GZXX*SAJO (K, INSlS2)) ZSlS2(K)=ZSlS2(K)+Sl-CI*S2 20 CONTINUE 28 IF (NCON.EQ.0) GO TO 24 GCONX1=0.0 GCONX2=0.0 GCONZ1=ER1*DLOG((1.DO-TM)/(1.DO+TM)) GCONZ2=ER1*PI IF (NCON.EQ.6) GO TO 29 GCONX1=GCONZ1/ER1 GCONX2=GCONZ2/ER1 29 CONTINUE GXXR=GCONX1*RX GXXX=GCONX2 *RX GZXR=GCONZ1*RZ GZXX=GCONZ2*RZ 25 CONTINUE NCON-0 GO TO 27 24 CONTINUE RETURN END This subroutine evaluates the integrand of the green's

Print file "SLOT. FTN" Page 4 4 CONTINUE AM (NMAX+2) =AK DM (NMAX+2) =AK IF (IFIRST.EQ.1) GO TO 5 DM(NMAX+1) =AM (NMAX+2) DO 6 I=1,NMAX DM(NMAX-I+1) =AM (NMAX-I+2) 6 CONTINUE GO TO 7 5 AM(NMAX+1) =DM (NMAX+2) DO 8 I=1,NMAX AM (NMAX- I+1) =DM (NMAX-I+2) 8 CONTINUE GO TO 7 C 3 DELTA= (AK-AKO) /FLOAT (NKOK) AM(1)=AKO DO 9 I=1,NKOK AM (I+1) =DELTA*FLOAT (NKOK) +AKO 9 CONTINUE 7 CONTINUE -----------------------------------------------— + Step 4: evaluation of vectors VZXE I -----—. --- —------------------------------------— + IF (IFIRST.EQ.2) GO TO 10 DO 11 I=1,NTE ARG=POLTE(I) VZXE (I) =HZXE (ARG) 11 CONTINUE 10 CONTINUE --------------------------------------------— + Step 5: evaluation of vector VXXM,VZXM I. --- —. --- —-----------------------------------— + DO 12 I=1,NTM ARG=POLTM(I) VXXM(I)=GXXM(ARG) VZXM (I) =GZXM (ARG) 12 CONTINUE RETURN END This subroutine evaluates the values of the integrand of the Green's function at different points in the interval [0,A]. Then it evaluetes the space integrals of the Bessel function at the same points and multiply these values with the corresponding values of the Green's function. Finally, it multiplies these products with known coeffic. and it adds them up. This way, the moments'-method space integrals of the first part of the Green's function are evaluated and are stored in the complex vectors ZS,ZS1S2 SUBROUTINE GREEN IMPLICIT REAL*8 (A-H,O-Z) COMPLEX ZS,ZS1S2,CI COMMON/COMP/ZS (50), ZS1S2 (350), NS, NS1S2 COMMON/MAT/PLI, AI, TI,V(3), IY COMMON/PUT/SSJO (250,7),SAJO (250,7), YSIN,YCOS COMMON/ADON/DIST(250,7, 10),RCOE(20,250,7, 10),AX,SERS(5),SERA(5), *DARG(10, 4), S (10,2), WREAL, NSER, NMAX(7)

Print file "SLOT. FTN" Page 5 COMMON/DAT/ER, H, T, DLX, XC, YC, ZC, ERC, RMC, XO, YO, A,TPI,TPI2, PI,W,El, *E2,EER, A AK, AKAKK,FA,OFFSET(7), WDELTA,OFFLIM, ERROR, THI,FI, *NS1,NS2,NSS2,NOFF C COMMON/DATT/COAL(20),POINT(20),CN(51),BM(51),POLTM(20), *POLTE(20),AM(41),DM(41),POLES(40),VXXM(20),VZXM(20),VZXE(20), *BPOINT(10),BCOAL(10),MPOINT,NPOINT,NKO,MA,NTM,NTE,NKOK,IFIRST C COMMON/COEF/RX, XX, RZ, XZ, FRX, FRZ, F1X, F1Z C COMMON/IOFF/INS, INS1S2 C C --- —--------- ---------------------------------— + C Evaluation of the coefficients for the C FF's functions C --- —------------------------------------------— + FlX=1.DO FlZ=(1.D0-ER)/( 1.D0+ER)*(1.D0+E2)*(1.D0+0.5D0*El)) CALL ARIS DO 1 I=1,NPOINT IY=I AI=COAL(I) TI=POINT(I) evaluation of intervals 1 and 2 IAD=1 DO 2 N=1,NKO AUP=CN(N+1) ALOW=CN(N) CALL FUNCT (IAD,AUP,ALOW,N) 2 CONTINUE evaluation of intervals 3 and 4 NTTM=NTM IF (IFIRST.EQ.2) NTTM=NKOK DO 3 IAD=3,4 IFD=0 DO 4 N=1,NTTM IFD=IFD+1 AUP=AM(IFD+1) ALOW=AM(IFD) CALL FUNCT(IAD,AUP,ALOW,N) IFD=IFD+1 4 CONTINUE 3 CONTINUE IF (IFIRST.EQ.2) GO TO 9 evaluation of the intervals 5 and 6,9,11 DO 5 IAD=5,6 IFD=0 DO 6 N=1,NTE IFD=IFD+1 AUP=DM(IFD+l) ALOW=DM(IFD) CALL FUNCT(IAD,AUP,ALOW,N) IFD=IFD+1 6 CONTINUE 5 CONTINUE 9 CONTINUE evaluation of the interval 7

Print file FVSLOTFTWPg Page 6 IAD=7 DO 7 N=1,MA AUP=BM(N+1) ALOW=BM (N) CALL FUNCT(IADAUPALOWN) 7 CONTINUE 1 CONTINUE evaluation of the intervals 8,10 IAD=8 IFD=0 DO 8 N=1,NTM IFD=IFD+1 AUP=AM(IFD+1) ALOW=AM(IFD) CALL FUNCT(IADAUPALOWN) IFD=IFD+1 8 CONTINUE RETURN END Functions GXXMGZXMHZXE These functions evaluate the residues from the different poles FUNCTION GXXM(X) IMPLICIT REAL*8 (A-HO-Z) COMMON/DAT/ERH,T,DLX, XC, YC, ZC, ERC, RMC, XO, YO, A, TPI, TPI2, PI, W, E1, *E2,EERAKO, AKjAKKFAOFFSET(7),WDELTAOFFLIMERROR, THIFI, *NS1,NS2,NSS2,NOFF X2=X*X AKO2=AKO *?AKO AK2=AK*AK RM=DSQRT (AK2-X2) RMO=DSQRT (X2-AKO2) RMH=RM*H RMOH=RMO *H RMT=RM* (-H+T) SXN=RM*DCOS (RMT) -ER*RMO *DS IN (RMT) SXD= (ER+RMOH) * (RM/RMO) *DCOS (RMH) + (1.DO+ER*RMOH) *DSIN( RMH ) GXXM= S XN / S XD RETURN END FUNCTION GZXM(X) IMPLICIT REAL*8 (A-HO-Z) COMMON/ DAT / ER, H, T, DLX, XC, YC, ZC, ERC, RMC, XO, Y O, A, TP I, TP I 2, PI, W, E 1, *E2,EERAAKO,fAKJAjK, FAOFFSET(7),WDELTAOFFLIMERRORTHIFI, *NS1,NS2,NSS2,NOFF X2=X*X AKO2-AKO *AK 0 AK2=AK*AK RM=DSQRT (AK2-X2) RMO=DSQRT (X2-AKO2) RMH=RM *H RMOH-RMO *H RMT-RM*T CST-DCOS(RMT)

Print file "SLOT.FTN" Page 7 CSH=DCOS (RMH) SNH=DSIN (RMH) SXN=RM*CST SXD= (RM*CSH+RMO*SNH) * ( (ER+RM0H) *CSH/PRMO+ (1.DO+ER*RM0H) *SNH/RM) GZXM=SXN/ SXD RETURN END C C C........................................................................ C FUNCTION HZXE(X) IMPLICIT REAL*8 (A-HO-Z) C COMMON/DAT/ERHTDLXXCYCZCERCRMC, XOYOATPITPI2,PI,W, El, *E2,EERrAKOAKJ,AKKFAOFFSET(7),WDELTAOFFLIM,ERROR,THI,FI, *NS1,NS2,NSS2,NOFF C X2=X*X AKO2=AKO *"AK AK2=AK*AK RM=DSQRT (AK2-X2) RMO=DSQRT (X2-AKO2) RMH=1RM*H RMT=RM*T RMOH=RMO *H CSH=DCOS (RMH) CST=DCOS (RMT) SNH=DSIN (RMH) SXN=RM*CST SXD= (ER*RMO*CSH-RM*SNH) * (2..DO+RMOH) * (SNH/RMO-CSH/RM) HZXE=SXN/SXD RETURN END C........................ C 1) This subroutine evaluates the integrand of the Green's C function at different points (subroutine Grei). 2) It evaluates the space integrals comming from the application of moments' method (subroutine adonis) 3) Multiply these two valueswith appropriate weighting coefficients and it adds them upZXX2*SAJO(K) SUBROUTINE FUNCT (IADAUPALOWN) IMPLICIT REAL*8 (A-HO-Z) REAL*4 S1,S2 COMPLEX ZSZS1S2,CI C"OMMON/COMP/ZS (50), ZSlS2 (350),NSNS1S2 COMMON/MAT/PLIAITI,V(3),IY COMMON/PUT/SSJO (250,7),SAJO (250,7),YSIN, YCOS COMMON/ADON/DIST(250,7,1O),RCOE(2f,25f,7,1O),AXSERS(5), *SERA(5), DARG ( 10, 4),S ( 10, 2), WREAL, NSER, 7NMAX) COMMON/DAT/ERHTDLXXCYCZCERCRMC, XOYOATPITPI2,PIWE1, *E2,EERfAKO0AKIAKKJAKFAOFFSET(7),WDELTAOFFLIMERRORTHIFI, *NS1,NS2,NSS2fNOFF COMMON/DATT/COAL(20),POINT(20),CN(51),BM(51),POLTM(20), *POLTE(20),AM(41),DM(41),POLES (40),VXXM(20),VZXM(20),VZXE(20), *BPOINT(10),BCOAL(10),MPOINTNPOINTNKOMANTMNTENKOK, IFIRST COMMON/COEF/RX, XX, RZ, XZ, FRX, FRZ, F X, F1Z

Print file "SLOT. FTM" Page I10 C function at different points C.......................................................................... SUBROUTINE GREI (X, XFM, XFE, IAD, TM) IMPLICIT REAL*8 (A-HO-Z) C COMMON/DAT/ERH,T,DLXXCYCZCERCRMCXOYOA,TPI,TPI2,PI,W,E1, *E2,EER, AK0QAK, AKK,FAOFFSET(7),WDELTAOFFLIM,ERRORTHIFI, *NSlNS2,NSS2,NOFF C COMMON/COEF/RX, XX, RZ, XZ, FRX, FRZ, FiX, F1Z C X2=X*X AK2=AK*AK AK02=AKO*AK0 RM=DSQRT (DABS (AK2-X2)) RMO=DSQRT (DABS (X2-AK02)) RMH=RM*H RMT=RM*T RMHT=RM* (-H+T) C CSH=DCOS (RMH) SNH=DSIN (RMH) CST=DCOS (RMT) SNT=DSIN (RMT) CSHT=DCOS (RMHT) SNHT=DSIN (RMHT) RM2=RM*RM RM02=RMO*RM0 CSH2=CSH*CSH ERMO=ER*RMO ERM02=ERMO*ERMO EXX=DEXP (-X*T/FA)/FA EXZ=DEXP (-X* (2.DO*H-T)/FA)/FA IF (IAD.NE.7) GO TO 100 EX=DEXP (RMH) TANH= (EX-1.DO/EX) / (EX+1.DO/EX) CSHH= (EX+1.DO/EX) /2.DO EX=DEXP (RMT) CSHT=0.5D0* (EX+1.DO/EX) SNHT=0.5D0* (EX-1.DO/EX) TANT=SNHT/CSHT EX=DEXP (RMHT) CSHHT=0.5D0* (EX+1.DO/EX) SNHHT=0.5D0* (EX-1.DO/EX) TANHT=SNHHT/CSHHT 100 IF (IAD.NE.1) GO TO 1 DEN=RM2+ (ERM02-RM2) *CSH2 RNOM=-RM2*SNT+ (RM2-ERM02) *CSH*SNHT XNOM=ER*RM*RMO*CST Cl=X/ RM RX=C1 *RNOM/DEN XX=C1 *XNOM/DEN FRX=F1X*EXX DEN=DEN* (RM02+AK02* (ER-i.DO) *CSH2) RNOM=CST* (RM2+ER*RM02) *CSH*SNH XNOM-CST*RM*RMO* (1.DO-(1.DO+ER) *CSH2) CI=X*RM RZ=-C1 *RNOM/DEN XZ=Cl*XNOM/DEN FRZ=F1Z*EXZ RETURN

?rint file FVSLoT.IF2FPe 1Page I I. IF (IAD.NE.3) GO TO 2 C1=X-XFM IF (DABS(AK-X).LT.1.D-6) GO TO 10 DEN=ERMO *CSH RM *SNH RNOM= (RMP*CSHT-ERMO *SNHT) C2=X/RM RX=C1 *C2 *RNOM/DEN DEN=DEN* (RM*CSH+*RMO *SNH) RNOM=CST C3=X*RM RZ=C1 *C3 *RNOM/DEN FRX=F1X*EXX FRZ=F1Z *EXZ RETURN 10 RNOM=1. DO -ERMO* (-H+T) RX=C1 *X*RNOM/ERMO FRX=F1X*EXX RZ=X*C1/ (ERMO* (1.DO+RMO*H)) FRZF1lZ*EXZ RETURN 2 IF (IAD.NE.5) GO TO 4 C 1=X-XFE IF (DABS (AK-X).LT.1.D-6) GO TO 13 RNOM=RM*CSHT-ERMO *SNHT DENERMO *CSH-RM* SNH RX= (X/RM) *RNOM/DEN FRXF 1X*EXX RNOM=RM*CST DEN=DEN* (RM4*CSH+RMO*SNH) RZ=X*C1 *RNOM/DEN FRZ=F1Z*EXZ RETURN 13 RX=X* (1.DO-ERMO *(-H+T) ) /ERMO FRX=F1X*EXX RZ=X*C1/ (ERMO* (1.DO+RMO*H)) FRZ=F1Z*EXZ RETURN 4 IF (IAD.NE.7) GO TO 6 IF (DABS (X-AK).LT.1.D-6) GO TO 15 DEN=ERMO +RM*TANH RNOM= (RM- ERMO*TANHT) * (CSHT-TANH*SNHT) RX= (X/RM) *RNOM/DEN FRX=F1X*EXX RNOM=X* (RM*CSHT)/ (CSHH*CSHH) DEN=DEN* (RM+RMO*TANH) RZ=RNOM/ DEN FRZ=F1Z *EXZ RETURN 15 RX=X* (1. D 0 -ERMO * (-H+T) ) /ERMO FRX=F1X*EXX RZ=(X/ERMO) / (1.DO+RMO*H) FRZ=F1Z*EXZ 6 CONTINUE RETURN END ARIS SUBROUTINE ARIS

Print file "SLOT. TN"ae Page 12 IMPLICIT REAL*8 (A-HO-Z) COMMON/DAT/ERHTDLXXCYCZCERCRMCXrYOA,TPFITPI2,PIWE1, *E2,EERAKOAKAKK, FAOFFSET(7),WDELTAOFFLIMERRORTHIFFI, *NS1,NS2,NSS2,NOFF COMMON/DATT/COAL(20),POINT (20)rCN(51),BM(51),POLTM(20), *POLTE (20),AM(41), DM(41),POLES (40)fVXXM(20),VZXM(20),VZXE (20), *BPOINT(10),BCOAL(10),MPOINTNPOINTNK0,MANTM,NTENK0KIFIRST COMMON/ADON/D IST(250,7, 10), RCOE (20,250,7, 10),, AX, SERS (5), *SERA(5),DARG(110 4),S(10 2) WREAL, NSER, NMAX (7) COMMON/COEF/RX, XX, RZ, XZ, FRX, FRZ, FiX, F1Z + --- —------------------------------------------------— + Formation of the matrices: DIST, DARG, RCOE + --- —------------------------------------------------— + W2=W/2.DO U=WREAL/W THMIN=DATAN (DSQRT(1.DO/(U*U)-1l.DO)) THMAX=P I -THMIN AX- (THMAX-THMIN) /2.DO BX= (THMAX+THMIN) /2.DO X=PI/4.DO DO 1 J=1,NOFF MAX=NMAX (J) LPOINT=MPOINT IF (OFFSET(J).LE.OFFLIM) LPOINT=NPOINT DO 2 I=1,LPOINT POIN=BPOINT(I) IF (OFFSET(J).LE.OFFLIM) POIN=POINT (I) FI=X*(POIN+1.DO) THETA=AX*PO IN+BX AS=DSIN (FI) AC=DCOS (FI) DARG(I, 1)=W2*AC DARG(I,2)=AC DARG(I, 3)=AS DARG (I,4)=X DO 3 N=1,MAX AXN=FLOAT (N-2) *DLX IF (OFFSET(J).GT.OFFLIM) GO TO 4 DIST (N, J, I)=AXN*AS GO TO 5 4 AXN2=AXN*AXN BXN=OFFSET(J) -W*DCOS (THETA) /2.DO BXN2=BXN*BXN DIST (N, J, I) =DSQRT (AXN2+BXN2) SIG-DIST(N, J, I) SIG2=SIG*SIG SIG3=S IG2 *" IG DSIG=DABS (AXN) /SIG DSIG2=BXN2 /SIG3 DSIG3=-3.DO*DSIG*DSIG2/SIG DSIG4=-3.DO*DSIG2*(DSIG2-4.DO*DSIG**2/SIG)/SIG DSIG5=-3.DO*(-15.DO*DSIG2**2*DSIG+(20.DO/SIG)* DSIG2*DSIG**3)/SIG2 DSIG6=-3.DO*(-15.DO*DSIG2**3+(180.DO/SIG)*DSIG2 **2*DSIG**2.(120.DO/SIG2)*DSIG2*DSIG**4)/ SIG2 DSIG7=-3.DO*(525.D0*DSIG2**3*DSIG-(2100.DO/SIG)* DSIG2**2*DSIG**3+(840.DO/SIG2)*DSIG2*DSIG * **5) /SIG3

Print file "StOT.F2J" DSIG8=-3.DO*(525.DO*DSIG2**4-(12600.DO/SIG)*DSIG2 **3*DSIG**2+(25200.DO/SIG2)DSG2 *DSIG2**2*DSIG**4 -(6720.DO/SIG3) *DSIG2*DSIG**6) /SIG3 Evaluation of the coefficients Gij G21=DSIG2 G22=DSIG**2 G4 1=DS 1G4 G42=4.DO*DSIG3*DSIG+3.DO*DSIG2**2 G43=6.DO*DSIG2*DSIG**2 G4 4=DS IG* * 4 G6 1=DS 1G6 G62=6.DO*DSIG5*DSIG+15.DO*DSIG4*DSIG2+10.DO*DSIG3**2 G63=15.DO*DSIG4*DSIG**2+6.DO*DSIG3*DSIG2*DSIG+15.DO **DSIG2**3 G64=20.DO*DSIG3*DSIG**3+45.DO*DSIG2**2*DSIG**2 G65=15.DO*DSIG2*DSIG**4 G6 6=DS IG* * 6 G8 1=DS 1G8 G82=8.DO*DSIG7*DSIG+28.DO*DSIG6*DSIG2+56.DO*DSTG5 *DSIG3+35.DO*DSIG4**2 G83=28.DO*DSIG6*DSIG**2+168.DO*DSIG5*DSIG2*DSIG+ 280.DO*DSIG4*DSIG3*DSIG+210.DO*DSIG4*DSIG2**2+ 28Q.DO*DSIG3**2*DSIG2 G84=56.DO*DSIG5*DSIG**3+420.DO*DSIG4*DSIG2*DSIG**2 +280.DO*DSIG3**2*DSIG**2+840.DO*DSIG3*DSIG2**2 *DSIG+105.DO*DSIG2**4 G85=70.DO*DSIG4*DSIG**4+560.DO*DSIG3*DSIG2*DSIG**3 +420.DO*DSIG2**3*DSIG**2 G86=56.DO*DSIG3*DSIG**5+210.DODSSIG2**2*DSIG**4 G87=28.DO*DSIG2*DSIG**6 G8 8=DS IG* * 8 RCOE (2, N, J, I) =-0. 5D0* (G22+SIG*G21) RCOE(1,NJ, I) =0. 5D0* (G22-SIG*G21) ----— ~ - - -- -- --- - - 1 ~ SX0.5D0*SIG*(G42-SIG*G41) S30=-O.5DO*SIG*(G42+SIG*G41) S31=0.25DO*(SX+3.DO*G43) S33=0.25D0*(SX-G43) RCOE (3,N, J, I) 0. 5D0 *(SIG*S33/3.DO+G44/4.DO) RCOE(4,N,J,fI)=0.5DO*(SIG*S31+SIG*S33/3.DOG44) RCOE(5,NJI)=0.5D0*(SIG*S31+3.DO*G44/4.DO) RCOE (6, N, J, I) =SIG*S30 SX=SIG*S33/3.DO+G64/4.DO ST=SIG*S31+SIG*S33/3.DO-G64 S5M3=SIG2 *S30 S5M1=O.5DO*SIG*(SIG*S31+3.DO*G64/4.DO) S51=0.25D0*(0.5DO*SIG*ST5.DO*G65/2.DO) S53=0.25DD0*(0.5DOD*SIG*ST+0.25D0*SIG*SX+0.5DG65/ 4.DO) S55=0.125D0*(0.5DO*SIG*SX-0.5*G65) RCOE (7, N, J, I) =0. 5D0 *(SIG*S55/5.DO+G66/16.DO) RCOE (8,N, J, I) =0. 5DO* (SIG*S53/3.DO+SIG*S5515.DO6.DO*G66/16.DO) RCOE (9,N, J, I) =0. 5DO* (SIG*S51+SIG*S53/3.DO+15.DO* G66/16.DO) RCOE(10, N, J, I)=0.5DO*(SIG*S5-101.DO*G66/16.DO) RCOE (1 1, N, J, I) -SIG* S 5M2 RCOE (12,r N, J, I) =SIG*S5M3 Page 13

Print file "SL~OT.FTNp Page 14 S7M5=SIG2*S5M3 S7M3=SIG2 *S5M1 S7M1=O.5DO0*SIG*(SIG*S51-10.DO*G86/16.DO) S71=0. 5D0 (0.25DO*SIG*(SIG*S51+SIG*S53/3.DO+ 15.D0OG86/16.DO)+35.DO*G87/32.DO) S73=0.5D(0.25DO(2 *SIG*(SIG*S51+SIG*S53/3.DO+15.DO *G86/16.DO)+0.125DO*SIG*(SIG*S53/3.DO+SIG* S55/5.DO-6.DO*G86/16.DO)-21.DO*G87/32.DO) S75=0. 5D0 (0.125D0*SIG* (SIG*S53/3.D+SIG*S55/5.D6.DO*G86/16.DO) +(SIG/12.DO)*(SIG*S55/5.DO+ G86/16.DO)+7.DO*G87/32.DO) S77=0.5DO((SIG/12.DO)*(SIG*S55/5.DO+G86/16.DO)G87/32.DO) RCOE(13,NJI)=0.5D0*(SIG*S77/7.DO+G88/64.DO) RCOE(14,N,JI)=0.5D0*(SIG*S75/5.DO+S77*SIG/7.DO -8.DO*G88/64.DO) RCOE(15,NJI)=0.5D0*(SIG*S73/3.DO+SIG*S75/5.DO +28.DO*G88/64.DO) RCOE(l6,NJIh=0.5D0*(SIG*S71+SIG*S73/3.DO-56.D *G88/64.DO) RCOE(17,N, J, I)=0.5DO* (SIG*S71+70.DO*G88/64.DO) RCOE (18,N, J, I) =SIG*S7M1 RCOE(19,N,J, I)=SIG*S7M3 RCOE (20, N, J, I) =SIG*S7M5 5 CONTINUE 3 CONTINUE 2 CONTINUE 1 CONTINUE Formation of the series s(dlx). Storage in vectors SERS(5),SERA(5) U1=2.DO*THMIN/FLOAT (NSER) DO 6 JN=1,NSER S2= (2.DO*FLOAT (JN) -1.DO) S2=S2/ (2.DO*FLOAT (NSER)) S3=DCOS (S2 *THMIN) S (JN, 2) S3*w/2. DO S (JN, 1)=U1 6 CONTINUE ADL=AKK*DLX ADL2=ADL*ADL ADL3=ADL2 *ADL ADL4=ADL3 *AD ADL5=ADL4 *AD ADL 6=ADL5 *AD YS IN=D S IN (AD L) YCOS=DCOS (ADL) SER1=(1.DO-YCOS) *2.DO/AKK SER2=-YSIN/3.D0+ADL*YCOS/4.D0+ADL2*YSIN/l0.DO ADL3*YCOS/36.DO -ADL4*YSIN/168.D0+ADL5*YCOS/960.D0+ADL6*YSIN/6480.DO SER3=YSIN/60.D0 ADL*5.D0*YCOS/360.D0 ADL2*YSIN/168.D0+ADL3 * *YCOS/560.0D+ADL4*YSIN/2592D0 -ADLL5*YCOS/12960.D0-ADL6 *ySIN/95040OD0 SER4=-YSIN/2520.DO+ADL*YCOS/2880.DO+ADL2 *YS IN/6480.DO-ADL3 * *YCOS/216000.D-ADL4*YSIN/95040.D0+ADL5*YCOS/518400.DO SER5=YSIN/18 1440.DO-ADL*YCOS/2O1600.DO-ADL2*YSIN/443520.DO+ ADL3*YCOS/1442775.9DO SERS (1)=SER1*SER1 SERS (2) =DLX*2.DO*SER1*SER2

Print file "~SLOT FTN ag 1 Page 15 SERS (3)-DLX* (DLX*SER2*SER2+2.DO*SER1*SER3) SERS (4)-=DLX* (2.DO*SER1*SER4+2.DO*DLX*SER2*SER3) SERS (5) =DLX* (DLX*SER3*SER3+2.D0*DLX*SER2*SER4) C SERA(1) =SER1 SERA (2) =DLX*SER2 SERA (3) =DLX*SER3 SERA(4) =DLX*SER4 SERA (5) =DLX*SER5 111 CONTINUE RETURN END ADONIS This subroutine evaluates the space integaris of the bessel function SUBROUTINE ADONIS IMPLICIT REAL*8 (A-H, O-Z) DIMENSION BJ(1O,2),DERIV(9,3) COMMON/ADON/DIST(250,7,10),RCOE(20,250,7,10),AXSERS(5),. *SERA(5), DARG ( 10, 4), S ( 10, 2) f WREALf NSER, NMAX (7) COMMON/PUT/SSJO(250,7),SAJO(250,7),YSINYCOS COMMON/DAT/ERHT,,DLX,,XC,,YC,,ZC,,ERC,,RMC,,XO,,YO,A,,TPI,fTPi2,,PI,,W,,El, *E2.EER, AKOQf Aj(,AKK, FAf OFFSET (7) fWDELTA, OFFLIM, ERROR, THIrFI, *NS1,NS2,NSS2,NOFF COMMON/DATT/COAL (20),POINT (20),CN(51),BM(51),POLTM(20), *POLTE (20),AM(41),DM(4l),POLES(40)rVXXM(20)fVZXM(20),VZXE (20), *BPOINT(10),BCOAL(10),MPOINTNPOINTNKOMANTMNTENKOKIFIRST C:OMMON/BSS/ARG(10),AARG COMMON/MAT/PLIIAIrTIrV (3), IY COMMON/COEF/RX, XX, RZ, XZr FRX, FRZr FlX, F1Z ARX=W*AX/2.DO Wl=2.DO*YCOS PR1=PLI *DLX PR2=PR1 *PR1 PR4=PR2 *PR2 PR6=PR4 *PR2 PR8=PR6 *PR2 DO 1 J=1,NOFF MAX=NMAX (J) DO 2 N=1,MAX SSJO (N,J)=0.DO SAJO (N, J) =0.DO 2 CONTINUE 1 CONTINUE DO 11 J=1,NOFF LPOINT=MPOINT IF (OFFSET(J).GT.OFFLIM) GO TO 12 LPOINT=NPOINT DO 13 I=1,NPOINT ARG (I) =P LI *DARG (I,, 1) 13 CONTINUE CALL BESS1(BJ) 12 DO 14 I=1,LPOINT DO 17 NK=1,5

Print file "SL~OT.F7l" DERIV(NK, 1)=O.DO DERIV(NK,2)=O.DO 17 CONTINUE AS INARX*BCOAL (I) IF (OFFSET(J).GT.OFFLIM) GO TO 15 AS IN=W*DARG (I, 4) *COAL (I) AROF=PLI*OFFSET (J) *DARG (I,2) COFF=DCOS (AROF) SSUM=O.DO DO 16 JN=1,NSER ARAF=PLI* S (JN, 2) *DARG (I, 2) CAFF=DCOS (ARAF) SSUM=SSUM+S (JN, 1) *CAFF 16 CONTINUE 15 CONTINUE KMAX=NMAX (J) DO 18 K=1,KMAX DO 20 NK=1,5 DERIV (NK, 1)=DERIV (NK, 2) DERIV (NK, 2) =DERIV (NK, 3) 20 CONTINUE IF (OFFSET(J).GT.OFFLIM) GO TO 21 SIN1-DARG(I, 3) SIN2=SIN1*SIN1 COS1=DCOS (PLI*DIST (K, J, I)) TERM=COFF* (BJ(I, 1) -SSUM/PI) *COS 1 DERIV(1, 3) =TERM S IN1=S IN2 DERIV(2, 3)=TERM*SIN1 SIN1=SIN1*SIN2 DERIV(3, 3) =TERM*SIN1 SIN1=SIN1*SIN2 DERIV(4, 3)=TERM*SIN1 SIN1=SIN1*SIN2 DERIV(5,3)=TERM*SIN1 GO TO 22 21 AARG=PLI *DIST (K, Jr I) ARG2=AARG*AARG ARG4=ARG2 *ARG2 ARG6=ARG4 *ARG2 CALL BESS2(BJ) DERIV(1, 3) =BJ(1, 2) DERIV(2, 3) =RCOE (1, K, J, I) *BJ(3,2) + RCOE(2, K, J,rI) *BJ(1,2) DERIV(3, 3) =RCOE (3, K,J, I) *BJ(5, 2) + RCOE (4,fK, J, I) *BJ(3,2) +(RCOE (5,K,J, I) +RCOE(6, K, J, I) /ARG2) *BJ(1,2) DERIV(4, 3)=RCOE(7,K,KJ, I) *BJ(7,2)+ RCOE(8, K, JrI)*BJ(5,2)+RCOE(9,K, JrI)I) BJ(3,2) +(RCOE (10,fK,rJ, I) +RCOE (11,rK, J, I) /ARG2+RCOE (12,rK,rJ, I) /ARG4) * BJ(1,2) DERIV(5,3)=RCOE(13, K, J, I)*BJ(9,2)+ RCOE(14,fK, J, I)*BJ(7,2)+RCOE(15, K, J, I)*BJ(5,2)+RCOE(16,fK, J, I)*BJ(3,2)+ (RCOE(17,fK,rJ, I)+RCOE(18,fK, J, I) /ARG2 +RCOE (19, K,rJ, I) /ARG4+RCOE (20,fK, J, I) /ARG6) *BJ(1,2) 22 IF (K.LT.3) GO TO 18 SUMS=SERS (1) *DERIV(1, 2) -PR2*SERS (2) *DERIV(2,2) * +PR4*SERS (3) *DERIV(3, 2) -PR6*SERS (4) *DERIV (4,2) +PR8*SERS (5) *DERIV(5,2) CH1=SERA(1)* (DERIV(1, 1)+DERIV(1,3) W1*DERIV (1,2)) CH2-SERA(2)* (DERIV(2, 1)+DERIV(2,3) -W1*DERIV Page 16

Print file "SLOT. F"Pa 1 Page 1 7 (2,2) )*PR2 CH3=SERA(3)* (DERIV(3, 1)+DERIV(3, 3) W1*DERIV (3,2)) *PR4 CH4=SERA(4) * (DERIV(4,.1) +DERIV(4, 3) W1*DERIV (4,2)) *PR6 CH5=SERA (5) * (DERIV (5, 1) +DERIV (5, 3) -W1 *DERIV (5,2)) *PR8 SUMA=CH1 -CH2 +CH3 -CH 4 +CH5 KJ=K-2 SSJO (KJ, J) =SSJ0 (KJ, J) +ASIN*SUMS SAJO (KJ, J) =SAJ0 (KJ, J) +AS IN* SUMA CCCC C IF (KJ.EQ.1)WRITE (6,665) KJJ,SSJ0(KJJ), C * SUMSSAJO(KJJ),SUMA C665 FORMAT(10X,'KJ=',14,2X,' J=',I4/10X,'SSJO=', C * E14.7,2X,'SUMS=',E14.7/10X,'SAJ='f,E14.7, C * 2X,'SUMA=',E14.7/) CCCC 18 CONTINUE 14 CONTINUE 11 CONTINUE RETURN END BESS1 This subroutine gives values for the zeroth order Bessel functions. It is used for small offsets SUBROUTINE BESS1(BJ) IMPLICIT REAL*8 (A-HO-Z) DIMENSION BJ(1 0,2) COMMON/COEF / RX, XX, RZ, XZ, FRX, FRZ, FiX, F 1 Z COMMON/ADON/DIST (250,7, 10), RCOE (20,250, 7, 10), AX, SERS (5), *SERA (5),DARG (10, 4),S (10, 2) WREAL, NSER, NMAX(7) COMMON/BSS/ARG(10),AARG COMMON/DATT/COAL (20),POINT (20),CN(51),BM(51),POLTM(20), *POLTE (20),AM(41),DM(4l),POLES(40), VXXM(20),VZXM(20),VZXE(20), *BPOINT(10),BCOAL(10),MPOINTNPOINTNKOMANTMNTENKOKIFIRST PI=3. 141592653589D0 DO 1 IJ=1,NPOINT X=ARG(IJ) IF (X.GT.0.001D0) GO TO 10 X3=X/3.DO X32=X3 *x3 X34=X32*X32 X3 6=X3 4 *X32 BJ0=1.DO-2.2499997D0*X32+1.2656208D0*X34-0.3163866D0 **X36 BJ(IJ, 3)=BJO GO TO 1 10 IF (X.GT.3.DO) GO TO 12 X3=X/3.DO X32=X3 *x3 X3 4=X32 *X32 X36=X34*X32 X3 8=X3 6 *X32 X310=X38*X32 X312= X310 *X32 BJO=1.DO-2.2499997D0*X32+1.2656208D0*X34-0.3163866D0 *X36+0.0444479D0*X38-.00039444D0*X310+000021000 DO*X312

Print file "SLOT.FTN" Page 18 BJ(IJ,1)=BJO GO TO 1 12 CONTINUE X3=3.DO/X x32=x3*X3 X33=X32 *x3 X34=X33 *x3 X35=X34*X3 X36=X35* x3 FJ0=0.79788456D-0.00000077D0*X3-0 00552740D0*X32-0.0000 9512D0*X33+0.'00137237D0*X34-0.00072805D0*X35+0.00014 476D0*X36 TJO=X-0. 78539816D0-0. 04166397D0*X3-0. 00003954D0*X32+0.00 262573D0*X33-0.00054125D0*X34-0.00029333D0*X35+0.000 13558D0*X36 WCON=DSQRT(1.DO/X) BJ (IJ, 1) =WCON*FJO*DCOS (TJO) 1 CONTINUE RETURN END TAIL This subroutine evaluates the tail contribution SUBROUTINE TAIL IMPLICIT REAL*8 (A-HO-Z) COMPLEX ZS,ZS1S2 DIMENSION MAX(8,2) COMMON/CTAIL/S1 (4,205,7),D1 (4,205,7),D2 (4,205,7), *T1 (4,205,7), T2 ( 4,205,7) ) T3 (4,205,7) f T4 (4,205,7) COMMON/COMP/ZS (50), ZS1S2 (350),NSN51S2 COMMON/DAT/ERHTDLXXCYCZCERCRMCXOYOATPITPI2,PIW,El, E2, EERAK fAKfAECfFA, OFFSET (7),WDELTAOFFLIM,ERRORTHIFI, *NS1,NS2,NSS2,NOFF COMMON/INT/XNS(40),CNS(40),XND(20,2),CND(20),XNT(40,3), *CNT(40),NDPNTPNSP COMMON/ADON/DIST(250,7,10),RCOE(20,250, 7,10),AXSERS(5), *SERA (5), DARG (10f 4), S (10f 2) f WREALrNSER4NMAX (7) COMMON/OUT/GS(50),GS1S2(350) COMMON! IOFF/INS, 1NS1S2 This vector contains the values of t in the integrals hO Z1=T Z2=2.D0*H-T This vector contains the values of the coefficient C in the integrals hO C1=FA This vector contains the number of elements of the matrices ZS,Z512,,.. MAX(1, 1)=NS MAX(2,1)=NS1S2 MAX(1, 2) =INS MAX(2,2) =INS1S2

'rirt file "SLOT.JFTN" Page 19 This vector contains the values of the coefficient A in the integrals hO AKK=TPI W2=W/2.DO THMIN=WREAL/W THMIN=DATAN(DSQRT(l.DO/THMIN**21.DO)) THMAX=P I-THMIN PI2=PI/2.DO PI4=PI/4.DO DLX2=DLX/2.DO DLX4=DLX2 *DLX2 YCOS=DCOS (AKK*DLX) CCS=DCOS (2.DO*AKK*DLX) YSIN=DSIN (AKK*DLX) SSN=DSIN(2.DO*AKK*DLX) + --- —------------------------------— + I Evaluation of S1,S2,S3,S4,S5,S6 I (Single Integrals) + --- —------------------------------— + DO 201 J=1,7 DO 202 K=l,205 DO 203 JK=l,4 Sl(JK, K, J) =0.D0 Dl(JKK,J)=O.D0 D2 (JK, K, J) =0.DO Ti(JKK,J)=O.DO T2(JKK,J)=O.DO T3(JK,K,J)=0.DO T4(JKK,J)=O.DO 203 CONTINUE 202 CONT INUE 201 CONTINUE ZPl=Zi/C1 ZP2=Z2/C1 ZP12=ZP1*ZP1 ZP22=ZP2*ZP2 DO 1 J=1,NOFF KMAX=NMAX(J) +2 IF (OFFSET(J).LT.i.D-6) THMAX=PI DSP= (THMAX-THMIN)/4.DO DDP=DSP *DLX2 DTP=DSP *DLX4 COEFi= (THMAX-THMIN) /2.DO IF (OFFSET(J).LT.1.D-6) COEF1=(PI/2.DO-THMIN)/2.DO COEF2=(THMAX+THMIN) /2.DO IF (OFFSET(J).LT.l.D-6) COEF2=(PI/2.DO+THMIN)/2.DO DO 10 I=iNSP THI=COEF1*XNS (I) +COEF2 Ci=DCOS(THI) C2=W2 *C1 C2=OFFSET(J)-C2 CW=C2 *C2 AASIN=CNS (I) *DSP DO 11 K=i,KMAX XN=FLOAT (K-3) *DLX RAD2=XN*XN+CW TRADl=DSQRT (RAD2+ZPi2)

Prin t file. "SLOT. FTN " ge2 Page 2 0 TRAD2=DSQRT (RAD2+ZP22) S1 (lrKrJ)=S1 (1, KrJ) +DLOG (2.DO* (TRAD1+XN) )*AASIN Si (2, K, J) =S1 (2, K, J) +DLOG (2.DO* (TRAD2+XN) ) *AASIN 11 CONTINUE 10 CONTINUE ----------------------------------------------------------------------— + I EVALUATION OF Dl,D2,D4,D5 + --- —---------------------------------------------------------------— + DO 20 I=1,NDP THI=COEF1*XND (I, 1) +COEF2 XI=DLX2* (XND (I,2) +1.DO) C1=DCOS (THI) C2=W2 *Cl C2=OFFSET (J) -C2 CW=C2*C2 AASIN=CND (I) *DDP SV1=DSIN(AKK* (DLX-XI)) SV2=-SV1 SV4=DSIN (AKK*XI) C2=DCOS (AKK* (DLX-XI)) DO 21 K=1,KMAX XPXI+FLOAT (K-2) *DLX XNM=-XI+FLOAT (K-2) *DLX RADP2=XNP *XNP+CW~ RAIDM2-XNM*XNM+CW TRAP1=DSQRT (RADP2+ZP12) TRAP2=DSQRT (RADP2+ZP22) TRAM1=DSQRT (RADM2+ZP12) TRAM2=DSQRT (RADM2+ZP22) XA1=AKK*XNP XA2=AKK*XNM XAP=DSIN (XA1) XAM=DS IN (XA2)/ SANP1=XAP*DLOG (2.DO* (TRAP1+XNP)) SANP 2=XAP *DLOG (2..DO* (TRAP2 +XNP) ) SANM1=XAM*DLOG (2. DO* (TRAM1+XNM) ) SANM2=XAM*DLOG (2.DO* (TRAM2+XNM)) XAP=DSIN(XAl/2.DO) XAM=D SIN (XA2 /2..DO) S ONP 1=XAP /TRAP 1 SONP2=XAP /TRAP2 SONM1=XAM/ TRAMi SONM2=XAM/ TRAM2 'Y1=-XNM/2.DO-DLX,Y2=-xNP/2.DO+DLX CY1=DCOS (AKK*Yl) CY2=DCOS (AKK*Y2) SY1=DSIN (AKK*Yl) S Y2=D SIN (AKK *Y2) Dl1 (1, K, J) =D 1 (l, KJ) + (SANP 1+SANM1) *SV2 *AS IN D2 (1, KJ) =D2 (1,K, J) +(CY1*SONPl-CY2*SONM1) *ASIN Dl1 (2,r Kr J) =D 1 (2,f Kf J) + (SANP 2+SANM2) * SV2 * AS IN D2 (2, KJ)=D2 (2, KJ) +(CY1*SONP2-CY2*SONM2) *ASIN 21 CONTINUE 20 CONTINUE evaluation of T1,T2,T3,T4

Print file "SLOT.FTN"P Page 21 C DO 30 I=1,NTP THI=COEFi*XNT (I, 1) +COEF2 XI=DLX2* (XNT (I,2) +1. DO) XIP=DLX2* (XNT (I, 3) +1.DO) C1=DCOS(THI) C2=W2 *Cl C2=OFFSET (J) -C2 CW=C2*C2 SVi=DSIN(AKK* (DLX-XI)) SV2=-SV1 SV3=DSIN(AKK* (DLX-XIP)) AASIN=DTP*CNT (I) DO 31 K=1,KMAX XNPP= (XI+XIP) +FLOAT (K-i) *DLX XNPM= (XI-XIP) +FLOAT (K-i) *DLX XNMP= (-XI+XIP) +FLOAT (K-i) *DLX XNMM= (-XI-XIP) +FLOAT (K-i) *DLX RADPP2=XNPP *tXNPP +CW RADPM2=XNPM*XNPM+CW RADMP 2 =XNMP * YNMP +CW RADMM2=XNMM*XNMM+CW TAPPi=DSQRT (RADPP2+ZP12) TAPP2=DSQRT (RADPP2+ZP22) TAPMi=DSQRT (RADPM2+ZPi2) TAPM2=DSQRT (RADPM2+ZP22) TAMPi=DSQRT (RADMP2+ZP12) TAMP2=DSQRT (RADMP2+ZP22) TAMMi=DSQRT (RADMM2+ZP12) TAMM2=DSQRT (RADMM2+ZP22) CST1=DCOS (AKK* (XNPM/2.DO+DLX) ) *DSIN((AKK*XNPP /2.DO) CST2=DCOS (AKK* (-XNMP/2.DO+DLX) ) *DSIN (AKK*XNMM /2.DO) CST3=DCOS (AKK* (XNMMI2.DO+DLX) ) *DSIN (AKK*XNMP /2.DO) CST4=DCOS (AKK* (-XNPP/2.DO+DLX) ) *DSIN (AKK*XNPM /2.DO) Ti (1, K, J) =Ti (1, K, J) +SV2*AASIN*CSTi/TAPP1 T2 (1,fKrJ) =T2 (1, K, J) +SV1 *AAS IN*CST2 /TAMM1 T3 (1, K, J) =T3 (1 f K, J) +SV1 *AAS IN*CST3 / TAMP 1 T4 (irK, J) =T4 (ifK, J) +SV2*AASIN*CST4/TAPMi Ti (2, K, J) =Ti (2, K, J) +SV2*AASIN*CST1/TAPP2 T2 (2,KJ)=T2 (2,KJ)+SVi*AASIN*CST2/TAMM2 T3(2,KJ)=T3(2,KJ)+SV1*AASIN*CST3/TAMP2 T4 (2, K, J) =T4 (2, K, J) +SV2 *AAS IN*CST4 /TAPM2 31 CONTINUE 30 CONTINUE 1 CONTINUE Evaluation of GSGSiS2 CZX=(1.DO-ER) / ((i.DO+ER)*(1.DO+E2)*(l.DO+0.5D0*El)) CXX=3.DO CS=TPI2*CZX/FA CAX=TP I *JCXX/ FA CAZ=TP I *CZX/FA DO 4 JM=lNOFF NJMAX=MAX (JM, 1) J=MAX(JM, 2) DO 62 N=I,NJMAX NP i=N+2 NO=N+i NM1=N

Prin t filI "SLOT.FM JF S T= — D1 (2,fNP 1,rJ) + 2. DO* YCO S *D1 (2,fNO, rJ) - Di (2,rNM 1,fJ) *+2.DO* (T1 (2, NJ)+T2 (2, NJ)-T3 (2f NJ)-T4 (2, NJ) MP2-N+4 MP i=N+3 MO=N+2 MMi=N+l1 t4M2=N SINP2=DSIN (AKK*FLOAT (N+i) *DLX) SINP1=DSIN (AKK*FLOAT (N) *DLX) SINO=DSIN (AKK*FLOAT (N-i) *DLX) SINMi=DSIN(AKK*FLOAT (N-2) *DLX) SINM2=DSIN (AKK*FLOAT (N-3) *DLX) ATX=SINP2*Si (1f MP2, J) -4.DO*YCOS*SINPi*Si (1fMPiJ) * ~+2.DO* (2.DO+CCS) *SINO*Si (i, MO, J) -4.DO*YCOS * ~~*SINM1*Sil(iMM1 J) +SINM2*S1l(iMM2, J) ATZ=SINP2*Sl(2,MP2,J)-4.DO*YCOS*SINPi*Sl(2,MPiJ) * ~+2.DO* (2.DO+CCS) *SINO*Sl (2, MO, J) -4.DO*YCOS * ~~*SINM1*Sil(2,MMi J) +SINM2*Sil(2,MM2, J) AAX=-2.DO* (D2 (i, NPi, J) 2.DO*YCOS*D2 (i, NO,J) * ~+D2 (iNMlJ)) AAZ=-2.DO* (D2 (2,NPi,, J) -2.DO*YCOS*D2 (2,,NO,, J) * ~+D2 (2,NM1,J)) AX=ATX+AAX AZ=ATZ +AAZ ZW=W* (CS*ST+CAX*AX-CAZ*AZ) IF (JM. EQ. 1) GS (N) =ZW IF (JM.EQ.2) GSiS2(N)=ZW 62 CONTINUE 4 CONTINUE RETURN END This subroutine evaluates the higher order bessel ucin sn the ascenting series expression or hankel's expansion. SUBROUTINE BESS2 (BJ) IMPLICIT REAL*8 (A-HO-Z) DIMENSION BJ(lO,2),U(4),RBJ(50,2) C"OMMON/Ba i/BJO, BJi CO:MMON/BSS/ARG(1O),X PI=3.141592653589 Evaluation of JO,Ji CALL BSJO(X) RBJ (1,2) =BJO RBJ(2,2) =BJ1 NC0N=1 N=IDINT (2. 4D0*X) IF (N.LT.iO) N=iO IF (X.LT.3.DO) GO TO 10 EVALUATION OF HIGHER ORDER BESSEL FUNCTIONS UP TO ORDER LESS THEN THE ARGUMENT NIMAX=IDINT (X) -1 IF (NIMAX.GT.9) NIMAX=9 DO 1 I=2,NIMAX NJi=I NJ2=I-i NB=I+l RBJ (NB, 2) =FLOAT (2 *NJ2) *RBJ (NJ1i 2) /X-RBJ (NJ2,,2) 1 CONTINUE IF (NTMAXEnQ19 GOn TO) 20 Page 22

?rint file "SLOT N FWae j.Page 23 NCON=NIMAX DEBYE'S ASYMPTOTIC EXPANSION-EVALUATION OF JN 10 DO 11 J=1,2 JN=N-J+ 1 XA=X/FLOAT (JN) XA=1.DO/XA XE=XA+DSQRT (XA*XA-1.DO) A=DLOG(XE) CTH=(XE+1.DO/XE)/(XE-i.DO/XE) CALL F(CTHU) TNH=1.DO/CTH R1=DEXP (FLOAT (JN) * (TNH-A)) R2=DSQRT (2.DO*PI*FLOAT (JN) *TNH) BN1=JN BN2=JN*JN BN3=BN2 *JN BN4=BN3 *JN RBJ(JN+1,2)=(R1/R2) * (1.DO+U(1) /BN1+U(2) /BN2+U(3) /BN3+ * U(4) /BN4) 11 CONTINUE EVALUATION OF HIGHER ORDER BESSEL FUNCTIONS WHEN X<10 NJMAX=N-2 -NCON DO 2 I=1,NJMAX NJB=N-I NJB1 NJB+1 NJB2=NJB1+1 RBJ(NJB,2)=2.DO*FLOAT(NJB) *RBJ(NJB1,2) /X-RBJ(NJB2,2) 2 CONTINUE 20 CONTINUE DO 3 1=1,9 BJ(I,2)=RBJ(I,2) 3 CONTINUE RETURN END SUBROUTINE BSJO (X) IMPLICIT REAL*8 (A-HO-Z) COMMON/BO1/BJO, BJ1 Evaluation of JO using the series expansion given in Abramowitz. P1=3. 141592653589D0 IF (X.GT.3.DO) GO TO 20 X3=X/3.DO X32=X3 *X3 X3 4=X32 *X32 X3 6=X32 *X34 X38-X32*X36 X310=X38*X32 X312=X310*X32 BJO=1.DO-2.2499997D0*X32+1.2656208DO*X34-0.3163866D0*X36+ 0 O.444479DO*X38-0.0039444DO*X310+0.0002100ODO*X312 BJ1=X* (0.5D0.56249985D0*X32+0.21O93573DO*X34O. 03954289D0 *X36+0.00443319D0*X38-..00031761D0*X310+000001109D0 * *X312) GO TO 21 20 X3=3.DO/X X32-X3*X3

Print file "SLOT.rFTN" Page 24 X33=X32 *X3 X34-X33*X3 X35=X34*X3 X36=X35*X3 FJO=0.79788456D0-0.00000077D0*X3-0.00552740D0*X32-0.0000 9512D *x33+0.00137237D0*X34-0.00072805D0*X35+000014476D0*X36 FJ1=0. 79788456D0+0. 00000156D0*X3+0. 01659667D0*X32+0. 00017105D0 * *X33-000249511D0*X34+0.00113653D0*X35-.000020033D0*X36 TJO=X-0.785398 16D0-0. 04166397D0*X3-0. 00003954D0*X32+0.002 62573D0 *x33-0.00054125D0*X34-.00029333D0*X35+000013558D0*X36 TJ1=X-2.35619449D0+0. 12499612D*X3+0.00005650D0*X320.00637879D *X33+0.00074348D0*X34+0.00079824D0*X35-0.00029166D0*X36 WCONtDSQRT (1.DO/X) BJO=WCON*FJO*DCOS(TJO) BJ1=WCON*FJ1*DCOS (TJ1) 21 CONTINUE RETURN END SUBROUTINE F(XU) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION U(4) X2=X*X X3=X2*X X4=X3*X X5=X4*X X6=X5 *X X7=X6*X X8=X7 *X X9=X8 *X< X10 "x9*X X11=X10'~tX X12=Xll1*X U(1) =(3.D0*X-5.D0*X3) /24.DO U (2) = (81.DO*X2-462.DO*X4+385.DO*X6) /1152.DO UJ(3)=(30375.DO*X3-369603.DO*X5+765765.DO*X7-425425.DO*X9)/ 414720.DO U (4)=(4465125.D0*X4-94121676.DO*X6+349922430.DO*X8-446185740.DO* x10+185910725.DO*Xl2)/39813120.DO RETURN END SUBROUTINE DATASLOT This subroutine gives all the data for integration used in subroutine SLOT.FTN SUBROUTINE DATASLOT IMPLICIT REAL*8 (A-H, O-Z) COMMON/DAT/ER, H, T, DLX, XC, YC, ZC, ERC, RMC, XOYOA, TPIfTPI2, PIW, El, *E2,EERfAKOAKfAJK, FAOFFSET(7),WDELTAOFFLIMERRORTHIFI, *NS1,NS2,NSS2,NOFF COMMON/DATT/COAL(20),POINT(20),CN(51),BM(51),POLTM(20), *POLTE (20)fAM(41),DM(41),POLES(40),VXXM(20)fVZXM(20),VZXE(20), *BPOINT(10),BCOAL(10),MPOINTNPOINTNKOfMANTMNTENKOKf IFIRST COMMON/INT/XNS(40),CNS(40),XND(20,2),CND(20),XNT(40,3), *CNT(40),NDPNTPNSP COMMON/ADON/DIST(250,7,10),RCOE(20,250,7,10),AXSERS(5), *SERA(5),DARG(10, 4) fS (10,2) fWREALf NSERNMAX(7) COMMON/ IOFF/INS, INS1S2

?rint file "SLOT.FTN" Page 25 PI=3.141592653589D0 TPI=2.DO*PI TPI2=TPI*TPI ------------------— + I ERROR FUNCTIONS I + --- —-------------— + A1=A*A-TPI2 A2=TPI2-TPI2/ER El=0.5D0*A2/A1 E2=ER*E1/(1.DO+ER) FA=DSQRT (1.D+TPI2/A1) AK0=2.DO*PI/DSQRT(EER) AKK=2.DO*PI AK=AKO*DSQRT (ER) ----------------------------------------------------— +: I Data for the poles I IFIRST- 0: dominant mode is TM wave (many poles) I I 1: dominant mode is TE wave (many poles) I I 2: only one TM surface wave + --- —---------------------------------------------—....+ + --- —----------------------- IData for the Integration I + --- —----------------------— + NKO=20 NKOK=1 MA=20 NSER=10 NPOINT=10 Vector COAL COAL(1)=0.0666713443D0 COAL(2)=0.14945134915D0 COAL (3) =0.21908636251D0 COAL(4)=0.26926671931D0 COAL(5) =0. 29552422471D0 COAL(6)=COAL(5) COAL(7)=COAL(4) COAL(8)=COAL(3) COAL(9)=COAL(2) COAL(10)=COAL(1) ------------ Vector POINT ------------ POINT(1)=0.973906528517D0 POINT(2)=0.865063366688D0 POINT(3)=0. 679409568299D0 POINT (4)=0.433395394129D0 POINT(5)=0.148874338981D0 POINT (6)=-POINT (5) POINT(7)=-POINT(4) POINT(8)=-POINT(3) POINT(9)=-POINT(2) POINT(10)=-POINT(1) MPOINT=5 Vector BCOAL Vector BCOAL BCOAL (1) =0.2369268851D0 BCOAL(2)=0.4786286705D0

Print file "SZOT.FTN" BCOAL(3)=0.5688888888DO BCOAL(4) =BCOAL(2) BCOAL(5) =BCOAL(1) Vector BPOINT BPOINT (1) =0. 9061798459D0 BPOINT (2) =0. 5384693101D0 BPOINT(3)=0.DO BPOINT (4) =-BPOINT (2) BPOINT (5) =-BPOINT (1) Single integration NSP=31 RS1=0.9970874818iDO RS2=0.9846859096EDO RS3=0. 96250392509D0 RS4-0.93075699789D0 RS5-0.88976002994D0 RS6-0.83992032014D0 RS7-0.78173314841D0 RS8-0.71577678458D0 RS900.64270672292D0 RS10.56324916140D0 RS1100.47819378204D0 RS1200.38838590160D0 RS13=0.29471806998D0 RS14=0.19812119933D0 RS150.09955531215D0 RS16=0.DO XNS(1)=RS1 XNS (2)=RS2 XNS (3) =RS3 XNS (4) =RS4 XNS (5) =RS 5 XNS (6)=RS6 XNS (7) =RS7 XNS (8) =RS8 XNS (9)=RS9 XNS (10)=RS1O XNS (11)=RS11 XNS (12)=RS12 XNS (13) =RS13 XNS (14)=RS14 XNS(15)=RS15 XNS (16)=RS16 XNS (17)=-RS15 XNS (18)=-RS14 XNS (19)=-RS13 XNS (20)=-RS12 XNS (21)=-RS11 XNS(22)=-RS1O XNS(23)=-RS9 XNS(24) —RS8 XNS (25) — RS7 XNS (26) =-RS6 XNS(27)=-RS5 XNS(28) =-RS4 XNS (29) =-RS3 XNS(30)=-RS2 XNS(31)=-RS1 CNS (1) =0.0074708315792D0 Page 26

Print file "SLOT. FT N Page 2 7 CNS (2)-=0.0173186207903D0 CNS (3)=0. 0270090191849D0 CNS (4) =0. 0364322739123D0 CNS (5)=0.0454937075272D0 CNS (6) =0. 0541030824249D0 CNS (7)=0. 0621747865610D0 CNS (8)=0. 0696285832354D0 CNS (9)=0. 0763903865987D0 CNS (10)=0.0823929917615D0 CNS (11)=0. 0875767406084D0 CNS(12)=0.0918901138936D0 CNS (13) =0. 0952902429123D0 CNS (14)=0.0977433353863D0 CNS (15)=0. 0992250112266D0 CNS (16)=0.0997205447934D0 CNS (17)=CNS (15) CNS (18) =CNS (14) CNS (19)=CNS (13) CNS (20) =CNS (12) CNS (21) =CNS (11) CNS (22) =CNS (10) CNS (23) =CNS (9) CNS (24) =CNS (8) CNS (25)=CNS (7) CNS (26)=CNS (6) CNS (27) =CNS (5) CNS (28) =CNS (4) CNS (29) =CNS (3) CNS (30)=CNS (2) CNS (31) =CNS (1) 2) Double Integration NDP=16 R1=DSQRT((15.DO-2.D0*DSQRT(30.DO))/35.D0) R2=-R1 Sl=DSQRT ((15.DO+2.D0*DSQRT (30.DO) )/35.DO) S2=-Si A1=4.D0* (59.D0+6.D0*DSQRT (30.DO) ) /864.DO A2=4.D0*(59.D0-6.D0*DSQRT(30.D0)) /864.DO A3=4.D0*49.D0/864.DO XND(1,1) =R1 XND(1,2)=R1 CND (1) =A1 XND(2,1)=R2 XND(2,2) =R1 CND (2) =A1 XND(3,1) =R1 XND(3,2) =R2 CND (3) =A1 XND (4, 1) =R2 XND(4,2) =R2 CND (4) =A1 XND(5,1) =Sl XND(5,2)=Sl CND (5) =A2 XND(6, 1)=S1

Print file "SLOT.FTW' ag 2 Page 28 XND (6,2)S2 CND (6) =A2 C XND (7,1) =S2 XND (7,2)=Sl CND (7) =A2 C XND (8,1) =S2 XND (8,2) =S2 CND (8) =A2 C XND (9,1) =R1 XND (9,2) =Sl OND (9) =A3 C XND (10,1) =Rl XND (10,,2) =S2 CND (10) =A3 C XND (11,1) =S1 XND (11,2) =R1 CND (11) =A3 C XND (12,1) =S2 XND (12,,2) =R1 CND (12)=A3 XND (13,1) =R2 XND (13,,2)=S1 'CND (13) =A3.XND (14,1) =R2 XND (14,2)=S2 CND (14) =A3 XND (15,1) =S1 XND (15,2)=R2 C'ND (15) =A3 XND (16,1) =S2 XND (1 6,r2) =R2 CND (16) =A3.3) Triple Integration NTP=34 RS1=0.9317380000D0 RS2=-RS1 UU1=0. 916744177 9D0 UU2 =-UUJ1 SS1=0.4086003800DQ SS2=-SS1 TT1=0.739852 9500D0 TT2=-TT1 Bl=8.DO*0.035581808 96D0 B2=8.DO*0.012478 92770D0 B3=8.DO*0.0528 6772 991D0 B4=8.DO*0.02672752182D0 XNT (1,1) =RS1 XNT (1,2) =0.DO XNT (1,3) =0.DO CNT (1) =B1 XNT (2,1) =RS2

Print flie ",szoT.FTN"~ XNT (2, 2) =0.DO XNT (2, 3) =0.DO CNT (2) =Bl C XNT (3,1) =0.DO XINT (3,r2) =RS 1 XNT (3,3) =0.DO CNT (3) =Bl C XNT (4,1) =0.DO XNT (4,2) =RS2 XNT (4, 3)=0.DO CNT (4) =B1 XNT (5,1) =0.DO XNT (5,2) =0.DO XNT (5,3) =RS1 CINT (5) =Bl XNT(6, 1)=0.DO XNT (6,2) =0.DO XNT (6, 3)=RS2 C-NT (6)=B1 XNT (7,1)=-UU1 XNT (7,2) =UT1 XNT (7,3)=0.DO CNT (7) =B2 XNT (8,1) =UU2 XNT (8,2) =TU1 XNT (8,3) =0.DO CNT (8) =B2 XNT(9, 1)=UU1 XNT (9,2) =UU2 XNT (9, 3) =0.DO CNT(9)=B2 XNT (10,1) =UU2 XNT (10,2) =tU2 XNT (10,f3) =0.DO CONT (1 0) =B2 XNT (11, 1)=UUJ XNT (1 1, 2) =0.DO XNT (11,3) =UU1 CNT (11) =B2 XNT (12,1) =tU1 XNT (12,2) =0.DO XNT (12,3) =IU2 CNT (12)=B2 XNT (13, 1) =UU2 XNT (1 3, 2) =0. DO XNT (13,3)-U1U1 CNT (13) =B2 XNT (14,1) =UU2 XNT (14,2) =0.DO XNT (14,3) =UU2 CNT (14)=B2 XNT (15,1) =0.DO XNT ( 15,2)9 =UUT1 Page 29

print file ""SLOT. FTN" Pge3 Page 30 XNT (15, 3)=UU1 CNT (15) =B2 XNT (16,1) =0.DO XNT (16,2)=tJU1 XNT (1 6, 3) =UU2 ONT (16) =B2 XNT (17,1) =0.DO XNT (17,2) =UU2 XNT (17,3) =UU1 CNT (17) =B2 XNT (18,1) =0.DO XNT (1 8, 2) =tJU2 XNT (18, 3)=tJU2 CNT (18) =B2 XNT (19,1) =SS1 XNT (19,,2) =SS1 XNT (19,3) =SS1 CNT (19) =B3 XNT (20,1) =SS1 XNT (20,,2) =SS1 XNT (20,,3) -SS2 CNT (20) =B3 XNT (21, 1)=SS1 XNT(21,,2)=SS2 XNT (21,,3) =SS1 CNT (21) =B3 XNT (22,1) =SS1 XNT (22, 2) =SS2 XNT (22,,3) =SS2 CNT (22) =B3 XNT (23,1) =SS2 XNT (23,,2) =SS1 XNT (23,3) =SS1 ONT (23) =B3 XNT (24,1) =SS2 XNT (24, 2) =SS1 XNT (24, 3) =SS2 CNT (24) =B3 XNT (25,1) =SS2 XNT (25,,2) =SS2 XNT (25,,3) =SS1 CNT (25) =B3 XNT (26,1) =SS2 XNT (26,,2) =SS2 XNT (26,,3) =SS2 CNT (26) =B3 XNT (27,1l)=TT1 XNT (27,2) =TT1 XNT (27,3) =TT1 CNT (27) =B4 XNT (28,1) =TT1 XNT (28,,2) =TT1 XNT (28, 3)=TT2

Prin t fil nSLOT. FTJ Pag 3 Page 31 CNT (2 8) =B 4 XNT (29,1) =TT1 XNT (29,2)=TT2 XNT (29,3) =TT1 ONT (29) =B4 XNT (30,1) =TT1.XNT (30,2) =TT2 XKNT (30,3) =TT2 CONT (30) =B4 XNT (31, 1) =TT2 XNT (31,2) =TT1 XNT (31,3) =TT1 C'-NT (31) =B4 XNT (32,1) =TT2 XNT (3 2, 2) =TT 1 XNT (32,3) =TT2 CNT (32) =B4 XNT (33,1) =TT2 XNT (33,2) =TT2 XNT (33,3) =TT1 CNT (33) =B4 XNT (34,1) =TT2 XNT (34,2) =TT2 XNT (34,3) =TT2 CNT (34) =B4 RETURN END

Print file "POLES.FTN" Page 1 n*********************************************************************** The name of this file is........... POLES.FTN...... SUBROUTINE SPOLES IMPLICIT REAL*8(A-H,O-Z) ER:....Dielectric constant H:... Height of the dielectric substrate NE:....Number of TE surface waves NM:....Number of tm surface waves XS:....Matrix of poles contributing to TE surface waves XR:.... Matrix of poles contributing to TM surface waves ERR:....Error in the computation of the poles DIMENSION XS(40),XR(40),LOR(40) COMMON/DAT/ER, H, T, DLX, XC, YC, ZC, ERC, RMC, XQ, YQ, A, TPI, TPI2, PI, W, E1, *E2,EER,BKO,BK,BKK,FA,OFFSET (7),WDELTA,OFFLIM,ERROR,THI,FI, *NS1,NS2,NSS2,NOFF COMMON/DATT/COAL(20),POINT(20),CN(51),BM(51),TMP(20),TEP(20), *AM(41),DM(41),TPO(40),VXXM(20),VZXM(20),VZXE(20),BPOINT(10), *BCOAL(10),MPOINT,NPOINT,NKO,MA,NM,NE,NKOK,IFIRST AER=DSQRT (ER) ER2=ER*ER PI2=PI*PI MAXE=5 ERR=0. 0000001DO DP=H PART I: TE MODES AK0=2.D0*PI AK=DSQRT(ER)*AKO X0=DP*DSQRT (AK**2-AKO**2) WRITE (6,300) AK0,AK,XO,PI 300 FORMAT(10X,'AKO=',E14.7,2X,'AK=',E14.7,2X,'X0=',E14.7, *2X,'PI=',E14.7/) AN=XO/PI+0.5D0 NE=AN IF (NE.EQ.0) GO TO' 310 DO 2 I=1,NE IF (X0-(2.D0*FLOAT(I)+1.D0)*PI/2.D0) 3,3,4 4 XS0=(2.DO*FLOAT(I) -1.D0) *PI/2.DO+ERR XS1=(2.DO*FLOAT(I)+1.DO) *PI/2.D0-ERR GO TO 5 3 XS0=(2.DO*FLOAT(I)-1.DO)*PI/2.DO+ERR XS1=X0 5 CONTINUE IF (DABS(XSO-XS1)-ERR) 22,7,7 7 XSA=(XSO+XS1)/2.DO Y=-DTAN(XSA)*DSQRT(X0**2-XSA**2)-XSA IF (Y) 8,9,10 9 XS(I)=XSA GO TO 222

Print file "POLES. FTN" Page 2 8 10 22 222 2 C C C301 C C XS1=XSA GO TO 5 XSO=XSA GO TO 5 XS (I) = (XSO+XS1) /2.D0 XS (I) =DSQRT (AK* *2-XS (I) **2/DP **2) CONTINUE WRITE (6,301) ER, H FORMAT(//1OX,' Dielectric Constant=',D16.9/10X,'Substrate ' * 'Thickness',D16.9///) 310 IF (NE.EQ.0) WRITE (6,304) 304 FORMAT(/////10X,'No TE waves excited in IF (NE.EQ.0) GO TO 312 IF (NE.GT.0) WRITE (6,305) NE 305 FORMAT(///1OX,'There are',14, *' TE waves excited in the substrate'//) DO 302 I=1,NE TEP(I)=XS(I) WRITE (6,303) I,TEP(I) 303 FORMAT (10X, I4,2X,D16.9) TEP (I)=TEP (I) /AER 302 CONTINUE 312 CONTINUE the substrate'//) C C C C C c 15 14 16 17 19 301 1 301 I END OF PART I PART II: TM MODES AN=XO/PI+1. DO NM=AN DO 13 I=1,NM IF (XO-(2.DO*FLOAT(I)+1.DO)*PI/2.DO) 14,14,15 XS1=FLOAT(I) *PI-PI/3.DO-0.01D0 GO TO 16 XS1=X0 XSO=FLOAT (I-1) *PI+ERR CONTINUE IF (DABS (XS0-XS1)-ERR) 113,19,19 XRA= (XSO+XS1) /2.DO WRITE (6,301) XRA FORMAT(10X, 'XRA=',E14.7/) Y=DSQRT (ER) **2* (1.DO/DTAN(XRA)) *DSQRT (X0**2-XRA**2) -XRA IF (Y) 20,21,24 21 XR(I)=XRA GO TO 333 20 XS1=XRA GO TO 17 24 XSO=XRA GO TO 17 113 XR(I)=(XSO+XS1)/2.DO 333 XR(I) =DSQRT (AK**2-XR (I) **2/DP**2) 13 CONTINUE WRITE (6,307) NM 307 FORMAT(///1OX,'There are',14,' TM waves excited in the substrate'/ */) DO 308 I=1,NM TMP (I)=XR(I) WRITE (6,306) I,XR(I) 306 FORMAT (10X, I4,2X, D16.9) TMP (I)=TMP (I)/AER 308 CONTINUE

Print file "POLES.FTN" Page 3 322 CONTINUE NK=NE+NM IF (NE.EQ.0) GO TO 350 DO 411 IQW=1,NE TPO (IQW) =TEP (IQW) LOR(IQW)=1 411 CONTINUE 350 CONTINUE DO 412 IQW=1,NM TPO (NE+IQW) =TMP (IQW) LOR(NE+IQW)=0 412 CONTINUE IF (NK.EQ.1) GO TO 416 NNK=NK- 1 DO 415 IIP=1,NNK IK=IIP+1 DO 413 IIF=IK,NK QWR=TPO(IIP) IIW=LOR (IIP) IF (TPO(IIP).LT.TPO(IIF)) GO TO 413 TPO (IIP) TPO (IIF) LOR(IIP) =LOR(IIF) TPO(IIF) =QWR LOR(IIF)=IIW 413 CONTINUE 415 CONTINUE IF (LOR(1).EQ.0) IFIRST-0 IF (LOR(1).EQ.1) IFIRST=1 GO TO 417 416 IFIRST=2 417 CONTINUE RETURN END

Print file "CAV INV.FTN"F The name of this file is.........CAV INV.FTN...............CAVINV.FTN. It finds the inverse matrix for the case of scatering one or two C slots on the gound of a dielectric substrate. The slots have cavities at the back for minimization of the reflected field. ** *************** ************ ************ * ** ******** *** ****** ****** SUBROUTINE CAV INV(NOR) IMPLICIT REAL*8 (A-H,O-Z) COMPLEX CUR,BMATR,CAV MAT(60) COMPLEX ZS,ZSlS2,CJ,CCAV,CAV,CINC,SUMC,DETA,CIN,RUN,CS1, *CADD,SUMN COMMON/COMP/ZS(50),ZS1S2(300),NS,NSlS2 COMMON/DAT/ER,H,T,DLX,XC,YC,ZC,ERC,RMC,X0,Y0,A,TPI,TPI2,PI,W,El, *E2,EER,AK0,AK,AKK,FA,OFFSET(7),WDELTA,OFFLIM,ERROR,THI,FI, *NS1,NS2,NSS2,NOFF COMMON/MAN/BMATR(260,260), IA(260), IB (260) COMMON/BESSEL/BJO (6000) COMMON/CAV/RUN(1000),IRUN COMMON/LOSS/TLOS-E,TLOSM DATA CJ=(0.0,1.0) NOEL1=NS1 NOR=NS1 IDEL=20...... First Diagonal Matrix......... ARG1=PI*W/(2.DO*YC) ARG2=PI*YO/YC CALL VBJO(ARG1,ARG2) CALL CAVITY SIN=DSIN(2.DO*PI*DLX) CS1=SNGL((16.DO/120.DO)*DSQRT(ER)/(PI**3*XC*YC*RMC*SIN*SIN)) / * (1.0-CJ*SNGL(TLOS M)) C:CAV=CJ*CS1 C1=PI/XC IMIN=1 IMAX=NOEL1 DO 4 I=IMIN, IMAX IXN=0 DO 5 KI=I,IMAX IXN=IXN+1 XJ= (IXN*DLX+XO) *C1 XI= (KI*DLX+XO) *C1 SUMN=(0.0,0.0) DO 200 IN=1,IRUN SINI=DSIN(FLOAT(IN)*XI) SINJ=DSIN(FLOAT(IN)*XJ) CADD=RUN(IN) *SNGL (SINI*SINJ) SUMN=SUMN+CADD 200 CONTINUE CAV=CCAV*SUMN IF (IXN.EQ.1) CAV MAT(KI)=CAV BMATR(IXN, KI) =-ZS (I) +CAV BMATR(KI, IXN)=BMATR(IXN,KI) 5 CONTINUE Page 1

Print file "CAV INV. FTN" 4 CONTINUE 300 CONTINUE C DO 400 I=1,N51 WRITE (6,401) ICAVMAT(I),ZS(I) 401 FORMAT(10X,'I=',I4,2X,'CAV=',E14.7,2X,E14.7,2X, 'zs=',E14.7,2XrE14.7/) 400 CONTINUE C CALL MINVCD (NOR,NOR,DETA) RETURN END C THIS SUBROUTINE INVERTS A SQUARE COMPLEX MATRIX SUBROUTINE MINVCD (IA,MADETA) IMPLICIT REAL*8 (A-HO-Z) COMPLEX A,PIV,DETATEMP,PIV1 COMMON/MAN/A(2 60, 2 60), IR (2 60),IC (2 60) DO 1 I=1,MA IR(I)=O 1 IC(I)=0 DETA=(1.00,0.00) S=0.00 R=MA 2 CALL SUBMCD(IAIAMAMAI,J) -P IV=A (I, J) DETA=P IV*DETA Y=CABS (P IV) IF (Y.EQ.0) GO TO 17 IR(I)=J IC(J)=1 PIV=(1.00,0.00) /PIV A(I,J) =PIV DO 5 K=lMA 5 IF (K.NE.J) A(I,K)=A(I, K)*PIV DO 9 K=1,MA IF (K.EQ.I) GO TO 9 PIV1=A(K, J) 6 DO 8 L=1,MA 8 IF (L.NE.J) A (K, L) =A (K, L) -P IV1 *A (I,,L) 9 CONTINUE DO 11 K=1,MA 11 IF (K.NE.I) A(K,J)=-PIV*A(K, J) S=S+1.00 IF (S.LT.R) GO TO 2 12 DO 16 I=1,MA K=IC (I) M=IR (I) IF (K.EQ.I) GO TO 16 DETA=-DETA DO 14 L=1,MA TEMP=A(K, L) A(K, L)=A(I, L) 14 A(I,L)=TEMP DO 15 L=1,MA TEMP=A(L,M) A(L,M)=A(L, I) 15 A(L,I)=TEMP IC (M) =K IR(K)=M 16 CONTINUE RETURN 17 WRITE (6,18)IJ 18 FORMAT (10X,'MATRIX IS SINGULAR'/1OX,'I=',14,5X,'J=',14) RETURN Page 2

Print file "CAV INV.FTN" Page 3 END C........................................................................ SUBROUTINE SUBMCD(IA,JA,MA,NA, I,J) IMPLICIT REAL*8 (A-H,O-Z) COMPLEX A COMMON/MAN/A(260,260), IR(260), IC(260) I=0 J=0 TEST=0.00 DO 5 K=1,MA IF (IR(K).NE.0) GO TO 5 DO 4 L=1,NA IF (IC(L).NE.0) GO TO 4 X=CABS(A(K,L) ) IF(X.LT.TEST) GO TO 4 I=K J=L TEST=X 4 CONTINUE 5 CONTINUE RETURN END

Prin t f ilIe "CAVITY FTN" Pae Page I C C SUBROUTINE CAVITY IMPLICIT REAL*8 (A-H,O-Z) COMPLEX RUNCRNERCK2,CICOEF, SERMCPRNMCSQRCARGCTERMCCOTr * SERN COMMON/CAV/RtJN(1000)fIRUN COMMON/BESSEL/BJ0O(6000),COMMON/DAT/ERHTDLXXCrYCZCrERCRMCX0,Y0,ATPI,TPI2,PI,WEl. *E2,,EERIAK0,A~(,AKKFAOFFSET (7),WDELTAOFFLIM,ERROR,THI,FI, *NS1,NS2,NSS2,NOFF C-OMMON/LOSS/TLOSETLOS M CI= (0.0,1.0) ERROR=1.D-4 PI=3. 14159265358 9D0 SERN= (0.0,0.0) NTEST=0 1)0 1 N=1, 1000 EN=1.DO RN=F LOAT (N) / (2. DO*XC) RN2=RN*RN RNP1=1.D0+RN RNM1=1.DO-RN RNPM=RNP 1 *RNM1 RNPM2=RNPM*RNPM CK2=(1.0-CI*SNGL(TLOSE) )*(1.0..CI*SNGL(TLOSM) )*SNGL(ERC*RMC * /ER) CRNER=CK2 -SNGL (RN2) SIN1=DSIN (PI*DLX*RNP1) SIN2=DSIN (PI*DLX*RNMl) SINP=SIN1 *SIN2 SINP2=SINP *SINP COEF=CRNER*SNGL (S INP2 *EN/RPpM2) SERM= (0.0,0.0) ITEST=0 M1=0 M1=M + 1 M=M1 -1 EM=0.5D0 IF(M.GT.0) EM=1.DO BJ=BJO (Ml) RM=FLOAT (M) /(2.DO*YC) RM2 =RM*RM CRNM=CRNER-SNGL (RM2) CABS 1=CABS (CRNM) IF (CABS1.LT.1.E-8) GO TO 100 CSQR=CSQRT (CRNM) CARG=SNGL (2.DO*PI*ZC) *CSQR IF (CABS1.LT.50.0) THEN CTERM=SNGL(EM*BJ) *(CCOS (CARG) /(CSQR*CSIN(CARG))) ELSE RARG=REAL (CARG) AARG=AIMAG (CARG) COT=DCOS (RARG) /DSIN (RARG) EX=DEXP (-2. DO *AARG) COTH= (1. DO0+EX) / (1. DO0-EX) CCOT= (-1. 0-CI*SNGL (COT*COTH) )/ (SNGL (COT) - * CI*SNGL(COTH)) CTERM=SNGL (EM*BJ) *CCOT/CSQR END IF SERM=SERM+CTERM RATIO=CABS (CTERM/SERM) 2

Print f I0 "CAvITY FTN"Page2 Page 2 IF (RATIO.LT.ERROR) GO TO 5 ITEST=0 GO TO 2 5 ITEST=ITEST+1 IF (ITEST.LT.4) GO TO 2 C C WRITE (6,10) N,M, SERM C 10 FORMAT (//6X,'Maximum M for a given N:',2X,'N=', C * 14,3X,' Mmax=', I4,2X,' Serm=', E14.7,1X,E14.7) C RUN (N) =COEF*SERM WRITE (6,40) NRUN(N) 40 FORMAT(10X,'N=',14, 3X,'RUN=',(E14.7,2X,E14.7)) SERN=SERN+RUN (N) RATIO=CABS (RUN (N) /SERN) IF (RATIO.LT.ERROR) GO TO 8 NTEST=0 GO TO 1 8 NTEST=NTEST+1 IF (NTEST.LT.4) GO TO 1 IRUN=N GO TO 101 1 CONTINUE GO TO 101 C 100 WRITE (6,20) 20 FORMAT(///1OX,'..........WARNING!........'/10X,'The eigenvalue Kz' became equal to 0'!!/) 101 CONTINUE RETURN END This function evaluates the zeroth order first kind Bessel Function JO SUBROUTINE VBJO(ARG1,ARG2) IMPLICIT REAL*8 (A-H, O-Z) COMMON/BESSEL/BJO (6000) 21=3.141592653589D0 1)0 1 M1=1, 6000 M=M1-1 X=FLOAT (M) *ARG1 X1=FLOAT (M) *ARG2 COS1=DCOS(Xl) COS2=COS 1 *COS 1 IF (X.GT.0.001D0) GO TO 10 X3=X/3.DO X32=X3 *X3 X3 4=X32 *X32 X3 6=X3 4 *x32 BSJO=1.DO2.2499997D0*X32+1.2656208D0*X34-0.3163866D0 **x36 BJO(Ml)=BSJO*COS2 GO TO 1 10 IF (X.GT.3.DO) GO TO 12 X3=X/3.DO X32=X3 *X3 X3 4=X32 *x32 X36=X34*X32 X38=X36*X32 X310=X38*X32 X312=X310*x32 BSJO=1.DO-2.2499997D0*X32+1.2656208DO*X340.3163866D0

Print file "CAVITY. FTN" * *X36+0.0444479DO*X38-0.0039444DO*X310+0.00021000 * DO*X312 BJO(M1)=BSJO*COS2 GO TO 1 12 CONTINUE X3=3.DO/X X32=X3*X3 X33=X32*X3 X34=X33*X3 X35=X34*X3 X36=X35*X3 FJ0=0.79788456D0-0.00000077D0*X3-0.00552740D0*X32-0.0000 * 9512D0*X33+0.00137237D0*X34-0.00072805D0*X35+0.00014 * 476D0*X36 TJO=X-0.78539816D0-0.04166397D0*X3-0.00003954D0*X32+0.00 * 262573D0*X33-0.00054125D0*X34-0.00029333D0*X35+0.000 * 13558D0*X36 WCON=DSQRT(1. DO/X) BSJO=WCON*FJO*DCOS(TJO) BJO(M1)=BSJO*COS2 1 CONTINUE RETURN END * *** * *********************** ****************** *** ******* ***** Page 3

Print file "MULT ^XC.FTN" The name of this file is: MULT-EXC.FTN This subroutine evaluates the excitation vector and multiplies it with the inverse of the impedance matrix to get the unknown field on the slot SUBROUTINE MULT EXC(NOR,R0) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION Y(260),PHASE(260) COMPLEX CUR,BMATR,CJ,CCAV,CAV,CINC, SUMC,CIN,R0 COMMON/DAT/ER,H,T,DLX,XC,YC,ZC,ERC,RMC,XO,Y0,A,TPI,TPI2,PI,W,E1, *E2,EER,KAK0,AK,AKK,FA,OFFSET(7),WDELTA,OFFLIM,ERROR,THI,FI, *NS1,NS2,NSS2,NOFF COMMON/MAN/BMATR(260,260), IA(260), IB(260) COMMON/PAT/CUR(260) DATA CJ-(0.0,1.0) RN0-120.D0*PI RNsRNO/DSQRT (ER) RK2.DO*PI RKO=RK/DSQRT (ER) SINI=DSIN (THI) COS IDCOS (THI) SINT=SINI/DSQRT (ER) SINT2=S INT*SINT COST=DSQRT (1.DO-SINT2) ARGT=RK*COST*T ARGH=RK*COST*H ARGL1=RK*DLX ARGL2=ARGL1 *S INT D1-4.DO*COSI*DCOS (FI) *DCOS (ARGT) * (DCOS (ARGL2) -DCOS (ARGL1)) * /(RK*(1.D0-SINT2)) DC=RN0*COST*DCOS (ARGH) DS=RN*COSI*DSIN (ARGH) CINC=SNGL (D1) / (SNGL (DC) +CJ*SNGL(DS)) ARG=ARGH DC=RNO*COST DS=RN*COSI*DTAN (ARG) R0=- (SNGL (DC) -CJ*SNGL (DS)) / (SNGL (DC) +CJ*SNGL (DS)) DO 70 IQ-1,NOR SUMC=(0.0, 0.0) 401 DO 170 JQ=1,NOR ARG=ARGL2*FLOAT (JQ) +X0 EC=DCOS (ARG) ES=DSIN (ARG) CIN=CINC* (SNGL (EC) -CJ*SNGL (ES)) WRITE (6,778) IQ,JQ,CIN,BMATR(IQ,JQ) 778 FORMAT (2X, I4,2X, I4, 2X, E14.7,2X, E14.7,2X, E14.7,2X, * E14.7) SUMC=SUMC-BMATR (IQ, JQ) *CIN 170 CONTINUE CUR (IQ) =SUMC 70 CONTINUE WRITE (6,71) 71 FORMAT (///10X,'Electric Field Distribution on the Slots',////) IMIN-1 IMAX-NOR Page 1

Print file MULLT _EXC.FTN" Page 2 Do 76 IQQ-IMINIMAX RECUR1-REAL(CUR(IQQ)) Y(IQQ)=CABS (CUR(IQQ)) AICUR1=AIMAG (CUR( IQQ)) PHCUR1=ATAN2 (AICURi, RECUR1) PHASE(IQQ)=180.OO*PHCURl/PI IF (IQQ.EQ.1) WRITE (6,77) 77 FORMAT(///1OX,'Electric Field on the first Slot',///) IF (IQQ.EQ.(NS1+1)) WRITE (6,78) 78 FORMAT(///1OX,'Electric Filed on the Second Slot',///) C WRITE (6,81) IQQCUR(IQQ),Y(IQQ),PHASE(IQQ) 81 FORMAT (lX, 'C(',14,' )=(',E14.7,', ',E14.7,')',2X, * ( El14. 7 ', ',El4. 7 ' ) ' C C WRITE (1,92) IQQCUR(IQQ),Y(IQQ),PHASE(IQQ) C 92 FORMAT(I5,4E15.3) C 76 CONTINUE RETURN END

?rint file "RErV7rC.F FPge PageI The name of this file is......REFLEC.FTN It computes the reflection coefficient for cavity-backed slots SUBROUTINE REFLEC(WREALRO, REFL) IMPLICIT REAL*8 (A-HO-Z) COMPLEX ETHEFICURREFL,RRRSRA COMMON/DAT/ER, H, TDLX, XC, YC, ZC, ERC, RMCXO, YOA, TPI, TPI2, PI, *WIDTH, El, E2, EER, AKO, AK, AKK, FA, OFFSET (7), WDELTA, OFFLIM, ERROR, *THIFI,NS1,NS2,NSS2,NOFF C:OMMON/PAR/NSLOTALONG(7),S(20,2),NSER COMMON/ PAT / CUR (2 6 0) NSLOT=NOFF NSER=10 U- (WREAL/WIDTH) U=DATAN(DSQRT(l.DO/ (U*U)-l.DO)) Ul=U/FLOAT (NSER) DO 3 JN=lNSER S2=2.DO*FLOAT (JN) -..DO S2=S2/ (2. DO*FLOAT (NSER)) S3=COS(S2*U) S (JN, 2) =S3*WIDTH/2.DO S (JN, 1)=Ul 3 CONTINUE ATH1=WREAL /WIDTH ATH1=DATAN (DSQRT(1. / (ATH1*ATHl)-1l.0)) ATH2=P I-ATH1. CURIN= (1.DO/PI) * (ATH2-ATH1) CALL EFIELD (ETH, EFI) RS=(ETH*SNGL(DCOS(THI)*DSIN (F))+EFI*SNGL (DCOS(F)))* *SNGL(CURIN) RA=2. 0*RS/SNGL (DSQRT (ER)) REFL=RO +RA WRITE (6,100) THIFIRORA 100 FORMAT(10X,'THI=',El4.7,2X,'FI=',El4.7/1OX,'RO=',El4.7,1X,E14.7 */10X,'RA=',E14.7,2X,E14.7) RETURN END The name of this subroutine is: EFIELD It evaluates the far field of a dipole SUBROUTINE EFIELD(ETH,EFI) IMPLICIT REAL*8(A-HO-Z) COMPLEX WETHEFIXNUMXDENF1XZNUMZDENFlZ,PTH,PFI, *SUMC, WEXP, CUR, COEF COMMON/DAT/ERHTDLXXC, YCZCERCRMCXOYOA,TPITPI2,PI, *WIDTHElE2,EERAKOAKAKKFAOFFSET(7),WDELTAOFFLIMERROR, *THIFINS1,NS2,NSS2,NOFF COMMON/PAR/NSLOT, ALONG (7),S (20, 2), NSER COMMON/PAT/CUR (260) COMMON/BO 1 /BJO, BJ1

Print file "REupY7C.FTN" Page 2 CKK=2.DO*PI CKO=CKK/DSQRT (ER) SINFT-DSIN(FI)*DSIN(THI) ARG=(CKO*WIDTH/2.DO) *SINFT CALL BSJ0(ARG) SSUM=0.DO DO 5 JN=1,NSER ARAF=CKO * S (JN, 2) *S INFT CAFF=DCOS (A.AF) SSUM=SSJM+S (JN, 1)*CAFF 5 CONTINUE TERMI= (BJO-SSUMM*2.DO/PI) ERTH=DSQRT(ER-DSIN(THI)**2) ERH=CK0 *ERTH*H SINH=DS IN (ERH) SINTH=DSIN(THI) SINFI=DSIN (FI) COSH=DCOS (ERH) COSTH=DCOS(THI) COSFI=DCOS(FI) W=(0 0, 1 0) XNUM=SNGL (COSTH) XDEN=SNGL (-ERTH*SINH) +W*SNGL (ER*COSTH*COSH) FlX=-SNGL((ER! (4.DO*PI))) )*XNUM/XDEN ZNUM=SNGL (SINH*SINTH*COSTH) ZDEN=XDEN* (-SNGL (COSTH*SINH) +W*SNGL (ERTH*COSH)) F1Z=SNGL((1.DO-ER) /(4.DO*PI)) *ZNUM/ZDEN PTH=SNGL (COSFI)*(FlX*SNGL(COSTH) FlZ*SNGL(SINTH)):PFI=SNGL(SINFI)*F1X For the single slot SUMC=(0.0,0.0) DO 1 I=1,NSLOT R8=CKO* (FLOAT (I) *DLX) *SINTH*COSFI SUMC=SUMC+CUR(I) *(SNGL(DCOS(R8))+SNGL(DSIN(R8) ) *W) 1 CONTINUE R8=CKO*ALONG(1)'*SINTH*COSFI R9=CKO*OFFSET(1)*SINFI*SINTH WEXP=(SNGL(DCOS(R8) )+ 3NGL(DSIN(R8) ) *W) *(SNGL(DCOS(R9) )+SNGL( *DSIN(R9) ) *W) COEF=WEXP * SNGL (TERMI) PTH=PTH *COEF PFI=PFI *COEF IF (ABS (FI).GT.1.E-4) GO TO 2 THER=ABS (ABS (THETA) -PI/2. 0) IF (THER.GE.l.E-4) GO TO 2 IF (ABS(EER-1.00).LT.1.E-6) GO TO 3 2 RNUM=DCOS (CKO*DLX*SINTH*COSFI) -DCOS (CKK*DLX) RDEN=DSIN (CKK*DLX) * (1.DO- (CKO*SINTH*COSFI/CKK) **2) RAT IO=RNUM/ RDEN ETH=PFI*SNGL (RATIO) EFI=PTH*SNGL (RATIO) ETH=ETH*SUMC EFI=EFI *StMC WRITE (6,4) ETHEFI 4 FORMAT (5X,'ETH=', (E14.8,1X,E14.8),5X,'EFI=', *(E14.8,1XE14.8)!) RETURN

file FFRBzE97C FTNr7W" IPage 3 RATIO=DLX*CKO/2 0 ETH=PFI*SNGL (RATIO) *SUMC EFI=-PTH*SNGL (RATIO) *SUMC RETURN END