NASA Technical Memorandum 101614 

CWG = A Fortran Program for 

Mutual Coupling in a Planar 

Array of Circular Waveguide-fed 

Apertures 

M. C. Bailey 

June 1989 

(lfASA-TI!l-'lOl614) C@G: A FORTBA8 PROG3Al FOR 

MIOIFUIL COWLIN6 I# A PLUlAR aRRAY OF 

Langley Research Center) 46 p CSCL 171 

CIRCULAR PAVEGUIDE-PBD APERTURES IHASL , 

National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23665-5225 

\ 889-23063 

Unclas 

G3/32 0224584 
. Introduction 

Circular waveguide-fed apertures or conical horns are often 

used as elements of planar array antennas. The purpose of this 

present paper is to document a computer code which calculates the 

mutual coupling between circular apertures in a conductive plane. 

The program is quite general in that the apertures do not have to 

be the same sizes, nor do they have to be polarized in the same 

direction. In addition, several waveguide modes (TE and/or TM) 

may be specified in the apertures and the mutual coupling between 

all combinations of apertures and modes will be calculated. The 

program also allows multiple layers of homogeneous dielectrics to 

be placed over the aperture array. Outside the layered region, 

one can specify either a homogeneous half-space or a perfect 

reflecting surface. The program is a revision of an earlier one 

(ref.1) in an attempt to overcome some difficulties encountered 

with the earlier version and to expand its versatility. 
Symbols 

ai 

a 

d 

J 

n 

EJ 

JJ 2 1 

siJ 

[SI 

TE 

TM 

XlYlZ 

Radius of ith aperture. 

Radius of jth aperture. 

Thickness of nth homogeneous dielectric layer. 

aperture electric field polarization vector. 

aperture electric field polarization vector. 

ith 

j th 

Bessel function of the first kind of order m and 

argument z. 

Derivative of J(z) with respect to 2. m 

d=T 

Wave propagation constant in free space, 2n/h. 

Wave propagation constant in nth layer. 

Wave propagation constant in homogeneous half-space. 

First index of waveguide mode in ith aperture. 

First index of waveguide mode in jth aperture. 

Second index of waveguide mode in ith aperture. 

Second index of waveguide mode in jth aperture. 

Number of homogeneous layers outside array plane. 

Center-to-center spacing between apertures. 

Reflection coefficient for plane wave with electric 

field polarized parallel to plane of incidence. 

Reflection coefficient for plane wave with electric 

field polarized perpendicular to plane of incidence. 

Plane wave reflection coefficient for parallel 

polarization for exterior scattering problem. 

Plane wave reflection coefficient for perpendicular 

polarization for exterior scattering problem. 

Mutual coupling between ith and jth apertures for 

one mode in each aperture. 

Complex square matrix whose elements are S . Denotes transverse (to z) electric field. 

Denotes transverse (to z) magnetic field. 

Cartesian coordinate variables. 

1J 

2 
xi 

y i 

YJ 

X 

J 

J 

cy1 

cy01 

13 

6 

e 

0 

n E 

%+1 

& 

P 

x-coordinate of ith aperture center. 

x-coordinate of jth aperture center. 

y-coordinate of ith aperture center. 

y-coordinate of jth aperture center. 

Mutual admittance between ith and jth apertures 

for one mode in each aperture. 

Complex square matrix whose elements are Y . Complex diagonal matrix whose elements are the 

waveguide modal characteristic admittances. 

Normalized radial propagation constant in 

cylindrical coordinate system. 

Distance from outer surface of Nth dielectric layer 

to exterior scattering boundary. 

Permittivity of free space. 

Permittivity of region just outside aperture plane. 

Permittivity of nth homogeneous layer. 

Permittivity of homogeneous half-space outside of 

Equivalent permittivity just inside of exterior 

scattering boundary. 

Wavelength in free space. 

Permeability of free space. 

Permeability of region just outside of aperture. 

Permeability of nth homogeneous layer. 

Permeability of homogeneous half-space outside of 

Equivalent permeability just inside of exterior 

scattering boundary. 

Polarization of ith aperture with respect to y-axis. 

Polarization of jth aperture with respect to y-axis. 

iJ 

layer. Nt h 

layer. Nth 

#j- + 

zero of J (2). 

zero of J'(z). 

nth m 

nth m 

3 
Theory 

The analytical formulation for the electromagnetic 

interaction between apertures in a planar array has been developed 

earlier (ref.1) and will not be repeated here. 

the apertures in the array is obtained from the appropriate 

coefficients of the scattering matrix, [SI, as follows: 

The mutual interaction between the various modes among all 

L AL 

Where [YO] is a diagonal matrix whose elements are the 

characteristic admittances of the waveguide modes and [Y] is the 

admittance matrix whose elements are the mutual admittances 

between each waveguide mode and all modes in each waveguide 

aperture. Therefore, for example, if the array consisted of 5 

apertures and each aperture field was assumed to be the 

superposition of 7 waveguide modes (TE and/or TM), then (1) would 

consist of square complex matrices each of size 35 by 35. The 

sizes of the matrices can expand very rapidly if one does not use 

some "engineering judgment" in selecting the modes to approximate 

the aperture fields. The expressions for calculating the elements 

of the admittance matrix (i.e., the complex mutual admittance 

between pairs of waveguide modes) are derived in reference 1. 

where , 

Tm = 1 for mi= 0 

i 

=2 otherwise. 

4 
Wl(@) and W2((3) are defined as 

where (for real angles) (3 can be interpreted as the sine of the 

angle of incidence for a plane wave with RI1(@) and R1(/3) being the 

reflection coefficients in the outward direction evaluated at the 

aperture boundary for parallel and perpendicular polarization at 

an air-dielectric interface. Notice from equation 2 that the 

mutual admittance is a weighted sum over all angles of incidence, 

including real space (Oq3sl) and invisible space (@>1). An 

alternate form for Wl(p) and W2(p) is 

where gl(p) and fl(@) are the solutions to the wave equations in 

the external region evaluated at the aperture plane, and g;(p) and 

f;(p) are the corresponding derivatives in the direction normal to 

the aperture plane. 

permeability of the external medium evaluated just outside the 

aperture plane. 

By applying the boundary conditions of continuity of 

tangential fields at the interface of each homogeneous dielectric 

layer, one can evaluate the mutual admittance between aperture 

modal fields with a multilayer dielectric medium outside (as shown 

in figure 1). 

c1 and 1-1, are the permittivity and 

5 
X 

J X 

1 

Figure la. Coordinate geometry for the i th and j th elements of 

a planar array of circular waveguide-fed apertures. 

E 

2 

p2 

d2 

tj 

... 

... 

... 

E 

N 

PN 

E 

N+l 

%+l 

E 

P 

pP 

Figure lb. Cross-section of i th and j th circular waveguide-fed 

apertures radiating into a multi-layer dielectric 

with an exterior arbitrary reflecting boundary. 

6 
Applying boundary conditions at layer interfaces, results in the 

following recursive equations, 

for n = N,(N-l),(N-2), ..., 3,2,1. where k is defined as 

Then starting with 

. 

equations 7 and 8 are evaluated successively from the outermost 

homogeneous layer to the layer just outside of the aperture plane. 

7 
In the special case of a perfect electric conductor of 

infinite extent on the outer surface of the Nth layer, fN+l(f3)=0 

and g;+,(@)=O; therefore, the starting conditions for equations 7 

and 8 become 

If the medium outside of the Nth dielectric layer is not a 

homogeneous half-space, the starting conditions (equation 9) are 

replaced by 

= -j K(P) 

E 

P 

‘N+1 

where , 

RI(p) and Rf;(P) are the plane wave reflection coefficients, 

evaluated at the external scattering boundary (shown dashed in 

figure lb) for perpendicular and parallel polarization of the 

electric field vector with respect to the plane of incidence. E 

and p are the equivalent permittivity and equivalent permeability 

just inside of the external scattering boundary. 

P 

P 

The starting 

8 
conditions in equation 11 allow for future modifications of the 

computer code to include problems for which the plane wave 

reflection is known from solutions to exterior scattering 

problems, such as objects in the vicinity of the array or an 

inhomogeneous ionized plasma. 

The other quantities in equation 2 are related to the Fourier 

transforms of the aperture modal fields as derived in reference 1. 

For the TE modes in the ith aperture: 

-1 -- - / 

For the TM modes in the ith aperture: 

For the jth aperture, the i subscripts are replaced with j in 

equations 12-15. 

If R=O, then the ithand jth apertures are coincident and 

equation 2 is the mutual admittance between modal fields in the 

same aperture. When R=O, the orthogonality of the modal functions 

results in U (p)=V (p)=O except when mi=m . u (p) and V (p) are as given in Table I. 

When R=O and mi=mJ; 

iJ iJ J 

iJ iJ 

9 
Table I 

mj ( R=O; m = 

I 

Modes 

TEI TEJ 

TMI TMJ 

TEl TMJ 

TMi TEJ 

Defining the following terms: 

[ 2/rm1) cos 

0 

0 

0 

R = J(XJ-xJ2 + (YJ-YJa 

U (p) and Vij(P) for R>O are given in Table 11. 

iJ 

10 
I 

Modes - 

TEi TEj 

TMI TMJ 

TEl TMj 

TMI TEj 

m 

(-1) ( s, B+(B) 

m 

(-1) ( s, B+(P) 

References 

0 

0 

1. M. C. Bailey: "Analysis of Finite-size Phased Arrays of 

Circular Waveguide Elementstt, NASA Technical Report, R-408, 

April 1974. 

11 
Appendix 

Computer Program 

+ STE I H BES JS 1 

~STEJ H BESJS I 

The computer code is written in FORTRAN IV and uses double 

precision variables throughout, except for integers. It was found 

that double precision was necessary in order to obtain sufficient 

accuracy on a VAX computer. The input data is read from logical 

unit 5. Results are output to logical unit 6. A flow chart 

indicating the subroutine calls is given and also a description of 

the input parameters with sample input and output data. The 

computer code has been verified against results in reference 1. 

Subroutine Flowchart 

..I.. ....................................................... )(. 

................ i . IF end Of ..i..... * .............. INPUT j input data! (stop) 

I ............................................................ 

I e MLAYER 

I 

A- 1 
Definition of Input Data 

NHOLE 

NMODE 

NUMTE 

NUMTM 

MIJ(1) 

NIJ(1) 

MIJP( I) 

NIJP( I) 

AIJ(1) 

XIJ(1) 

YIJ(1) 

PHIJP( I) 

F 

FCON 

ER 

CEP 

CUP 

NLAY 

Number of apertures in array. 

Number of modes per aperture. 

Number of TE-modes per aperture. 

Number of TM-modes per aperture. 

(NOTE: same modes assumed for all apertures). 

First index of I-th TE mode. 

Second index of I-th TE mode. 

First index of I-th TM mode. 

Second index of I-th TM mode. 

Radius of I-th aperture. 

x-coordinate of center of I-th aperture. 

y-coordinate of center of I-th aperture. 

Polarization of I-th aperture (degrees CCW). 

Frequency (Hertz). 

Factor for converting input to centimeters. 

Dielectric constant of waveguide interior. 

Relative epsilon outside of layered region. 

Relative mu outside of layered region. 

Number of homogeneous dielectric layers. 

Thickness of I-th homogeneous layer. 

Relative epsilon of I-th layer. 

Relative mu of I-th layer. 

A-2 
Input Data 

( 1 ) Read NHOLE, NMODE, NUMTE, NUMTM 

Check for end-of-file on unit 5. 

(2) If NUMTE=O, skip to (3). 

Read MIJ(I),NIJ(I) ; 14 To NUMTE 

(3) If NUMTM=O, skip to (4). 

Read MIJP(I),NIJP(I) ; 1=1 To NUMTM 

(4) Read AIJ(I),XIJ(I),YIJ(I),PHIJP~I) ; 14 To 

(5) Read F,FCON,ER 

(6 ) Read CEP, CUP 

(7) Read NLAY 

(8) If NLAY=O, skip to (9). 

Read D(I),CE(I),CU(I) ; I=l To 

(9) Perform calculations. 

Read next set of input data (1 

NLAY 

A-3 

3LE 
Sample Input Data 

2110 

11 

0.75 0.0 0.0 0.0 

0.75 0.0 2.5 0.0 

6.OD09 2.54 1.0 

1 

0.18 (2.6, -0.0156) (1 

(1.0,O. 0) (1.0,O. 0) 

0,o. 0) 

Sample Output Data 

........................................................................ 

MUTUAL COUPLING OF CIRCULAR APERTURES 

RADIATING INTO MULTI-LAYERS 

NUMBER OF LAYERS = 1 

LAYER TH I CKNESS COMPLEX EPSILON COMPLEX MU 

1 0.1800 ( 2.60000, -0.01560) ( 1.00, 0.00) 

OUTSIDE HALF-SPACE EPSILON = 1.00000 0.00000 

OUTSIDE HALF-SPACE MU - 1.00000 0.00000 

**** INPUT DIMENSIONS IN INCHES   

** ** FREQUENCY = .600000E+10 HERTZ   

++++++++ APERTURE ARRAY GEOMETRY +++++++++ 

HOLE RADIUS X Y POL. 

1 0.7500 0.000 0.000 0.000 

2 0.7500 0.000 2.500 0.000 

MODE 1 = TE-11 

CY( 1, 1)=( 0.3415E-02, 0.1691E-02) BB= 17.00 IH= I IM= 1 JH= 1 JM= 1 

CY( 1, 2)=( 0.3443E-04,-0.3158E-03) BB= 15.00 IH= 1 IM= 1 JH= 2 JM= 1 

YMN ( 1 )=( 0.1695E-02, O.OOOOE+OO ) IHOLE= 1 IMODE= 1 

YMN( 2)=( 0.1695E-02, O.OOOOE+OO) IHOLE= 2 IMODE= 1 

SCATTERING MATRIX 

S( 1, l)=(-O. 4036E+00, -O.l964E+OO) -6.9570 DB -154.0525 DEG 

S( 1, 2)=( 0.1871E-01, 0.3199E-01) -28.6224 DB 59.6795 DEG 

S( 2, 1 )=( 0.1871E-01, 0.3199E-01) -28.6224 DB 59.6795 DEG 

S( 2, 2)=(-0.4036E+00,-0.1964E+00) -6.9570 DB -154.0525 DEG ........................................................................ 

A-4 
c** ***ff*********** CWG ............................................... C Mutual admittance and scattering matrix calculation for * 

C circular waveguide array with TE-MN and TM-MN modes. * 

C NOTE: Field equations of Marcuvitz show * 

C TM modes polarized 90 degrees to TE modes. * 

C program modified to rotate polarization of TM modes * 

C by PI/2 such that TE-11 and TM-11 are polarized the same.* 

C * 

C Program can handle multiple homogeneous layers over apertures. * 

C * 

C Hampton, VA (1989) * C*********************************************~*********ff #** ********* 

102 

C (INPUT on unit #5) (OUTPUT on unit #6) * 

C PROGRAM: M.C. Bailey * 

FORMAT ( lX, 79 ( lH+ 1 1 

WRITE(6,102) 

CALL CONSTS 

100 WRITE(6,lOZ) 

CALL INPUT 

CALL YMAT 

CALL SMAT 

GO TO 100 

END 

A-5 
SUBROUTINE CONSTS 

IMPLICIT REAL*8 (A, B, D-H, 0-Z) 

IMPLICIT COMPLEX*16 (C) 

COMMON /ZEROS/ XMNP(9,6),XMN(9,6) 

COMMON /CONST/ C J , PI, TWOP I, DTOR, RTOD, BIG, SMALL, BMAX, YO, MAX, LAYM 

DATA XMNP ( l,l), XMNP ( 1,2 1, XMNP ( 1,3 )/3.831706,7.0155867,10.173468/ 

DATA XMNP (2.1 1, XMNP(2,2 1, XMNP (2,3 )/1.84118,5.33144,8.53632/ 

DATA XMNP(3,1),XMNP(3,2),XMNP(3,3)/3.05424,6.70613,9.96947/ 

DATA XMNP(4,1),XMNP(4,2),XMNP(4,3~/4.20119,8.01524,11.34592/ 

DATA XMNP(5,1),XMNP(5,2),XMNP(5,3~/5.31755,9.28240,12.68191/ 

DATA XMNP(6,l). XMNP(6,2), XMNP(6,3)/6.41562,10.51986,13.98719/ 

DATA XMNP(7,1),XMNP(7,2),XMNP(7,3~/7.50127,11.73494,15.26818/ 

DATA XMNP(8,1), XMNP(8,2), XMNP(8,3)/8.57784,12.93239,16.52937/ 

DATA XMNP(9,1), XMNP(9,2), XMNP(9,3)/9.64742,14.11552,17.77401/ 

DATA XMNP ( 1,4 1, XMNP ( 1,5 1, XMNP ( 1,6 )/13.323692,16.47063,19.615859/ 

DATA XMNP(2,4),XMNP(2,5),XMNP(2,6)/11.70600,14.86359,18.01553/ 

DATA XMNP(3,4), XMNP(3,5), XMNP(3,6)/13.17037,16.34752,19.51291/ 

DATA XMNP(4,4), XMNP(4,5), XMNP(4,6)/14.58585,17.78875,20.97248/ 

DATA XMNP(5,4), XMNP(5.5). XMNP(5,6)/15.96711,19.19603,22.40103/ 

DATA XMNP(6,4),XMNP(6,5),XMNP(6,6)/17.31284,20.57551,23.80358/ 

DATA XMNP(7,4),XMNP(7,5),XMNP(7,6)/18.63744,21.93172,25.18393/ 

DATA XMNP(8,4),XMNP(8,5),XMNP(8,6)/19.94185,23.26805,26.54503/ 

DATA XMNP ( 9,4 1, XMNP (9,5 1, XMNP ( 9,6 ) /21.22906,24.58720,27.88927/ 

DATA XMN ( 1, 1 ) , XMN ( 1,2), XMN ( 1,3)/2.4048256,5.5200781,8.6537279/ 

DATA XMN(2, l), XMN(2,2), XMN(2,3)/3.8317060,7.0155867,10.1734681/ 

DATA XMN(3,1), XMN(3,2), XMN(3,3)/5.1356223,8.4172441,11.6198412/ 

DATA XMN(4,1),XMN(4,2),XMN(4,3)/6.3801619,9.7610231,13.0152007/ 

DATA XMN(5,1),XMN(5,2),XMN(5,3)/7.5883427,11.0647095,14.3725367/ 

DATA XMN(6, l), XMN(6,2), XMN(6,3)/8.7714838,12.3386042,15.7001741/ 

DATA XMN(7,1),XMN(7,2),XMN(7,3)/9.93611,13.58929,17.00382/ 

DATA XMN (8,l) , XMN (8,2 1, XMN (8,3 ) /11.08637,14.82127,18.28758/ 

DATA XMN(9,1),XMN(9,2),XMN(9,3)/12.22509,16.03777,19.55454/ 

DATA XMN ( 1,4), XMN ( 1,5), XMN ( 1,6 )/11.7915344,14.9309177,18.071064/ 

DATA XMN(2,4),XMN(2,5),XMN(2,6)/13.32369,16.47063,19.61586/ 

DATA XMN (3,4), XMN (3,5), XMN (3,6 ~/14.79595,17.95982,21.11700/ 

DATA XMN(4,4),XMN(4,5),XMN(4,6)/16.22347,19.40942,22.58273/ 

DATA XMN(5,4),XMN(5,5),XMN(5,6)/17.61597,20.82693,24.01902/ 

DATA XMN(6,4),XMN(6,5),XMN(6,6)/18.98013,22.21780,25.43034/ 

DATA XMN(7,4), XMN(7,5), XMN(7,6)/20.32079,23.58608,26.82015/ 

DATA XMN(8,4), XMN (8,5), XMN (8.6 )/21.64154,24.93493,28.19119/ 

DATA XMN(9,4),XMN(9,5),XN(9,6)/22.94517,26.26681,29.54566/ 

DATA CJ,BIG,SMALL/(O.ODO, l.ODO), l.OD+38,1.OD-38/ 

DATA BMAX, MAX, LAYM/200.0,40,600/ 

PI=Z.O*DASIN(l. OD01 

Y0=1.0/(120.O*PI) 

DTOR=P I/ 180.0 

RTOD= 180.0/P I 

TWOPI=2. O*PI 

RETURN 

END 

A-6 
SUBROUTINE INPUT 

c....................................................................... 

c: Reads input data 

c......................................................................: 

IMPLICIT REAL*8 (A, B, D-H, 0-2) 

IMPLICIT COMPLEX*16 (C) 

LOG I CAL LPOL, LS IZE, LSAME, LCON , LPLAZ 

COMMON DIA, MMM, NHOLE, NMODE, NUMTE, NUMTM, F, FCON, ER, LPOL, LSIZE, LSAME 

COMMON /LAYER/ CEP, CUP, NLAY, D(6011, CEC6011, CU(6011, CEUP, LCON, LPLAZ 

COMMON /CONST/ CJ, PI, TWOPI, DTOR, RTOD, BIG, SMALL, BMAX, YO, MAX, LAYM 

COMMON /ARRYS/ MI J (20 1, NI J (20 1, MI JP(20 1, NI JP(20 1, 

A 

B 

CY (40,40 1, CA (40.40 1, CB (40,40 1, IPIV (40 1, INDX (40,2 1, 

AI J (40 I, XI J (40 1, YI J (40 1, PHI JP(40 1 

101 FORMAT( 1x1 

102 

c*****************************************ff***************************** 

C NHOLE = Number of apertures in array. * 

C NMODE = Number of modes per aperture. * 

C NUMTE = Number of TE-modes per aperture. * 

C NUMTM = Number of TM-modes per aperture. * 

C NOTE: (NUMTE+NUMTM)=NMODE * 

C The same modes are assumed in all apertures. * 

C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* 

FORMAT ( lX, 79 ( lH+ 1 1 c*ff****ff***************** INPUT DATA ................................. 

READ ( 5, * , END=9999 1 NHOLE, NMODE, NUMTE, NUMTM 

MMM=NHOLE*NMODE 

IF(MMM.LE.MAX1 GO TO 120 

WRITE(6,llO) 

STOP 

FORMAT(1X’MUTUAL COUPLING OF CIRCULAR APERTURES’) 

IF(NUMTE.GT.20) GO TO 500 

IF(NUMTM. GT. 20) GO TO 500 

IF(NUMTE.GT.NMODE) GO TO 500 

IF(NUMTM.GT.NMODE) GO TO 500 

IF(NMODE.GT.40) GO TO 500 

IF(NUMTE.EQ.0) GO TO 150 

110 FORMAT(lX’(NHOLE*NMODE) EXCEEDS DIMENSION OF CY’) 

120 WRITE(6,130) 

130 

IF((NUMTE+NUMTM).NE.NMODE) GO TO 500 

............................................................................ 

C MIJ(I),NIJ(I) = M,N indices of I-th TE mode. * 

C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* 

140 FORMAT(20(211,1X)) 

150 IF(NUMTM.EQ.0) GO TO 170 ........................................................................ 

C MIJP(I),NIJP(I) = M,N indices of I-th TM mode. # 

C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* 

160 FORMAT(20(211,1X)) 

READ (5,140 1 ( (MI J ( I 1, NI J ( I 1 I, 1=1 , NUMTE) 

READ(5,160) ( (MI JP( I), NI JP ( I) 1, I=l, NUMTM) 

A-7 
170 CONTINUE 

C AIJ(1) = Radius of I-th aperture. 

C XIJ(1) = X-coordinate of center of I-th aperutre. 

C YIJ(1) = Y-coordinate of center of I-th aperture. 

C PHIJP(1) = Angular rotation of XI-axis WRT X-axis (degrees CCW). * 

c******************************ff**************************************** 

* 

* 

* 

c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* 

C****************t***********************************************ff******** 

READ ( 5, * 1 ( (AI J ( I 1 , XI J ( I 1 , Y I J ( I 1, PHI JP ( I 1 1, I =1, NHOLE 

* 

* C F = Frequency (Hertz). 

C FCON = Factor for converting input dimensions to centimeters. 

C (2.54 for input in inches) (1.0 for input in centimeters).* 

C ER = Relative dielectric constant of material filling waveguides 

c...................................................................... 

C 

C F=3. OElO 

FCON=l . 0 * C 

........................................................................ READ(5, *IF, FCON, ER 

C****************************************************************** 

C CEP = Complex (relative) Epsilon of half-space outside of layers.* 

C CUP = Complex (relative) Mu of half-space outside of layers. 

C 

Input format is: (A,-B) (C,-D) * C 

C For free-space outside: (1.0,O. 0) (1.0,O. 0) C For a metal sheet approximation: (-l.OE10,O. 0) (1.0,O. 0) * 

C or CABS(CEP*CUP) greater than 1.OE09 * 

c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* 

C*********************************************************************** 

NLAY = Number of homogeneous layers over apertures. * C 

....................................................................... 

C NLAY=O (half -space 1. 

c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* 

* 

* 

* For input dimensions in wavelengths set: 

* 

* 

* 

READ(5, *)CEP,CUP 

* 

* 

READ ( 5, * 1 NLAY 

IF(NLAY. GT. LAYMISTOP 4444 

IF(NLAY.EQ.0) GO TO 210 . 

C********************ff****************#********************************* 

C NOTE: If NLAY=O, do not inter the following data. 

C 

C 

C CE(1) = Complex Epsilon of I-th layer (relative to free space). * 

C 

* 

* 

D(I) = Thickness of I-th layer. * 

C CU(1) = Complex Mu of I-th layer (relative to free space). * 

C NOTE: Values of CE(1) must be input with format: (A,-B) * 

C Values of CU(1) must be input with format: (C,-D) * 

C where dielectric loss tangent = B/A * 

* 

C * 

A-8 
C CAUT I ON : * 

C An input value of 0.0 for B could produce erroneous results. 

C A value of 0.0001 for the dielectric loss tangent should yield * 

C results which approximate a lossless dielectric in most cases. * 

C (Suggest varying B to determine sensitivity to small values) 

C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++* 

* 

* 

180 

190 

200 

210 

220 

230 

240 

250 

260 

270 

280 

290 

DO 180 I=l,NLAY 

READ (5, * ID ( I 1. CE ( I 1, CU ( I 1 

CONTINUE 

WRITE(6,190)NLAY 

FORMAT(1X’RADIATING INTO MULTI-LAYERS’/lX’NUMBER OF LAYERS = ’141) 

WRITE(6,200) 

FORMAT( 1X’ LAYER’ 3X’ THICKNESS’ 5X’ COMPLEX EPSILON’ 1OX’ COMPLEX MU’ 

IF(DABS(ER).LT. l.OD-04)ER=l.O 

IF(DABS(FCON).LT. l.OD-O4)FCON=l. 0 

XLAM= ( 30. IFCON I/ (F* 1.OD-09 1 

IF(DABS(XLAM-l.O).LT.O.OOOl~XLAM~l.O 

OMEGA=TWOPI*F 

IF(NLAY.EQ.0) GO TO 280 

DO 270 I=l,NLAY 

WR I TE ( 6,220 ) I , D ( I ) , CE ( I ) , CU ( I ) 

FORMAT(lXI4,2XF9.4,3X’ (’F10.5’’ F9.5’ )’3X’ (’F6.2’’F6.2’ 1’ 1 

D(I)=D(I)/XLAM 

EPR=DREAL(CE(I)) 

EPI=DIMAG(CE(I)) 

IF(EP1. LE. -O.OOOl*EPRIGO TO 270 

IF(EP1. LT. 0. OIGO TO 240 

EP I=-EP I 

WRITE(6,230) I 

FORMAT(1X’NOTE: The imaginary part of the dielectric constant’/ 

A8X’ for the layered region must be negative. ’/ 

BSX’This has been corrected for layer ’12/) 

IF(EP1. LE. -O.OOOl*EPR)GO TO 260 

TANL=DABS(EPI/EPR) 

WRITE(6,250)I,TANL 

FORMAT(1X’CAUTION: The dielectric loss tangent for the layered’/ 

AllX’region must not be zero. Erroneous results could also be’/ 

BllX’obtained for very small values of dielectric loss.’/ 

C11X’The loss tangent for layer ’12’ is ’E11.4/) 

CE(I)=DCMPLX(EPR,EPI) 

CONTINUE 

WRITE(6,lOl) 

CONT I NUE 

LCON= . FALSE. 

CEUP=CEP*CUP 

IF(CDABS(CEUP). GE. l.OD+9) LCON=. TRUE. 

ERC=DREAL(CEP) 

EIC=DIMAG(CEP) 

NS IGN=1 

IF(ERC. LT. O.O)NSIGN=-l 

A-9 
300 

310 

320 

330 

340 

350 

360 

370 

380 

390 

400 

410 

420 

430 

ERCP=NSIGN*ERC 

IF(ERCP. GT. l.OD+lO)ERCP=l. OD+10 

ERC=NSIGN*ERCP 

CEP=DCMPLX ( ERC , E I C 1 

IF(NLAY. EQ. 0. AND. (.NOT. LCON). AND. (.NOT. LPLAZ) )WRITE(6,300) 

FORMAT(1X’RADIATING INTO HALF-SPACE’/) 

IF (LCON )WRITE (6,320 1 

CONTINUE 

FORMAT(lX’**** PERFECT CONDUCTOR OUTSIDE  ’/I 

IF(. NOT. LCON)WRITE(6,330)CEP, CUP 

FORMAT(1X’OUTSIDE HALF-SPACE EPSILON = ’2F18.5/ 

IF(FC0N. EQ. 2.54)WRITE(6,340) 

FORMAT( 1X’ ** ** INPUT DIMENSIONS IN INCHES  ’ 

IF( (FCON. EQ. 1.0). AND. (F. EQ. 3.OE10) )WRITE(6,350) 

FORMAT(lX’** ** INPUT DIMENSIONS IN WAVELENGTHS  ’ 

IF( (FCON. EQ. 1.0). AND. (F. NE. 3.OE10) )WRITE(6,360) 

FORMAT( 1X’ ** ** INPUT DIMENSIONS IN CENTIMETERS  ’ 

IF( (FCON. NE. 1.0). AND. (FCON. NE. 2.54) IWRITE(6.370)FCON 

FORMAT( 1X’ ** ** INPUT CONVERSION FACTOR = ’ F10.5’ ****’ 1 

IF( (FCON. NE. 1.0). OR. (F. NE. 3.OE10) )WRITE(6,380)F 

WRITE(6,390) 

FORMAT(/lX’++++++++ APERTURE ARRAY GEOMETRY +++++++++’/ 

LPOL=. TRUE. 

LS I ZE= . TRUE. 

DI A=2. O*AI J ( 1 1 

DO 410 I=l,NHOLE 

WRITE(6.400) I, AIJ( I), XIJ( I 1, YIJ( I), PHIJP( I) 

FORMAT( lX, 13, lX, F9.4, lX, F8.3, lX, F8.3, lX, F8.3) 

AIJ(I)=AIJ(I)/XLAM 

XIJ(I)=XIJ(I)/XLAM 

YIJ(I)=YIJ(I)/XLAM 

PHIJP(I)=PHIJP(I)*DTOR 

IF(DABS(AIJ(I)-AIJ(1)).GT.O.OOOl)LSIZE=.FALSE. 

IF (DABS (PHI JP ( I ) -PHI JP ( 1 

CONTINUE 

LSAME=. FALSE. 

IF(LS1ZE. AND. LPOL)LSAME=. TRUE. 

WRITE(6,lOl) 

IF(NUMTE.EQ.0) GO TO 440 

DO 430 I=l,NUMTE 

WRITE(6,420) I MI J ( I), NIJ ( I) 

FORMAT(1X’MODE’ 12’ = TE-’ 11’11) 

IF((l+MIJ(I)).GT.9) GO TO 480 

IF(NIJ(I).GT.6) GO TO 480 

IF(MIJ(I).LT.O) GO TO 480 

IF(NIJ(I).LT.l) GO TO 480 

CONT I NUE 

+ 1X’OUTSIDE HALF-SPACE MU = ’2F18.5/) 

FORMAT(lX’**** FREQUENCY = ’Ell. 6’ HERTZ ****’ 1 

A1X’ HOLE’ ’3X’RADIUS’ ,6X, ’ X’ ,8X, ’ Y’ ,7X, ’ POL. ’ ) 

1. GT. 0.0001 LPOL=. FALSE. 

A-10 
440 IF(NUMTM.EQ.0) GO TO 470 

DO 460 I=l,NUMTM 

I P= I +NUMTE 

WRITE (6,450 ) IP, MI JP ( I 1, NI JP ( I 

450 FORMAT(lX’MODE’12’ = TM-’11,111 

IF((l+MIJP(I)).GT.9) GO TO 480 

IF(NIJP(I).GT.6) GO TO 480 

IF(MIJP(II.LT.0) GO TO 480 

IF(NIJP(I).LT. 1) GO TO 480 

460 CONTINUE 

470 WRITE(6,lOl) 

RETURN 

480 WRITE(6,490) 

490 FORMAT(1X’RANGE OF MODE INDICES EXCEEDED’/ 

A1X’ (M, N) VALUES ALLOWED: (0,l) THROUGH (8.6)’ ) 

STOP 

FORMAT(1X’NUMBER OF MODES EXCEEDS ARRAY DIMENSIONS’// 

STOP 

STOP 

END 

500 WRITE(6,51O)NHOLE, NMODE, NUMTE, NUMTM 

510 

A1X’ “OLE=’ 13,5X’ NMODE=’ 13,5X’ NUMTE=’ 13,5X’ NUMTM=’ 13) 

9999 WRITE(6,102) 

A-1 1 
SUBROUTINE YMAT ........................................................................ 

c: Calculates coefficients CY(I1,JJ) of complex admittance matrix [Y] : 

........................................................................ 

IMPLICIT REAL*8 (A, B, D-H, 0-2) 

IMPLICIT COMPLEX*16 (C) 

LOG I CAL TE I, TEJ , TM I, TMJ , LPOL , LS I ZE, LSAME, LCON , LPLAZ 

DIMENSION Y (21, YA(2) 

COMMON D I A, MMM, NHOLE , NMODE , NUMTE , NUMTM, F, FCON , ER, LPOL, LS I ZE , LSAME 

COMMON /APERS/ IT, MI, NI, MJ, NJ, AI, AJ, TEI, TEJ, TMI, TMJ, PHI, PHIP, R 

COMMON /CONST/ C J , PI, TWOP I, DTOR, RTOD, BIG, SMALL, BMAX , YO, MAX, LAYM 

COMMON /LAYER/ CEP, CUP, NLAY, D(6011, CE(6011, CU(6011, CEUP, LCON, LPLAZ 

COMMON /ARRYS/ MI J (20 I, NI J (20 I, MI JP (20 1, N I JP (20 1, 

A CY(40,4O),CA(40,40),CB(40,40), IPIV(40), INDX(40,2), 

B AI J (40), XI J (401, Y I J (401, PHI JP (40) 

EXTERNAL FCTCWG 

DO 1000 IHOLE=l, "OLE 

DO 1000 JHOLE=l , "OLE 

DO 1000 IMODE=l, NMODE 

DO 1000 JMODE=l , NMODE 

101 FORMAT(1X) 

II=(IHOLE-1)*NMODE+IMODE 

JJ=(JHOLE-l)*NMODE+JMODE 

IF(II.LE.JJ) GO TO 100 

CY(I1, JJ)=CY(JJ, 11) 

GO TO 1000 

IF(IHOLE.EQ. 1) GO TO 110 

IF( IHOLE. NE. JHOLE) GO TO 110 

IF(.NOT.LSIZE) GO TO 110 

CY ( I I, J J 1 =CY ( IMODE, JMODE 1 

GO TO 1000 

110 TEI=. TRUE. 

TEJ=. TRUE. 

TMI=. FALSE. 

TMJ=. FALSE. 

IF(IMODE.LE.NUMTE) GO TO 120 

TEI=. FALSE. 

TMI=. TRUE. 

TE J=. FALSE. 

TMJ=. TRUE. 

I T=O 

IF(TEI.AND. TEJ)IT=l 

IF(TMI.AND. TMJIIT=2 

IF(TEI.AND.TMJIIT=3 

IF(TM1. AND. TEJ 1 IT=4 

AI=AIJ(IHOLE) 

AJ=AIJ(JHOLE) 

100 CONTINUE 

120 IF(JM0DE.LE.NUMTE) GO TO 130 

130 CONTINUE 

A-12 
140 

150 

160 

170 

180 

190 

200 

210 

PHPI=PHIJP(IHOLE) 

PHPJ=PHIJP(JHOLE) 

IF(TMI)PHPI=PHPI-0.5*PI 

IF(TMJ)PHPJ=PHPJ-0.5*PI 

PHIPzPHPJ-PHPI 

XJI=XIJ(JHOLE)-XIJ(IH0LE) 

YJI=YIJ(JHOLE)-YIJ(IH0LE) 

R=DSQRT(XJI*XJI+YJI*YJI) 

IF( IHOLE. EQ. JHOLE)R=O. 0 

IF(DABS(R).LT. l.OD-20)R=O.O 

IF(DABS(PH1P). LT. l.OD-OI)PHIP=O. 0 

PHI=O. 0 

IF(R.LT. 1.OD-04) GO TO 160 

IF(DABS(XJ1). LT. 1.OD-04) GO TO 140 

PHI=DATANZ(YJI,XJI)-PHPI 

GO TO 150 

PHI=O. 5*PI 

IF(YJ1. LT. O.O)PHI=PHI+PI 

PHIzPHI-PHPI 

CONTINUE 

IF(TM1) GO TO 170 

IF(.NOT.TEI) STOP 1111 

MI=MIJ(IMODE) 

NI=NIJ(IMODE) 

GO TO 180 

IDEM=IMODE-NUMTE 

MI=MIJP(IDEM) 

NI=NIJP(IDEM) 

IF(TMJl GO TO 190 

IF(. NOT. TEJ 1 STOP 2222 

MJ=MIJ(JMODE) 

NJ=NIJ(JMODE) 

GO TO 200 

JDEMzJMODE-NUMTE 

MJ=MI JP (JDEM) 

NJ=NIJP(JDEM) 

CONT I NUE 

IF(IH0LE. NE. JHOLE) GO TO 210 

IF(MJ.EQ.MIl GO TO 210 

CY(I1, JJ)=(O.O,O.O) 

GO TO 1000 

CONTINUE 

IF(IH0LE. EQ. JHOLE) GO TO 230 

IF(IHOLE.EQ. 1) GO TO 230 

IF(.NOT.LSIZE) GO TO 230 

IF(. NOT. LPOL) GO TO 230 

IHMl=IHOLE-l 

JHMl=JHOLE-l 

DO 220 IKH=l, IHMl 

DO 220 JKH=1, JHMl 

A-13 
220 

230 

240 

250 

260 

270 

280 

IF (DABS (AI -A I J ( IKH ) ) . GT. 0.00 1 ) GO TO 220 

IF(DABS(AJ-AIJ(JKH)).GT.O.OOl) GO TO 220 

TX=DABS(DABS(XJI)-DABS(XIJ(1KH)-XIJ(JKH))) 

IF(TX.GT.0.001) GO TO 220 

TY=DABS(DABS(YJI)-DABS(YIJ(IKH)-YIJ(JKH))) 

IF(TY.GT.0.001) GO TO 220 

PT=DABS(DABS(PHIP)-DABS(PHIJP(IKH)-PHIJP(JKH))) 

IF(PT.GT.0.001) GO TO 220 

IK=(IKH-l)*NMODE+IMODE 

JK=(JKH-l)*NMODE+JMODE 

CY ( I I, J J 1 =CY ( IK, JK 1 

GO TO 1000 

CONTINUE 

CONT I NUE 

EMI=2.0 

EMJ=2.0 

IF(M1. EQ. O)EMI=l. 0 

IF(MJ. EQ. O)EMJ=l. 0 

CALL FACS(MI,NI,MJ,NJ,AI,AJ,R,PHI,PHIP,EMI) 

YTEST=O. 0001 

Y ( 1 )=O. 0 

Y(2)=0.0 

NDEL=50 

BDEL=O. 25 

AA=O. 0001 

B B = B D E L 

TOL=O. 0001 

CALL INTG2 ( AA, BB, FCTCWG, NDEL, TOL, YA, Y ) 

Y(l)=Y(l)+YA(l) 

Y (2 1 =Y (2 1 +YA (2 1 

IF(NLAY.EQ.0) GO TO 270 

DO 250 N=l,NLAY 

EUREAL=DREAL(CE(N)*CU(N)) 

IF ( EUREAL . LT . 0 . 0 EUREAL=-EUREAL 

IF( (BB*BB). LE. EUREAL) GO TO 260 

CONT I NUE 

GO TO 270 

CONT I NUE 

IF(DABS(BB-1.0). LT. 0.01 )BB=l. 0 

AA=BB 

BB=BB+BDEL 

IF(DABS(AA-l.O).LT.O.Ol~AA=l.OOOl 

IF(DABS(BB-l.O),LT.O. 01)BB=0.9999 

GO TO 240 

IF(BB.LT.2.0) GO TO 260 

BDEL=l . 0 

NBDEL=O 

AA=BB 

BB=BB+BDEL 

CALL I NTG2 ( AA, BB , FCTCWG, NDEL , TOL, YA , Y 1 

A-14 
290 

300 

310 

320 

330 

340 

Y(l)=Y(l)+YA(l) 

Y(Z)=Y(Z)+YA(Z) 

IF(DABS(Y(l)).LE.YTEST) GO TO 290 

IF(DABS(YA(l)/Y(l)).GT.YTEST) GO TO 310 

IF(DABS(Y(Z)).LE.YTESTl GO TO 300 

IF(DABS(YA(Z)/Y(Z)).GT.YTEST) GO TO 310 

NBDEL=NBDEL+l 

IF(NBDEL.GE.5) GO TO 330 

GO TO 320 

NBDEL=O 

IF(BB.GE.BMAX) GO TO 330 

IF(BB. GE. ZO.O)BDEL=S. 0 

AA=BB 

BB=BB+BDEL 

GO TO 280 

CYIJ=-YO*DSQRT(EMI*EMJ)*DCMPLX(Y(l),Y(2)) 

CY ( I I, JJ )=CY I J 

WRITE(6,340)II, JJ,CY(II, JJ),BB, IHOLE, IMODE, JHOLE, JMODE 

FORMAT(lX’CY(’IZ’, ’12’ )=(’E11.4’,’E11.4’ )’3X’BB=’F6.2. 

A4X’ IH=’ I2,2X’ IM=’ 12,3X’JH=’ I2,2X’JM=’ 12) 

1000 CONTINUE 

WRITE(6,lOl) 

RETURN 

END 

A-15 
SUBROUTINE SMAT ........................................................................ 

c: Calculates coefficients of complex scattering matrix [SI (see eq.1): 

........................................................................ 

IMPLICIT REAL*8 (A, B, D-H, 0-2) 

IMPLICIT COMPLEX*16 (C) 

LOGICAL TEI, TEJ, TMI, TMJ 

LOG I CAL LPOL, LS I ZE , LSAME 

COMMON DIA, MMM, NHOLE, NMODE, NUMTE, NUMTM, F, FCON, EFt, LPOL, LSIZE, LSAME 

COMMON /APEX/ IT,MI,NI,MJ,NJ,AI,AJ,TEI,TEJ,TMI,TMJ,PHI,PHIP,R 

COMMON /ZEROS/ XMNP(9,6),XMN(9,6) 

COMMON /CONST/ C J , PI, TWOP I, DTOR, RTOD , BIG, SMALL, BMAX, Y 0, MAX, LAYM 

COMMON /ARRYS/ MI J (201, NI J (201, MI JP(ZO), NI JP(ZO), 

A CY(40,40),CA(40,40),CB(40,40), IPIV(401, INDX(40,2), 

B AI J (40 1, XI J (40 1, Y I J (40 1, PHI JP (40 1 

DO 130 IHOLE=l , NHOLE 

DO 130 JHOLE=l , NHOLE 

DO 130 IMODE=l,NMODE 

DO 130 JMODE= 1. NMODE 

TMI=. FALSE. 

IF( IMODE. GT. NUMTE)TMI=. TRUE. 

II=(IHOLE-1)*NMODE+IMODE 

JJ=(JHOLE-l)*NMODE+JMODE 

AKOI=AIJ(IHOLE)*TWOPI 

IF(TM1)GO TO 100 

NI=NIJ(IMODE) 

MIP=MIJ(IMODE)+l 

FSQI=(XMNP(MIP,NI)/AKOI)**Z 

FTSQ=ER-FSQI 

FCON-DABS(FTSQ1 

IF (FCON. LT. SMALL 1 FCON=O . 0 

IF(FC0N. GT. BIG)FCON=BIG 

IF (FTSQ. GE. 0.0 1 CYMNO=YO*DSQRT ( FCON 1 

I F ( FTSQ . LT . 0 . 0 CYMNO=-C J *YO*DSQRT ( FCON 

GO TO 110 

NI=NIJP(IDEM) 

MIP=MIJP(IDEM)+l 

FSQI=(XMN(MIP,NI)/AKOI)**Z 

FTSQzER-FSQI 

FCON=DABS(FTSQ) 

IF (FCON . LT. SMALL 1 FCON=SMALL 

IF(FC0N. GT. BIG)FCON=BIG 

IF ( FTSQ . GE. 0.0 ) CYMNO=YO/DSQRT ( FCON ) 

IF(FTSQ. LT. O.O)CYMNO=YO/(-CJ*DSQRT(FCON) ) 

100 IDEMaIMODE-NUMTE 

110 CONTINUE 

CA(I1, JJ)=CY(II, JJ) 

CB(I1, JJ)=-CY(II, JJ) 

IF( I I. EQ. JJ )CA ( I I, JJ )=CYMNO+CY ( I I, JJ 1 

IF ( I I. EQ. J J CB ( I I, J J ) =CYMNO-CY ( I I, J J ) 

A-16 
120 

130 

140 

150 

160 

170 

180 

190 

200 

CY(I1, JJ)=CB(II, JJ) 

IF ( I I. EQ. JJ )WRITE (6,120 1 I I, CYMNO, IHOLE, IMODE 

FORMAT(lX’YMN(’I2’ )=(’E11.4’,’E11.4’ )’3X’IHOLE=’I2,3X’IMODE=’I2) 

CONT I NUE 

WRITE(6,140) 

FORMAT ( / 1 OX’ SCATTER I NG MATRIX ’ 1 

CALL CINV(CA, MMM, IPIV, INDX, MAX) 

DO 160 I=l,MMM 

DO 160 J=l,MMM 

CB(1, J)=(O. 0,O. 0) 

DO 150 K=l,MMM 

CB( I, J )=CB( I, J )+CY ( I, K)*CA(K, J 1 

CONT I NUE 

CONTINUE 

DO 200 I=l,MMM 

DO 200 J=l,MMM 

xIsoL=-3oo.o 

XDX=CDABS ( CB ( I, J 1 

IF(XDX. LT. 1.OD-15)GO TO 180 

XISOL=20.*DLOGlO(XDX) 

XXl=DIMAG(CB(I, J)) 

XX2=DREAL(CB(I, J)) 

PHASE=RTOD*DATAN2(XXl,XXZ) 

WRITE (6,170 1 I, J, CB ( I, J 1 , XISOL, PHASE 

FORMAT ( 1X’ S ( ’ I2 ’ , ’ I2 ’ 1 = ( ’ El 1.4’ , ’ El 1.4’ ) ’ 3XF9.4’ DB’ 3XF9.4’ DEG’ 

GO TO 200 

WRITE(6,lgO 1 I, J, CB( I, J 1 

FORMAT(lX’S(’ 12’, ’ 12’ )=(’Ell. 4’, ’Ell. 4’ 1’3X’BELOW -300 DB’ 

CONTINUE 

RETURN 

END 

A-17 
SUBROUTINE CINV(CA, N, IPIV, INDX, MAX) 

c............................................................... ........ c: Performs matrix inversion of (N by N) complex matrix CA. 

c: On return, CA contains inverted matrix. CA must be dimensioned at : 

c: least N by N in calling routine. IPIV and INDX are used internally: 

c: by CINV and must be dimensioned at least N and N by 2 in calling : 

c: routine. MAX is row dimension of CA and INDX in calling routine. : 

c......................................................................: 

100 

110 

120 

130 

140 

150 

IMPLICIT REAL*8 (A, B, D-H, 0-Z) 

IMPLICIT COMPLEX*16 (C) 

DIMENSION CA(MAX,N), IPIV(N), INDX(MAX, 2) 

DATA C0,Cl /(O.ODO,O.ODO), (l.ODO,O.ODO)/ 

DATA BIG, SMALL /l. OD+38,1.OD-38/ 

ISCALE=O 

RL=B I G 

RS=SMALL 

CDET=Cl 

ADM=1 . 0 

DO 100 J=l,N 

IPIV (J 1 =O 

DO 220 I=l,N 

AVM=O. 0 

DO 130 J=l,N 

IF(IPIV(J).EQ.l) GO TO 130 

DO 120 K=l,N 

IF(IPIV(K)-l)llO,l20,260 

CONT I NUE 

AVA=CDABS(CA(J,K)) 

IF(AVM.GE.AVA1 GO TO 120 

IROW=J 

I COL=K 

AVM=AVA 

CONTINUE 

CONT I NUE 

IF(AVM.EQ.O.0) GO TO 250 

IPIv(IcoL)=IPIv(IcoL~+l 

IF(IROW.EQ. ICOL) GO TO 150 

CDET=-CDET 

DO 140 L=l,N 

CSWAP=CA(IROW,L) 

CA(IROW, L)=CA( ICOL, L) 

CA( ICOL, L)=CSWAP 

CONTINUE 

CONT I NUE 

INDX( I, 1 )=IROW 

INDX(I,Z)=ICOL 

CPIV=CA( ICOL, ICOL) 

APV=CDABS(CPIV) 

IF(APV.EQ.O.0) GO TO 250 

CPIVI=CPIV 

A-18 
ADM=CDABS(CDET) 

IF(ADM.LT.RL1 GO TO 160 

CDET=CDET/RL 

ADM=CDABS(CDET) 

ISCALE=ISCALE+l 

IF(ADM.LT.RL) GO TO 170 

CDET=CDET/RL 

ISCALE=ISCALE+l 

GO TO 170 

IF(ADM.GT.RS1 GO TO 170 

CDET=CDET*RL 

ADM=CDABS(CDET) 

ISCALE=ISCALE-1 

IF(ADM.GT.RS) GO TO 170 

CDET=CDET*RL 

ISCALE=ISCALE-1 

APV=CDABS(CPIVI) 

IF(APV.LT.RL) GO TO 180 

CPIVI=CPIVI/RL 

APV=CDABS(CPIVI) 

ISCALE=ISCALE+l 

IF(APV.LT.RL) GO TO 190 

CPIVI=CPIVI/RL 

ISCALE=ISCALE+l 

GO TO 190 

IF(APV.GT.RS1 GO TO 190 

CPIVI=CPIVI*RL 

APV=CDABS(CPIVI) 

IF(APV.GT.RS1 GO TO 190 

CPIVI=CPIVI*RL 

ISCALE-ISCALE-1 

CDET=CDET*CPIVI 

CA( ICOL, ICOL)=Cl 

DO 200 L=l,N 

DO 220 L1=1,N 

IF(L1.EQ. ICOL) GO TO 220 

CSWAP=CA(Ll,ICOL) 

CA(L1, ICOL)=CO 

DO 210 L=l,N 

160 CONTINUE 

170 CONTINUE 

180 CONTINUE 

ISCALE=ISCALE-1 

190 CONTINUE 

200 CA(ICOL,L)=CA(ICOL,L)/CPIV 

210 CA(Ll,L)=CA(Ll,L)-CA(ICOL,L)*CSWAP 

220 CONTINUE 

DO 240 I=l,N 

L=N+ 1 -I 

IF(INDX(L, 1).EQ. INDX(L,2)) GO TO 240 

A-19 
IROW=INDX(L, 1) 

ICOL=INDX(L, 2) 

DO 230 K=l,N 

CSWAP=CA (K, IROW 1 

CA(K,IROW)=CA(K,ICOL) 

CA (K, ICOL 1 =CSWAP 

230 CONTINUE 

240 CONTINUE 

250 CDET=CO 

260 RETURN 

GO TO 260 

I SCALE=O 

END 

A-20 
SUBROUTINE INTG2 (AA, BB, FX, NT, TOL, SUM, Y 1 ................................................. ....................... c: Integrates 2 real functions (may be real and imaginary parts of : 

c: of a complex function) over limits of AA to BB by using a Simpson : 

c: integration procedure with interval halving until relative error : 

c: TOL is reached on each interval. NT is number of intervals. 

c: SUM is a 2 element array of the resultant integrals over AA to BB. : 

c: Y is the cumulative values of SUM from previous calls to INTG2. 

c: Y must be set to zero prior to first call to INTG2 

c: and is used when integrating FX over wide limits ............... c: by successive calls to INTGZ. FX(X,WK) is user : May 1989 : 

c: supplied subroutine to evaluate integrands WK. : M.C.Bailey : 

c: FX must be declared EXTERNAL in calling routine. : Hampton, VA : 

c:..... .................................................. :.............: 

100 

110 

200 

IMPLICIT REAL*8 (A-H,O-Z) 

DIMENSION SUM (2 1, Y (2 1, WK (2 1, SUMA (2 1, SUMB (2 1, ESUM (2 9 BASE (2 

DIMENSION A(401,2), B(401,Z) 

DATA NMAX, N IMAX/400,2/ 

IF(NT. GT. NMAX) STOP 

DEL=(BB-AA)/NT 

DEL02=O. 5DO*DEL 

DO 100 I=l,NT+l 

X=AA+(I-l)*DEL 

CALL FX(X,WK) 

A(I, l)=WK(l) 

A ( I, 2 1 =WK (2 1 

CONTINUE 

BASE(l)-DELO2*(A(l,l)+A(NT+l,l)) 

BASE(2)=DELOZ*(A(1,2)+A(NT+l,2)) 

DO 110 I=Z,NT 

BASE( 1 )=BASE( 1 )+DEL*A( I, 1 

BASE(2)=BASE(2)+DEL*A(I,2) 

CONT I NUE 

DENOM=l . D-28 

TEST=DABS(BASE(l)+Y(l)) 

IF ( TEST. GT . DENOM DENOM=TEST 

TEST=DABS(BASE(2)+Y(Z)) 

IF ( TEST. GT . DENOM 1 DENOM=TEST 

SUM(l)=O.ODO 

SUM(2)-0. OD0 

DO 300 K=l,NT 

ND= 1 

FF=AA+(K-l)*DEL 

X=FF+DEL 

CALL FX(X,WK) 

A(2,l )=WK(l) 

A(2,2)=WK(2) 

SUMA ( 1 )=DEL02* (A( 1,l) +A(2,1) 1 

SUMA(2)=DEL02*(A(1,2)+A(2,2) 1 

ND=ND+ND 

A-21 
IF(ND.GT.NMAX) GO TO 270 

NB=ND/2 

NA=NB+ 1 

NP=ND+ 1 

DELN=DEL/ND 

IF(DELN.LT. l.D-6) GO TO 270 

M) 210 J=l,NB 

X=FF+DELN*(J+J-l) 

CALL FX(X,WK) 

B(J, l)=WK(l) 

B(J, 2)=WK(2) 

210 CONTINUE 

NAP=NA+ 1 

NBP=NB+ 1 

DO 220 J=1, NA 

JL=NAP-J 

JP=JL+JL-l 

A(JP, l)=A(JL, 1) 

A(JP,Z)=A(JL,Z) 

220 CONTINUE 

DO 230 J=1, NB 

JL=NBP- J 

JP=JL+ JL 

A(JP, l)=B(JL, 1) 

A(JP,Z)=B(JL,Z) 

SUMB( 1 )=O. 5DO* (A( 1,l )+A(NP, 1 1 1 

SUMB(2)=0.5DO*(A(l,2)+A(NP,Z)) 

DO 240 NS=2, ND 

SUMB(l)=SUMB(l)+A(NS,l) 

SUMB (2 1 =SUMB (2 1 +A (NS, 2 1 

SUMB(l)=DELN*SUMB(l) 

SUMB(2)=DELN*SUMB(Z) 

ESUM(l)=DABS(SUMB(l)-SUMA(l))/DENOM 

ESUM(2)=DABS(SvMB(Z)-SUMA(2))/DENOM 

230 CONTINUE 

240 CONTINUE 

SUMA(l)=SUMB(l) 

SUMA(ZI=SUMB(Z) 

IF(ESUM(l).GT.TOL) GO TO 200 

IF(ESUM(2). GT. TOL) GO TO 200 

GO TO 280 

270 ND=ND/2 

280 CONTINUE 

SUM(l)=SUM(l)+SUMA(l) 

SUM(Z)=SUM(2)+SUMA(Z) 

A(l,l)=A(NP, 1) 

A(l,Z)=A(NP, 2) 

300 CONTINUE 

RETURN 

END 

A-22 
SUBROUTINE FCTCWG(BETA,FI) 

c... .................................................................... c: Evaluates integrand (real and imaginary parts) of equation 2. 

c:................... .................................................. : 

10 

15 

20 

25 

30 

40 

50 

900 

90 1 

IMPLICIT REAL*8 (A, B, D-H, 0-2) 

IMPLICIT COMPLEX*16 (C) 

LOGICAL TEI, TEJ, TMI, TMJ, LCON, LPLAZ 

DIMENSION FI(2) 

COMMON /APERS/ IT, MI, N I, MJ, NJ, AI, AJ, TEI , TEJ, TMI , TMJ , PHI, PHIP, R 

COMMON /LAYER/ CEP,CUP,NLAY,D(601~,CE(6Ol~,CU(6Ol),CEUP,LCON,LPLAZ 

CALL MLAYER (BETA, CW1, CW2 ) 

CALL UVIJ(BETA,UIJ,VIJ, IT) 

GO TO (10,20,30,40) IT 

GO TO 900 

CALL STEI (BETA, SSI, TTI 1 

IF(M1. NE. MJ) GO TO 15 

IF(NI.NE.NJ) GO TO 15 

IF(DABS(A1-AJ). GT. 1. D-4) GO TO 15 

SS J=SS I 

TTJ=TTI 

GO TO 50 

CALL STEJ (BETA, SSJ, TTJ 1 

GO TO 50 

CALL STMI (BETA, SSI, TTI 1 

IF(MI.NE.MJ) GO TO 25 

IF(NI.NE.NJ) GO TO 25 

IF(DABS(A1-AJ).GT. l.D-4) GO TO 25 

SSJ=SSI 

TTJ=TT I 

GO TO 50 

CALL STMJ (BETA, SSJ, TTJ 1 

GO TO 50 

CALL STEI (BETA, SSI, TTI 1 

CALL STMJ (BETA, SSJ, TTJ 

GO TO 50 

CALL STMI (BETA, SSI, TTI ) 

CALL STEJ(BETA, SSJ, TTJ) 

Fl=+SSI*SSJ*UIJ 

FZ=-TTI*TTJ*VIJ 

CCW=(Fl*CWl+FZ*CW2)*BETA 

FI(l)=DREAL(CCW) 

FI(Zl=DIMAG(CCW) 

RETURN 

WRITE(6,901) IT 

FORMAT(lX’### ERRORS IN FCTCWG ###’ 

A/lX’(IT must be between 1 and 4)’ 

B/lX’ IT=’ 15) 

STOP 

END 

A-23 
SUBROUT I NE MLAYER ( BETA, CW1, CW2 1 ........................................................................ 

c: Evaluates complex functions in equations 5 and 6. 

........................................................................ 

IMPLICIT REAL*8 (A, B, D-H, 0-Z) 

IMPLICIT COMPLEX*16 (C) 

LOGICAL LCON , LPLAZ 

COMMON /LAYER/ CEP, CUP, NLAY, D(6011, CE(6011, CU(6011, CEUP, LCON, LPLAZ 

DATA XKO, BIG, CSMAL/6.283185307179586DO, l.OD+28, (1.OD-28,O. ODO)/ 

DATA COY C1, CJ/ (0. ODO, 0. ODO) , ( 1. ODO, 0. ODO) , (0. ODO, 1. ODO)/ 

CU(NLAY+l)=CUP 

CE(NLAY+l)=CEP 

BSQ=BETA*BETA 

CPF=CO 

CGPOG=CO 

CFPOF=C 1 

IF(LC0N) GO TO 35 

CPF=Cl 

CBSQPXEUP-BSQ 

IF(LPLAZ)CBSQP=Cl-BSQ 

IF(DREAL(CBSQP))20,10,10 

10 CKZ=XKO*CDSQRT(CBSQP) 

GO TO 30 

20 CK2-CJ*XKO*CDSQRT(-CBSQP) 

30 CFPOFz-CJ*CKZ 

CGPOG=CFPOF 

35 IF(NLAY.EQ.0) GO TO 200 

N=NLAY 

40 CARG=CE(N)*CU(N)-BSQ 

IF(DREAL(CARG1. LT. 0.0) GO TO 50 

CKZ=XKO*CDSQRT(CARG) 

GO TO 60 

50 CKZ=-CJ*XKO*CDSQRT(-CARG) 

60 CKZD=CKZ*D(N) 

A=DREAL(CKZD) 

B=DIMAG(CKZD) 

DSINA=DSIN(A) 

DCOSA=DCOS(A) 

NSIGN=l 

IF(B. LT. O.O)NSIGN=-l 

BP=NSIGN*B 

IF(BP.GT. 1.OD-10) GO TO 70 

CSINKD=DSINA 

CCOSKD=DCOSA 

GO TO 100 

DCOSHB=O. 5* (1.6948892D+28 1 

DSINHB=NSIGN*DCOSHB 

GO TO 90 

70 IF(BP. LT. 65) GO TO 80 

80 EP=DEXP ( BP 1 

A-24 
. EM=l.ODO/EP 

DCOSHB=O.SDO*(EP+EM) 

DSINHB=NSIGN*DSINHB 

DSINHB=O.SDO*(EP-EM) 

90 CSINKD=DSINA*DCOSHB+CJ*DCOSA*DSINHB 

CCOSKD=DCOSA*DCOSHB-CJ*DSINA*DSINHB 

100 CONTINUE 

CcV=CFPoF*(cV(N)/cV(N+l))/CKZ 

CCE=CGWG*(CE(N)/CE(N+l))/CKZ 

CFOFl=CPF*CCOSKD-CCU*CSINKD . 

CGOGl=CCOSKD-CCE*CSINKD 

CFPOFl=CKZ*(CPF*CSINKD+CCU*CCOSKD) 

CGPOGl=CKZ*(CSINKD+CCE*CCOSKD) 

CFPOF=CFPOFl/CFOFl 

CGPOG=CGPOGl/CGOG1 

CPF=C 1 

N-N- 1 

IF(N.GT.0) GO TO 40 

CW2=CJ*CFPOF/(XKO*cV(l)) 

RETURN 

END 

200 CWl=-CJ*XKO*CE(l)/CGPOG 

A-25 
SUBROUTINE FACS(MII,NI,MJJ,NJ,AI,AJ,R,PHI,PHIP,EMI) 

IMPLICIT REAL*8(A-H,O-Z) 

COMMON /FACTS/ MI,MJ,MP,MM,MPP,MMP,MIl,MJl,MXP,MPI,MPJ,MPPI,MPPJ, 

A 

B 

C 

EM1, EM2, AIKO, AJKO , RKO , XEI , XEJ , XM I , XMJ , 

XESI,XESJ,XMSI,XMSJ,SFI,SFJ, TFI, TFJ, TFEI, TFEJ, 

XDS, XDC, XDSP, XDCP, XDSM, XDCM 

COMMON /ZEROS/ XMNP(9,6), XMN(9,6) 

PI=2.ODO*DASIN(l. OD01 

TWOPI=2.ODO*PI 

MI=MI I 

MJ=MJJ 

MP=M J+M I 

MM=M J-M I 

MPP=MP+ 1 

MMP=MM+ 1 

MIl=(-l)**MI 

MJl=(-l)**MJ 

MXP=-MM+ 1 

MP I =MI + 1 

MP J=M J+ 1 

MPPI=MPI+l 

MPP J=MP J+ 1 

EMl=EMI-l. 0 

EM2=2.0/EMI 

AIKO=AI*TWOPI 

AJKO=AJ*TWOPI 

RKO=R*TWOP I 

PHIPMI=MI*PHIP 

PHIPMJ=MJ*PHIP 

XDS=DSIN(PHIPMI) 

XDC=DCOS(PHIPMI) 

AGMP=MP*PHI-PHIPMJ 

AGMM=MM*PHI-PHIPMJ 

XDSP=DSIN(AGMP) 

XDCP=DCOS(AGMP) 

XDSM=DSIN(AGMM) 

XDCM=DCOS(AGMM) 

XEI=XMNP(MPI, NI 1 

XEJ=XMNP(MPJ,NJ) 

XMI=XMN(MPI, NI 1 

XMJ=XMN(MPJ,NJ) 

XESI=XEI*XEI 

XESJ=XEJ*XEJ 

XMS I =xM I *xM I 

XMS J=XM J * XM J 

DEI=XESI-MI*MI 

DEJzXESJ-MJ*MJ 

FEI=AIKO/DSQRT(DEI) 

FEJ=AJKO/DSQRT(DEJ) 

SFI=MI*FEI 

A-26 
SFJ=MJ*FEJ 

TFI=XESI*FEI 

TFJ=XESJ*FEJ 

TFEI=-0.5DO*MI*MI/(XI*DEI) 

TFEJ=-0.5DO*MJ*MJ/(XJ*DEJ) 

RETURN 

END 

A-27 
SUBROUTINE UVI J (BETA, UI J, VI J, IT 1 

c.......................,........................................... .... c: Evaluates functions in table I or table 11. 

........................................................................ .- 

IMPLICIT REAL*8(A-H,O-Z) 

DIMENSION BJ(21) 

COMMON /FACTS/ MI,MJ,MP,MM,MPP,MMP,MIl,MJl,MXP,MPI,MPJ,MPPI,MPPJ, 

A 

B 

C 

EM1, EM2, AIKO, AJKO,RKO, XEI, XEJ, XMI,XMJ, 

XESI, XESJ, XMSI, XMSJ, SFI, SFJ, TFI, TFJ, TFEI, TFEJ, 

XDS , XDC, XDSP , XDCP , XDSM, XDCM 

UI J=O. OD0 

VI J=O. OD0 

IF(RKO.GT. 1.OD-04) GO TO 150 

IF(MI.NE.MJ) RETURN 

GO TO (110,120,130,140)IT 

110 UIJ=-EMl*XDC 

VIJ=EMZ*XDC 

RETURN 

120 UIJ=-EM2*XDC 

RETURN 

130 UIJ=EMl*XDS 

RETURN 

140 UIJ=EMl*XDS 

RETURN 

150 ARG=RKO*BETA 

CALL BESJS(ARG, BJ, MPP) 

GO TO (160,190,220,250) IT

IF(MI.GT.MJ) GO TO 170 

TZ=MIl*BJ(MMP)*XDCM 

GO TO 180 

170 TZ=MJl*BJ(MXP)*XDCM 

180 UIJ=MJl*(Tl-T2) 

VIJ=MJl*(Tl+TZ) 

RETURN 

190 Tl=BJ(MPP)*XDCP 

IF(MI.GT.MJ1 GO TO 200 

TZ=MIl*BJ(MMP)*XDCM 

GO TO 210 

200 TZ=MJl*BJ(MXP)*XDCM 

210 UIJ=-MJl*(Tl+T2) 

220 Tl=BJ(MPP)*XDSP 

160 Tl=BJ(MPP)*XDCP 

RETURN 

IF(M1. GT. MJ) GO TO 230 

T2=MIl*BJ(MMP)*XDSM 

GO TO 240 

230 T2=MJl*BJ(MXP)*XDSM 

240 UIJ=MJl*(Tl-TZ) 

250 Tl=BJ(MPP)*XDSP 

RETURN 

. A-28 
IF(MI.GT.MJ) GO TO 260 

TZ=MIl*BJ(MMP)*XDSM 

GO TO 270 

260 TZ=MJl*BJ(MXP)*XDSM 

270 UIJ=MJl*(Tl-TZ) 

RETURN 

END 

A-29 
SUBROUT I NE STE I (BETA, SS, TT 1 

c....................................................................... 

c: Evaluates equations 12 and 13 for i-th aperture. 

c:.....................................................................: 

IMPLICIT REAL*8(A-H,O-Z) 

DIMENSION BJ(21) 

COMMON /FACTS/ MI, MJ, MP, MM, MPP, MMP, MI 1, MJ1, MXP, MPI, MPJ, MPPI, MPPJ, 

A 

B XESI,XESJ,XMSI,XMSJ,SFI,SFJ,TFI,TFJ,TFEI,TFEJ, 

C 

EM1, EM2, AIKO, AJKO,RKO, XEI, XEJ, XMI, XMJ, 

XDS, XDC, XDSP , XDCP , XDSM, XDCM 

ARG=AIKO*BETA 

CALL BESJS(ARG, BJ, MPPI 1 

BJA=BJ(MPI)/ARG 

SS-SFI'BJA 

DNzXESI-ARG*ARG 

IF(DABS(DN).GT. 1.OD-05) GO TO 10 

CALL BESJS(XE1, BJ, MPI 

TT=TFEI*BJ(MPI) 

RETURN 

RETURN 

END 

10 TT=TFI*(MI*BJA-BJ(MPPI))/DN 

SUBROUT I NE STM I (BETA, SS, TT 1 

c...................................................................... . 

c: Evaluates equations 14 and 15 for i-th aperture. 

c:.....................................................................: 

IMPLICIT REAL*8(A-H,O-Z) 

DIMENSION BJ(21) 

COMMON /FACTS/ MI, MJ, MP, MM, MPP, MMP, MI1, MJ1, MXP, MPI, MPJ, MPPI, MPPJ, 

A 

B XESI,XESJ,XMSI,XMSJ,SFI,SFJ,TFI,TFJ,TFEI,TFEJ, 

C 

EM1, EM2, AIKO, AJKO, RKO, XEI, XEJ, XMI, XMJ, 

XDS, XDC , XDSP , XDCP , XDSM, XDCM 

TT=O. OD0 

ARG=AIKO*BETA 

DN=XMSI-ARG*ARG 

IF(DABS(DN).GT. 1.OD-05) GO TO 10 

CALL BESJS(XM1, BJ, MPPI 

SS=AIKO*0.5DO*BJ(MPPI) 

RETURN 

CALL BESJS(ARG, BJ, MPI 1 

SS=AIKO*ARG*BJ(MPI)/DN 

RETURN 

END 

10 

A-30 
. SUBROUT I NE STE J ( BETA, SS , TT 1 ........................................................................ 

c: Evaluates equations 12 and 13 for j-th aperture. ........................................................................ 

IMPLICIT REAL*8(A-H,O-Z) 

DIMENSION BJ(21) 

COMMON /FACTS/ MI,MJ,MP,MM,MPP,MMP,MIl,MJl,MXP,MPI,MPJ,MPPI,MPPJ, 

A 

B 

C 

EM1, EM2, AIKO, AJKO,RKO, XEI,XEJ, XMI,XMJ, 

XESI , XESJ, XMSI , XMSJ, SFI , SFJ , TFI , TFJ , TFEI , TFEJ, 

XDS, XDC, XDSP, XDCP, XDSM, XDCM 

ARG=AJKO*BETA 

CALL BESJS(ARG,BJ,MPPJ) 

BJA=BJ(MPJ)/ARG 

SS=SFJ*BJA 

DN=XESJ-ARG*ARG 

IF(DABS(DN). GT. 1.OD-05) GO TO 10 

CALL BESJS(XEJ,BJ,MPJ) 

TT=TFEJ*BJ(MPJ) 

RETURN 

RETURN 

END 

10 TT=TFJ*(MJ*BJA-BJ(MPPJ))/DN 

SUBROUTINE STMJ(BETA,SS,TT) 

c....................................................................... 

c: Evaluates equations 14 and 15 for j-th aperture. 

c: ..................................................................... : 

IMPLICIT REAL*8(A-H,O-Z) 

DIMENSION BJ(21) 

COMMON /FACTS/ MI,MJ,MP,MM,MPP,MMP,MIl,MJl,MXP,MPI,MPJ,MPPI,MPPJ, 

A 

B 

C 

EM1, EM2, AIKO, AJKO,RKO, XEI, XEJ, XMI, XMJ, 

XESI, XESJ, XMSI, XMSJ, SFI, SFJ, TFI, TFJ, TFEI, TFEJ, 

XDS, XDC, XDSP, XDCP, XDSM, XDCM 

TT=O. OD0 

ARG=AJKO*BETA 

DNzXMSJ-ARG*ARG 

IF(DABS(DN). GT. 1.OD-05) GO TO 10 

CALL BESJS(XMJ,BJ,MPPJ) 

SS=AJKO*0.5DO*BJ(MPPJ) 

RETURN 

10 CALL BESJS(ARG, BJ, MPJ) 

SS=AJKO*ARG*BJ(MPJ)/DN 

RETURN 

END 

A-31 
SUBROUTINE BESJS(XX, BJ, MP) ........................................................................ 

c: Evaluates Bessel function of first kind or orders 0 through MP-1 : 

c: for argument XX, and stores values in BJ(1) through BJ(MP). 

c: BJ must be dimensioned at least (MP+l) in calling program. 

........................................................................ 

100 

110 

120 

130 

140 

150 

160 

170 

180 

190 

200 

210 

220 

IMPLICIT REAL*8 (A-H,O-Z) 

DIMENSION BJ(l),B(130) 

DATA PI,PI02 /3.14159265358979D0,1.5707963267949ODO/ 

DATA BIG, SMALL, PI04 /l. OD+38,l.OD-38,0.78539816339745DO/ 

DATA TWOTRD /0.666666666666667DO/ 

X=DABS(XX) 

BJ(l)=l. OD0 

MT=MP- 1 

N=MT+ 1 

IF(X.LT.80. )GO TO 110 

DO 100 I=l,N 

BJ(I)~DSQRT(2.ODO/(PI*X~~*DCOS(X-PIO4-PIO2*~1-1~~ 

GO TO 310 

CONT I NUE 

DO 120 I-2,N 

BJ(I)=O.ODO 

IF (X. LT. 1. D-06 )RETURN 

IF(X-15. )130,130,140 

NTEST=ZO.ODO+ 10. ODO*X-X**TWOTRD 

GO TO 150 

NTEST=90.ODO+O. 5DO*X 

IF(MT-NTEST)l70,160,160 

NzNTEST- 1 

GO TO 180 

N=MT 

BPREV=O. 0 

Nl=N+l 

F=2.ODO/X 

D=l . OD-06 

IF(X-5. ) 190,200,200 

MA=X+6.ODO 

GO TO 210 

MA=1.4DO*X+60.ODO/X 

I F I XX=X 

MB=N+IFIXX/4+2 

MZERO=MAXO (MA, MB 1 

MMAX=NTEST 

DO 280 M=MZERO, MMAX, 3 

FM 1 =SMALL 

FM=O. OD0 

ALPHA=O. OD0 

IF(M-(M/2)*2)230,220,230 

JT=- 1 

GO TO 240 

A-32 
230 

240 

250 

260 

270 

280 

290 

300 

310 

320 

330 

JT= 1 

M2=M-2 

DO 260 K=l,M2 

MK=M-K 

B(MK)=F*MK*FMl-FM 

IF(B(MK)-BIG)Z50,310,310 

FM=FMl 

FMl=B(MK) 

JT=-JT 

S=l+JT 

ALPHA=ALPHA+B(MK)*S 

B(l)=F*FMl-FM 

ALPHA=ALPHA+B(l) 

BTEST=B(Nl)/ALPHA 

IF(DABS(BTEST-BPREV)-DABS(D*BTEST))290,290,270 

BPREV=BTEST 

CONT I NUE 

DO 300 I=l,Nl 

BJ(I)=B(I)/ALPHA 

IF(XX. LT. 0. ODOIGO TO 320 

RETURN 

N=MT+ 1 

DO 330 I=l,N 

NN=I - 1 

BJ(I)=BJ(I)*(-l.O)**NN 

RETURN 

END 

A-33 
Report Documentation Page 

1. Report No. 2. Government Accession No. 1 3. Recipient's Catalog No. 

NASA TM- 101614 

4. Title and Subtitle 5. Report Date 

CWG - A Fortran Program for Mutual Coupling in a 

P1 anar Array of Ci rcul ar Wavegui de-fed Apertures June 1989 

6. Performing Organization Code 

16. Abstract 

7. Author(s) 

This paper is a documentation of a Fortran program which calculates the mutual 

coupling between circular apertures in a conductive plane. 

quite general in that the apertures do not have to be the same sizes, nor do they 

have to be polarized in the same direction. In addition, several waveguide 

modes (TE and/or TM) may be specified in the apertures and the mutual coupling 

between all combinations of apertures and modes will be calculated. The program 

also allows multiple layers of homogeneous dielectrics to be placed over the 

aperture array. Outside the layered region, one can specify either a homogeneous 

half-space, or a perfect reflecting surface. 

The program is 

8. Performing Organization Report No. 

TWork Unit No. 

9. Performing Organization Name and Address 

NASA Langley Research Center 

Hampton, VA 23665-5225 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 

Washington, D.C. 20546-0001 

506-44-21-03 

11. Contract or Grant No. 

13. Type of Report and Period Covered 

Technical Memorandum 

14. Sponsoring Agency Code 

17. Key Words (Suggested by Author(s)) 18. Distribution Statement 

antenna 

apertures 

arrays 

coup1 i ng 

Unclassified - Unlimited 

Subject Category 32 

19. Security Classif. (of this report1 20. Security Classif. (of this pagel 21. No. of pages 

unclassified uncl assi f i ed 45 

22. Price 

A03