c this is a standard Fortran 77 program for calculating two apertures 
c on an infinite extention conducting plane backed waveguides and network. 
c In the program, the lines which begin with 'c*' are for user to get some 
c interested parameters. Only take out 'c*' of the corresponding lines and 
c compile the program. Run that again, then you can get the interesting 
c parameters from the corresponding output files.   
	IMPLICIT NONE
	COMPLEX*8 YMNA
	COMPLEX*8 YMNB(100,100),YMN(100,100),AIM(100),AI(50),BI(50),
     1  VM(100),CMAT11(50),CMAT12(50),CMAT21(50),CMAT22(50),ABMOD(50)
     2  ,BMOD(50),YL(2,2),YL11(50),YL12(50),YL21(50),YL22(50),CMAT(2,2)
     3  ,YMNAA(100,100)
        COMPLEX*8 HMAG(5000),HSC(5000),HTAT(5000),HSCT(5000)
	COMPLEX*8 AG(5000),HSCN(5000)
	COMPLEX*8 HSCRR,HFIELD,SMAG,ZJ,FMAG,FARH
	REAL*4  AIR(100),AII(100),AIA(100),AVR(100),AVI(100)
        REAL*4 HMGR(5000),HMGI(5000),HMGA(5000),HTTR(5000),HTTI(5000)
        REAL*4 HTTA(5000),HSCR(5000),HSCI(5000),HSCA(5000),ANG(5000)
	REAL*4 AVA(100),AGA(5000),HSNR(5000)  
	REAL TIME(2)  
	REAL DTIME,TOTAL,STIME  
	REAL*4 UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,WL,R0,ANG0,FAI,AA1
     1  ,AA2,BB1,BB2,X0,Y0
	INTEGER I,J,M,N,KI,KJ,IJ,KJJ,KII,II,IA,IK,JK
	COMMON/T1/ I,J,M,N
        COMMON/T2/ UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,WL,FAI
C OPEN BAPIN.DAT ON 10 FOR IMPUT.
C OPEN BAPOUT.DAT ON 11 FOR OUTPUT  
	STIME=DTIME(TIME)  
	OPEN (30,FILE='emin1.dat',STATUS='OLD')
	OPEN (12,FILE='emin2.dat',STATUS='OLD')
	OPEN (14,FILE='emin3.dat',STATUS='OLD')
	OPEN (16,FILE='emin4.dat',STATUS='OLD')
c* 	OPEN (9,FILE='emout0.dat',STATUS='UNKNOWN')
	OPEN (11,FILE='emsout1.dat',STATUS='UNKNOWN')
c*	OPEN (13,FILE='emout2.dat',STATUS='UNKNOWN')
c*	OPEN (15,FILE='emout3.dat',STATUS='UNKNOWN')
c*	OPEN (17,FILE='emout4.dat',STATUS='UNKNOWN')
c*	OPEN (19,FILE='emout5.dat',STATUS='UNKNOWN')
c*	OPEN (21,FILE='emout6.dat',STATUS='UNKNOWN')
c*	OPEN (55,FILE='emout7.dat',STATUS='UNKNOWN')
c*	OPEN (39,FILE='emout51.dat',STATUS='UNKNOWN')
c*	OPEN (31,FILE='emout61.dat',STATUS='UNKNOWN')
c*	OPEN (65,FILE='emout71.dat',STATUS='UNKNOWN')
	OPEN (49,FILE='emsout2.dat',STATUS='UNKNOWN')
c*	OPEN (41,FILE='emout62.dat',STATUS='UNKNOWN')
c*	OPEN (75,FILE='emout72.dat',STATUS='UNKNOWN')
c*	OPEN (92,FILE='emout8.dat',STATUS='UNKNOWN')
c*	OPEN (96,FILE='emout9.dat',STATUS='UNKNOWN')
 1	FORMAT(4I3)
 2	FORMAT(9(E12.5))
 3	FORMAT(2(E12.5,5X,E12.5))
 4	FORMAT(2X,'KI=',I2,5X,'KJ=',I2,5X,'Ymnb=',E12.5,2X,'+j',2X,E12.5)
 5 	FORMAT(2X,'KI=',I2,2X,'kj=',I2,5X,'Ymn=',E17.10,2X,'+j',2X,E17.10)
 6   	FORMAT(2X,'M=',I2,5X,'Vm=',E12.5,2X,'+j',2X,E12.5)
 7	FORMAT(4(E12.5,E12.5))
 8	FORMAT(2X,'HMAG=',E12.5,2X,'+J',2X,E12.5)
 9	FORMAT(2X,'HSC=',E12.5,2X,'+J',2X,E12.5)
 10	FORMAT(2X,'HSCAT=',E12.5,2X,'+J',2X,E12.5)
 11	FORMAT(2X,'HTAT=',E12.5,2X,'+J',2X,E12.5)
C READ IN ORIGINAL DATA.
	READ(30,*) I,J,M,N
	READ(12,*) UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2
	READ(14,*) WL,R0,ANG0,FAI
	PRINT*,'I=',I,'J=',J,'M=',M,'N=',N
	PRINT*,'UM=',UM1,UM2,'EM=',EM1,EM2,'Habd=',HA,HB,HD,'L=',AL1,AL2
	PRINT*,'WL=',WL,'R0=',R0,'ANG0=',ANG0,'FAI=',FAI
	WRITE(11,12) I,J,M,N
	WRITE(11,13) UM1,UM2,EM1,EM2
	WRITE(11,14) HA,HB,HD
	WRITE(11,15) AL1,AL2,WL
	WRITE(11,16) R0,ANG0,FAI
 12	FORMAT(10X,'I=',I3,5X,'J=',I3,5X,'M=',I3,5X,'N=',I3)
 13     FORMAT(2X,'UM1=',E12.5,2X,'UM2=',E12.5,2X,'EM1=',E12.5
     1,2X,'EM2=',E12.5)
 14	FORMAT(2X,'A=',E12.5,5X,'B=',E12.5,5X,'D=',E12.5)
 15	FORMAT(2X,'L1=',E12.5,5X,'L2=',E12.5,5X,'WL=',E12.5)
 16     FORMAT(2X,'R0=',E12.5,5X,'ANG0=',E12.5,5X,'FAI=',E12.5)
c calculate eguivalent microwave network parameters
	DO 600 KI=1,I+J
	READ(16,*) YL11(KI),YL12(KI),YL21(KI),YL22(KI)
	YL(1,1)=YL11(KI)
	YL(1,2)=YL12(KI)
	YL(2,1)=YL21(KI)
	YL(2,2)=YL22(KI)
	CALL CNETWORK(YL,CMAT,KI)
	CMAT11(KI)=CMAT(1,1)
	CMAT12(KI)=CMAT(1,2)
	CMAT21(KI)=CMAT(2,1)
	CMAT22(KI)=CMAT(2,2)
 600	CONTINUE
	DO 606 KI=1,I+J
	DO 650 KJ=1,I+J
	IF((KI.EQ.1).AND.(KJ.EQ.1)) YMNB(1,1)=CMAT11(1)
	IF((KI.EQ.1).AND.(KJ.EQ.I+1)) YMNB(1,I+1)=CMAT12(1)
	IF((KI.EQ.I+1).AND.(KJ.EQ.1)) YMNB(I+1,1)=CMAT21(1)
	IF((KI.EQ.I+1).AND.(KJ.EQ.I+1)) YMNB(I+1,I+1)=CMAT22(1)
	IF((KI.EQ.KJ).AND.(KI.LT.I+1)) YMNB(KI,KJ)=CMAT11(KI)
	IF((KI.EQ.KJ).AND.(KI.GT.I+1)) YMNB(KI,KJ)=CMAT22(KI)
	WRITE(11,4)KI,KJ,YMNB(KI,KJ)
 650	CONTINUE
 606	CONTINUE
	DO 101 KI=1,I
	DO 102 KJ=1,I
	CALL APULSE(KI,KJ,M,M,HA,HA,YMNA,HD,HD)
	YMNAA(KI,KJ)=YMNA
	YMN(KI,KJ)=YMNAA(KI,KJ)+YMNB(KI,KJ)
	WRITE(11,5) KI,KJ,YMN(KI,KJ)
 102	CONTINUE
 101	CONTINUE
	IJ=I+J
	DO 201 KI=1,I
	DO 202 KJ=I+1,IJ
	KJJ=KJ-I
	CALL APULSE(KI,KJJ,M,N,HA,HB,YMNA,HD,-HD)
	YMNAA(KI,KJ)=YMNA
	YMN(KI,KJ)=YMNAA(KI,KJ)+YMNB(KI,KJ)
	WRITE(11,5) KI,KJ,YMN(KI,KJ)
 202	CONTINUE
 201	CONTINUE
	DO 301 KI=I+1,IJ
	KII=KI-I
	DO 302 KJ=1,I
	CALL APULSE(KII,KJ,N,M,HB,HA,YMNA,-HD,HD)
	YMNAA(KI,KJ)=YMNA
	YMN(KI,KJ)=YMNAA(KI,KJ)+YMNB(KI,KJ)
	WRITE(11,5) KI,KJ,YMN(KI,KJ)
 302	CONTINUE
 301	CONTINUE
	DO 401 KI=I+1,IJ
	KII=KI-I
	DO 402 KJ=I+1,IJ
	KJJ=KJ-I
	CALL APULSE(KII,KJJ,N,N,HB,HB,YMNA,-HD,-HD)
	YMNAA(KI,KJ)=YMNA
	YMN(KI,KJ)=YMNAA(KI,KJ)+YMNB(KI,KJ)
	WRITE(11,5) KI,KJ,YMN(KI,KJ)
 402	CONTINUE
 401	CONTINUE
	CALL CURRENT(AI,BI)
	DO 501 KI=1,I
	AIM(KI)=AI(KI)
	WRITE(11,21) KI,AIM(KI)
 21	FORMAT(2X,'I=',I3,5X,'I(i)=',2X,E12.5,2X,'+I'2X,E12.5)
 501	CONTINUE
	DO 502 KJ=I+1,IJ
	KJJ=KJ-I
	AIM(KJ)=BI(KJJ)
	WRITE(11,22) KJ,AIM(KJ)
 22	FORMAT(2X,'I=',I3,5X,'I(i)=',2X,E12.5,2X,'+I'2X,E12.5)
 502	CONTINUE
	CALL GAUSS(IJ,YMN,AIM,VM)
	DO 171 II=1,IJ
	WRITE(11,6) II,VM(II)
 171	CONTINUE
	DO 123 II=1,IJ
	AIR(II)=REAL(AIM(II))
	AII(II)=AIMAG(AIM(II))
	AIA(II)=ABS(AIM(II))
	AVR(II)=REAL(VM(II))
	AVI(II)=AIMAG(VM(II))
	AVA(II)=ABS(VM(II))
	WRITE(13,*) II,AIR(II),AII(II),AIA(II)
	WRITE(15,*) II,AVR(II),AVI(II),AVA(II)
 123	CONTINUE
	ZJ=(0.,1.)
	AA1=HD
	BB1=HD+HA
	AA2=-HD-HB
	BB2=-HD
	ABMOD(1)=VM(1)/SQRT(HA)
	BMOD(1)=VM(I+1)/SQRT(HB)
	DO 306 KI=2,I
	ABMOD(KI)=VM(KI)*SQRT(2./HA)
 306	CONTINUE
	DO 405 KJ=2,J
	BMOD(KJ)=VM(I+KJ)*SQRT(2./HB)
 405	CONTINUE
	HTAT(1)=(0.,0.)
	HSC(1)=(0.,0.)
	HSCN(1)=(0.,0.)
	DO 777 IA=1,901
	ANG(IA)=(IA-1)*3.49065852E-3+ANG0
	X0=R0*COS(ANG(IA))
	Y0=R0*SIN(ANG(IA))
	HFIELD=(0.0,0.0)
	FARH=(0.,0.)
	DO 666 IK=1,I
	CALL HF(IK,WL,X0,Y0,AA1,BB1,HA,HD,SMAG)
	HFIELD=HFIELD+ABMOD(IK)*SMAG
	CALL AFHF(IK,WL,R0,ANG(IA),AA1,BB1,HA,HD,FMAG)
	FARH=FARH+ABMOD(IK)*FMAG
 666	CONTINUE
	DO 888 JK=1,J
	CALL HF(JK,WL,X0,Y0,AA2,BB2,HB,-HD-HB,SMAG)
	HFIELD=HFIELD+BMOD(JK)*SMAG
	CALL AFHF(JK,WL,R0,ANG(IA),AA2,BB2,HB,-HB-HD,FMAG)
	FARH=FARH+BMOD(JK)*FMAG
 888	CONTINUE
	HMAG(IA)=2.*3.141592654*3.E8/WL*8.841941283E-12*HFIELD
	AG(IA)=2.*3.141592654*3.E8/WL*8.841941283E-12*FARH
	HSCRR=-ZJ*2.*3.141592654/WL*(X0*COS(FAI)+Y0*SIN(FAI))
	HSC(IA)=EXP(HSCRR)
	HSCN(IA)=EXP(-1.*HSCRR)
	HSNR(IA)=ABS(HSCN(IA))
	HTAT(IA)=HMAG(IA)+HSC(IA)+HSCN(IA)
	HSCT(IA)=HSCN(IA)+HMAG(IA)
	HMGR(IA)=REAL(HMAG(IA))
	HMGI(IA)=AIMAG(HMAG(IA))
	HMGA(IA)=ABS(HMAG(IA))
	HTTR(IA)=REAL(HTAT(IA))
	HTTI(IA)=AIMAG(HTAT(IA))
	HTTA(IA)=ABS(HTAT(IA))
	HSCR(IA)=REAL(HSCT(IA))
	HSCI(IA)=AIMAG(HSCT(IA))
	HSCA(IA)=ABS(HSCT(IA))
	AGA(IA)=ABS(AG(IA))
c*	WRITE(9,*) HMAG(IA),HSC(IA)
c*	WRITE(17,*) HSCT(IA),HTAT(IA)
c*	WRITE(19,*) HMGR(IA)
c*	WRITE(39,*) HMGI(IA)
	WRITE(49,*) HMGA(IA)
c*	WRITE(21,*) HTTR(IA)
c*	WRITE(31,*) HTTI(IA)
c*	WRITE(41,*) HTTA(IA)
c*	WRITE(55,*) HSCR(IA)
c*	WRITE(65,*) HSCI(IA)
c*	WRITE(75,*) HSCA(IA)
c*	WRITE(92,*) HSNR(IA)
c*	WRITE(96,*) AGA(IA)
	HSC(IA+1)=(0.,0.)
	HSCN(IA+1)=(0.,0.)
	HTAT(IA+1)=(0.,0.)
 777	CONTINUE  
	TOTAL=DTIME(TIME)-STIME  
	PRINT*,'Time=',TOTAL   
	STOP
778	END
C
C
C     SUBPROGRAM
C     TO SOLVE A LINEAR COMPLEX SYSTEM
C     N:    THE DIMENSION OF THE SYSTEM
C     UL:   THE INPUT COMPLEX MATRIX
C     B:    THE INPUT COMPLEX VECTOR
C     X:    THE OUTPUT COMPLEX VECTO
C            UL*X = B 
C     ***********************************************************
C     TO DO LU-DECOMPOSITION 
C                                                                               
      SUBROUTINE DECOMP(N,IPS,UL)
        IMPLICIT NONE
        INTEGER NC
	PARAMETER (NC=100)
	COMPLEX*8 UL(NC,NC),PIVOT,EM
      REAL*4 SCL(NC)
	INTEGER IPS(NC),N,I,J,K,NM1,IP,IPV,IPP,KP1
	REAL*4 RN,ULM,BIG,SIZE
	DO 5 I=1,N
      IPS(I)=I
      RN=0.
      DO 2 J=1,N
      ULM=ABS(REAL(UL(I,J)))+ABS(AIMAG(UL(I,J)))
      IF (RN-ULM) 1,2,2
1     RN=ULM
2     CONTINUE
      SCL(I)=1./RN
5     CONTINUE
      NM1=N-1
      DO 17 K=1,NM1
      BIG=0.
      DO 11 I=K,N
      IP=IPS(I)
      SIZE=(ABS(REAL(UL(IP,K)))+ABS(AIMAG(UL(IP,K))))*SCL(IP)
      IF(SIZE-BIG) 11,11,10
10    BIG=SIZE
      IPV=I
11    CONTINUE
      IF(IPV-K) 14,15,14
14    J=IPS(K)
      IPS(K)=IPS(IPV)
      IPS(IPV)=J
15    IPP=IPS(K)
      PIVOT=UL(IPP,K)
      KP1=K+1
      DO 16 I=KP1,N
      IP=IPS(I)
      EM=-UL(IP,K)/PIVOT
18    UL(IP,K)=-EM
      DO 16 J=KP1,N
      UL(IP,J)=UL(IP,J)+EM*UL(IPP,J)
16    CONTINUE
17    CONTINUE
      RETURN
      END       
C     **************************************************************** 
C     TO SOLVE    UL*X = B
C
C
      SUBROUTINE GAUSS(N,UL,B,X)
        IMPLICIT NONE
        INTEGER NC
	PARAMETER (NC=100)
      COMPLEX*8 UL(NC,NC),B(NC),X(NC),SUM
      INTEGER IPS(NC)
	INTEGER N,IP,I,J,IM1,IBAOK,IP1,NP1
      CALL DECOMP(N,IPS,UL)
      NP1=N+1
      IP=IPS(1)
      X(1)=B(IP)
      DO 2 I=2,N
      IP=IPS(I)
      IM1=I-1
      SUM=(0.0,0.0)
      DO 1 J=1,IM1
1     SUM=SUM+UL(IP,J)*X(J)
2     X(I)=B(IP)-SUM
      IP=IPS(N)
      X(N)=X(N)/UL(IP,N)
      DO 4 IBAOK=2,N
      I=NP1-IBAOK
      IP=IPS(I)
      IP1=I+1
      SUM=(0.0,0.0)
      DO 3 J=IP1,N
3     SUM=SUM+UL(IP,J)*X(J)
4     X(I)=(X(I)-SUM)/UL(IP,I)
      RETURN
      END
C 
C 
C                                                                               
C 
	SUBROUTINE CURRENT(AI,BI)
C this is a subroutine for calculating the excitation vector
	IMPLICIT NONE 
	COMPLEX*8 AI(50),BI(50)
	COMPLEX*8 BX,BID,ZJ,AX,AAI,BBI,AID
	INTEGER I,J,M,N,KI,KJ,KN
	REAL*4 UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,WL,FAI,DELT1,DELT2
     1,HK,AKM,AKN,TES,EFC1,EFC2,RXM,CXKI,AAX,RXN,CXKJ,AIM,BIM,BBX
	COMMON/T1/ I,J,M,N
        COMMON/T2/ UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,WL,FAI
	ZJ=(0.,1.)
	DELT1=HA/M
	DELT2=HB/N
	HK=6.283185307/WL
	AKM=HK*DELT1
	AKN=HK*DELT2
	TES=ABS(COS(FAI))
	IF(TES.LT.1.E-10) GOTO 987
	EFC1=2.*SIN(AKM*COS(FAI))/(HK*COS(FAI))
	EFC2=2.*SIN(AKN*COS(FAI))/(HK*COS(FAI))
	GOTO 988
 987	EFC1=2.*DELT1
	EFC2=2.*DELT2
 988	DO 1000 KI=1,I
	AAI=(0.0,0.0)
	DO 100 KJ=1,M
	RXM=(KJ+0.5)*DELT1
	CXKI=KI*3.141592654*RXM/HA
	IF(KI.EQ.1) AIM=1./SQRT(HA)
	IF(KI.GT.1) AIM=SQRT(2./HA)*COS(CXKI)
	AAX=HK*(RXM+HD)*COS(FAI)
	AX=-ZJ*AAX
	IF(AAX.LE.-7.8E1) GOTO 10
	AID=AIM*EXP(AX)*EFC1
	GOTO 15
 10	AID=(0.0,0.0)
 15	AAI=AAI+AID
 100	CONTINUE
	AI(KI)=AAI
 1000	CONTINUE
	DO 2000 KJ=1,J
	BBI=(0.,0.)
	DO 200 KN=1,N
	RXN=-(KN+0.5)*DELT2
	CXKJ=KJ*3.141592654*RXN/HB
	IF(KJ.EQ.1) BIM=1./SQRT(HB)
	IF(KJ.GT.1) BIM=SQRT(2./HB)*COS(CXKJ)
	BBX=HK*(RXN-HD)*COS(FAI)
	BX=ZJ*BBX
	IF(BBX.LE.-7.8E1) GOTO 20
	BID=BIM*EXP(BX)*EFC2
	GOTO 25
 20	BID=(0.,0.)
 25	BBI=BBI+BID
 200	CONTINUE
	BI(KJ)=BBI
 2000	CONTINUE
	RETURN 
	END
C
C
	SUBROUTINE APULSE(KI,KJ,KN,KM,HW1,HW2,YMNA,START1,START2)
C this is a subroutine for calculating the 'A' region equivalence admittances 
C by using the pulse basis mode function expansion.  
	IMPLICIT NONE
	COMPLEX*8 YMND,YMNA,SOUT,SOUT1
	REAL*4 HW1,HW2,XM,XN,GMN,UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,
     1WL,FAI,UM0,EM0,A,B,C,AD,YMNDI,SGN,START1,START2
	REAL*4 YMNI(550),YMNJ(550)
	INTEGER I,J,M,N,MN,MM,KI,KJ,KN,KM
	COMMON/T1/ I,J,M,N
        COMMON/T2/ UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,WL,FAI
	EXTERNAL SGN
C this is for calculating mode basis coefficients
	UM0=1.256637061E-6
	EM0=8.841941283E-12
	DO 3300 MN=1,KN	
	XN=(MN+0.5)*HW1/KN*SGN(START1)
	IF(KI.EQ.1) YMNI(MN)=1./SQRT(HW1)
	IF(KI.GT.1) YMNI(MN)=SQRT(2./HW1)*COS((KI-1)*3.141592654*XN/HW1)
 3300	CONTINUE
	DO 3000 MM=1,KM
	XM=(MM+0.5)*HW2/KM*SGN(START2)
	IF(KJ.EQ.1) YMNJ(MM)=1./SQRT(HW2)
	IF(KJ.GT.1) YMNJ(MM)=SQRT(2./HW2)*COS((KJ-1)*3.141592654*XM/HW2)
 3000	CONTINUE
C this is for calculating admittance of pulse basis
	YMNA=(0.0,0.0)
	DO 2200 MN=1,KN
	DO 2020 MM=1,KM
	IF((START1.EQ.START2).AND.(MM.EQ.MN)) GOTO 1100
	A=START1+(MN-1)*HW1/KN*SGN(START1)
	B=START1+MN*HW1/KN*SGN(START1)
	C=START2+(MM-1)*HW2/KM*SGN(START2)
	AD=START2+MM*HW2/KM*SGN(START2)
	CALL CGAUS(A,B,C,AD,WL,SOUT)
	SOUT=SOUT*SGN(START1)*SGN(START2)
	GOTO 1001
 1100	GMN=SQRT(UM1*EM1)*2.*3.141592654/WL*HW1/KN*0.1638049574
	YMNDI=0.6366197725D0*ALOG(GMN)
	SOUT1=CMPLX(1.0,YMNDI)
	SOUT=SOUT1*HW1*HW2/KM/KN
 1001	YMND=YMNI(MN)*YMNJ(MM)*SOUT/WL*1.66663E-2
	YMNA=YMNA+YMND
 2020	CONTINUE
 2200	CONTINUE
	RETURN
	END
C
C
      SUBROUTINE CGAUS(A,B,C,AD,WL,SOUT)
C      IMPLICIT NONE      
       COMPLEX*8 H20,SS,ZZ,SOUT
	REAL*4 A,B,C,AD,XM,XR,ZM,ZR,DELTAX,DELTAZ
        REAL*4 X(5),W(5)
	REAL*4 A1,A2,A3,A4,WL
	EXTERNAL H20 
      DATA X/0.1488743389,0.4333953941,
     1          0.6794095682,0.860633666,0.9739065285/
      DATA W/0.2955242247,0.2692667193,
     2          0.2190863625,0.1494513491,0.666713443E-1/
      XM=0.5*(B+A)
      XR=0.5*(B-A)
      ZM=0.5*(AD+C)
      ZR=0.5*(AD-C)
      SS=CMPLX(0.0,0.0)
      DO 11 J=1,5
       DELTAZ=ZR*X(J)
       ZZ=CMPLX(0.0,0.0)
      DO 12 I=1,5
        DELTAX=XR*X(I)
	A1=6.283185307/WL*ABS(XM-ZM+DELTAX+DELTAZ)
	A2=6.283185307/WL*ABS(XM-ZM+DELTAX-DELTAZ)
	A3=6.283185307/WL*ABS(XM-ZM-DELTAX+DELTAZ)
	A4=6.283185307/WL*ABS(XM-ZM-DELTAX-DELTAZ)
	ZZ=ZZ+W(I)*(H20(A1)+H20(A2)+H20(A3)+H20(A4))
12      CONTINUE
	SS=SS+W(J)*ZZ*XR
11    CONTINUE
      SOUT=ZR*SS
      RETURN
      END
C
C
        COMPLEX*8 FUNCTION H20(A)
	IMPLICIT NONE      
	REAL*4 BJ,BY,F0,T0,A
	IF(A.GT.3.0) GO TO 5
      BJ=(A/3.0)**2
      BJ=1.0+BJ*(-2.2499997+BJ*(1.2656208+BJ*
     1(-0.3163866+BJ*(4.44479E-2+BJ*(-3.9444E-3+BJ*
     22.1E-4)))))
      BY=(A/3.0)**2
      BY=2.0/3.14159265*ALOG(A/2.0)*BJ+0.36746691
     1+BY*(0.60559366+BY*(-0.74350384+BY*(0.25300117+BY*(
     2-4.261214E-2+BY*(4.27916E-3-BY*2.4846E-4)))))
      GO TO 10
5     BJ=3.0/A
      F0=0.79788456+BJ*(-7.7E-7+BJ*(-5.52740E-3+BJ*(-9.512E-5
     1+BJ*(1.37237E-3+BJ*(-7.2805E-4+BJ*1.4476E-4)))))
      T0=A-0.78539816+BJ*(-4.166397E-2+BJ*(-3.954E-5+BJ*(2.62573E-3
     1+BJ*(-5.4125E-4+BJ*(-2.9333E-4+BJ*1.3558E-4)))))
      BY=SQRT(A)
      BJ=F0*COS(T0)/BY
      BY=F0*SIN(T0)/BY
10    CONTINUE
      H20=CMPLX(BJ,-BY)
      RETURN
      END
C
C
	REAL*4 FUNCTION SGN(X)
	IMPLICIT NONE 
	REAL*4 X
	IF(X.EQ.0.) SGN=0.
	IF(X.LT.0.) SGN=-1.
	IF(X.GT.0.) SGN=1.
	RETURN
	END
C
C
C
	SUBROUTINE CNETWORK(YL,CMAT,NK)
	IMPLICIT NONE
	COMPLEX*8 YL(2,2),CMAT(2,2)
	COMPLEX*8 SINH,COSH,YL11,YL12,YL21,YL22,Z,GB1R,GB1,GB2R,GB2
     1,GA1,GAMA1,GA2,GAMA2,BAMT1,BAMT2,BAMT3,BAMT4,DELT1,DELT2,
     2BMT1,BMT2,BMT3,BMT4,DMAT11,DMAT12,DMAT13,DMAT14,DMAT1,DMBT11,
     1DMBT12,DMBT13,DMBT14,DMBT1,DMAT21,DMAT22,DMAT23,DMAT24,DMAT2,
     1DMBT21,DMBT22,DMBT23,DMBT24,DMBT2,ZJ
	REAL*4 UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,WL,FAI,HK1,HK2,ETA11
     1,ETA22,CB
	INTEGER I,J,M,N,NN,NK
	COMMON/T1/ I,J,M,N
        COMMON/T2/ UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,WL,FAI
	SINH(Z)=0.5*(EXP(Z)-EXP(-Z))
	COSH(Z)=0.5*(EXP(Z)+EXP(-Z))
	NN=NK-1
	YL11=YL(1,1)
	YL12=YL(1,2)
	YL21=YL(2,1)
	YL22=YL(2,2)
	HK1=6.283185307/WL*SQRT(EM1*UM1)
	HK2=6.283185307/WL*SQRT(EM2*UM2)
	ETA11=SQRT(UM1/EM1)*377
	ETA22=SQRT(UM2/EM2)*377
	IF(NN.EQ.0) GOTO 1600
	IF(NN.EQ.I) GOTO 1700
	ZJ=(0.,1.)
	IF(NN.LT.I) CB=NN*WL/2.
	IF(NN.GT.I) CB=(NN-I)*WL/2.
	GB1R=CMPLX(1.-((HA/CB)**2)*(UM1*EM1),0.0)
	GB1=1./SQRT(GB1R)
	GB2R=CMPLX(1.-((HB/CB)**2)*(UM2*EM2),0.0)
	GB2=1./SQRT(GB2R)
	CMAT(1,1)=1./ETA11*GB1*ZJ
	CMAT(2,2)=1./ETA22*GB2*ZJ
	CMAT(1,2)=(0.,0.)
	CMAT(2,1)=(0.,0.)
	GOTO 1700
 1600	GA1=CMPLX(HK1*HK1-NN*NN*9.869604401/(HA*HA),0.0)
	GAMA1=SQRT(GA1)
	GA2=CMPLX(HK2*HK2-NN*NN*9.869604401/(HB*HB),0.0)
	GAMA2=SQRT(GA2)
	BAMT1=COSH(GAMA1*AL1)+ETA11*YL11*SINH(GAMA1*AL1)
	BAMT2=YL12*ETA22*SINH(GAMA2*AL2)
	BAMT3=ETA11*YL21*SINH(GAMA1*AL1)
	BAMT4=COSH(GAMA2*AL2)+YL22*ETA22*SINH(GAMA2*AL2)
	DELT1=BAMT1*BAMT4-BAMT2*BAMT3
	BMT1=COSH(GAMA1*AL1)+YL11*ETA11*SINH(GAMA1*AL1)
	BMT2=YL12*ETA22*SINH(GAMA2*AL2)
	BMT3=YL21*ETA11*SINH(GAMA1*AL1)
	BMT4=COSH(GAMA2*AL2)+YL22*ETA22*SINH(GAMA2*AL2)
	DELT2=BMT1*BMT4-BMT2*BMT3
	DMAT11=YL11*COSH(GAMA1*AL1)+SINH(GAMA1*AL1)/ETA11
	DMAT12=YL12*ETA22*SINH(GAMA2*AL2)
	DMAT13=YL21*COSH(GAMA1*AL1)
	DMAT14=COSH(GAMA2*AL2)+YL22*ETA22*SINH(GAMA2*AL2)
	DMAT1=DMAT11*DMAT14-DMAT12*DMAT13
	DMBT11=YL12*COSH(GAMA2*AL2)
	DMBT12=YL12*ETA22*SINH(GAMA2*AL2)
	DMBT13=SINH(GAMA2*AL2)/ETA22+YL22*COSH(GAMA2*AL2)
	DMBT14=COSH(GAMA2*AL2)+YL22*ETA22*SINH(GAMA2*AL2)
	DMBT1=DMBT11*DMBT14-DMBT12*DMBT13
	DMAT21=COSH(GAMA1*AL1)+YL11*ETA22*SINH(GAMA1*AL1)
	DMAT22=YL11*COSH(GAMA1*AL1)+SINH(GAMA1*AL1)/ETA11
	DMAT23=ETA11*YL21*SINH(GAMA1*AL1)
	DMAT24=YL21*COSH(GAMA1*AL1)
	DMAT2=DMAT21*DMAT24-DMAT22*DMAT23
	DMBT21=COSH(GAMA1*AL1)+YL11*ETA11*SINH(GAMA1*AL1)
	DMBT22=YL12*COSH(GAMA2*AL2)
	DMBT23=YL21*ETA11*SINH(GAMA1*AL1)
	DMBT24=SINH(GAMA2*AL2)/ETA22+YL22*COSH(GAMA2*AL2)
	DMBT2=DMBT21*DMBT24-DMBT22*DMBT23
	CMAT(1,1)=DMAT1/DELT1
	CMAT(1,2)=DMBT1/DELT2
	CMAT(2,1)=DMAT2/DELT1
	CMAT(2,2)=DMBT2/DELT2
 1700	RETURN
	END
C
	SUBROUTINE HF(K,WL,XX,YY,AA,BB,HW,HD1,SMAG)
C this is for calculate scattering field strength in space.
	IMPLICIT NONE
	COMPLEX*8 CVAL,FIELD,SMAG
	REAL*4 WL,XX,YY,AA,BB,HW,HD1,XX0,YY0,WLT,HHW,HHD,AB
	INTEGER K,KK
	COMMON/T13/ KK
	COMMON/T3/ XX0,YY0,WLT,HHW,HHD,AB
	EXTERNAL FIELD	
	XX0=XX
	YY0=YY
	WLT=WL
	HHW=HW
	HHD=HD1
	AB=AA
	KK=K
	CALL CGQ1(FIELD,AA,BB,10,CVAL)
	SMAG=CVAL
	RETURN
	END
C
	COMPLEX*8 FUNCTION FIELD(X)
	IMPLICIT NONE
	COMPLEX*8 H20
	REAL*4 XX0,YY0,WLT,HHW,HHD,AB,SGN,RR,X1,X
	INTEGER KK
	COMMON/T13/ KK
	COMMON/T3/ XX0,YY0,WLT,HHW,HHD,AB
	EXTERNAL H20,SGN
	RR=2.*3.141592654/WLT*SQRT((X-XX0)**2+YY0*YY0)
	IF(KK.EQ.1) GOTO 200
	X1=KK*3.141592654*(X-HHD*SGN(AB))/HHW
	FIELD=H20(RR)*COS(X1)
	GOTO 100
 200	FIELD=H20(RR)
 100	RETURN
	END
C
	SUBROUTINE AFHF(K,WL,RR,ANG,AA,BB,HW,HD1,FMAG)
c this subroutine is for calculating the approximate scattering field strength.
	IMPLICIT NONE
	COMPLEX*8 APROX,ZJ,ZW,ZEFF,ZK,ZEFFE,CVAL,FMAG
	REAL*4 WL,RR,ANG,AA,BB,HD1,WLT,ANGT,HW,HWT,HDT,RRT
	INTEGER K,KK
	COMMON/T5/ WLT,ANGT,HWT,HDT,RRT
	COMMON/T51/ KK
	EXTERNAL APROX
	WLT=WL
	ANGT=ANG
	HWT=HW
	HDT=HD1
	RRT=RR
	KK=K
	ZJ=(0.,1.)
	ZW=ZJ*WL/RR
	ZEFF=SQRT(ZW)/3.141592654
	ZK=-ZJ*2.*3.141592654/WL*RR
	ZEFFE=EXP(ZK)
	CALL CGQ1(APROX,AA,BB,10,CVAL)
	FMAG=CVAL*ZEFF*ZEFFE
	RETURN
	END
C
	COMPLEX*8 FUNCTION APROX(X)
	IMPLICIT NONE
	COMPLEX*8 ZJ,ZX,ZZX,ZM,ZZM
	REAL*4 WLT,ANGT,HWT,HDT,RRT,X
	INTEGER KK
	COMMON/T5/ WLT,ANGT,HWT,HDT,RRT
	COMMON/T51/ KK
	ZJ=(0.,1.)
	ZX=ZJ*2.*3.141592654*X/WLT*COS(ANGT)
	ZZX=EXP(ZX)
	IF(KK.EQ.1) GOTO 101
	ZM=KK*3.141592654*(X-HDT)/HWT
 	ZZM=COS(ZM)
	APROX=ZZX*ZZM
	GOTO 202
 101	APROX=ZZX
 202	RETURN
	END
C
C----------SUBROUTINE--CGQ1-------------------------------------------
C                                                                     
      SUBROUTINE CGQ1(CF,XL,XU,N,CVAL)                                
C*********************************************************************
C
C
C        PERFORMS INTEGRATION OF A FUNCTION OF A SINGLE VARIABLE
C        BY GAUSSIAN QUADRATURE.
C
C        N = ORDER OF GUASSIAN QUADRATURE APPROXIMATION
C            (2,4,8,10,12,16,32)
C        CF = EXTERNALLY SUPPLIED FUNCTION....MUST BE FUNCTION OF ONE
C             VARIABLE FOR CQG1.
C        XL = LOWER BOUND OF VARIABLE
C        XU = UPPER BOUND OF VARIABLE
C        CVAL = RESULTING VALUE OF THE INTEGRATION
C
C        CF MUST BE LISTED IN AN EXTERNAL STATEMENT
C
C       PREPARED BY MICHAEL G. HARRISON  E.E. DEPT    JUNE 22, 1972
C
C
C***********************************************************************
 
 
      IMPLICIT NONE
	COMPLEX*8 CF,CVAL
	INTEGER N,L,NP,NE,J,NS(8),NQ(8)
	REAL*4 AX,BX,DX,XL,XU
      REAL*4 Q1(52),Q2(24),Q3(32),QG(108)
      EQUIVALENCE (Q1(1),QG(1)),(Q2(1),QG(53)),(Q3(1),QG(77))
	EXTERNAL CF
      DATA Q1
     $               /.288675134594812882,0.50,.43056815579702629,
     $.17392742256872693,.16999052179242813,.32607257743127307,  
     $ 0.48014492824876812,.50614268145188130E-1,.39833323870681337,
     $.11119051722668724,.26276620495816449,.15685332293894364,
     $.9171732124782490E-1,.18134189168918099,.48695326425858586,
     $.3333567215434407E-1,.43253168334449225,.747256745752903E-1,
     $.3397047841496122,.10954318125799102,.2166976970646236,
     $.13463335965499818,.74437169490815605E-1,.14776211235737644,
     $0.49078031712335963,.23587668193255914E-1,.45205862818523743,
     $.53469662997659215E-1,.38495133709715234,.8003916427167311E-1,
     $.29365897714330872,.10158371336153296,.18391574949909010,
     $.11674626826917740,.62616704255734458E-1,.12457352290670139,
     $.49470046749582497,.13576229705877047E-1,.47228751153661629,
     $.31126761969323946E-1,.43281560119391587,.47579255841246392E-1,
     $.37770220417750152,.62314485627766936E-1,.30893812220132187,
     $.7479799440828837E-1,.22900838882861369,.8457825969750127E-1,
     $.14080177538962946,.9130170752246179E-1,.47506254918818720E-1,
     $.9472530522753425E-1 /
      DATA Q2
     $              / 0.49759360999851068   ,0.61706148999935998E-2 ,
     *            0.48736427798565475   ,    0.14265694314466832E-1 ,
     *            0.46913727600136638   ,    0.22138719408709903E-1, 
     *            0.44320776350220052   ,    0.29649292457718890E-1 ,
     *            0.41000099298695146   ,    0.36673240705540153E-1 ,
     *            0.37006209578927718   ,    0.43095080765976638E-1 ,
     *            0.32404682596848778   ,    0.48809326052056944E-1 ,
     *            0.27271073569441977   ,    0.53722135057982817E-1 ,
     *            0.21689675381302257   ,    0.57752834026862801E-1 ,
     *            0.15752133984808169   ,    0.60835236463901696E-1 ,
     *            0.95559433736808150E-1   ,    0.62918728173414148E-1 ,
     *            0.32028446431302813E-1   ,    0.63969097673376078E-1 /
      DATA Q3
     $              /0.49863193092474078 ,   0.35093050047350483E-2 ,
     *            0.49280575577263417   ,    0.81371973654528350E-2 ,
     *            0.48238112779375322   ,    0.12696032654631030E-1 ,
     *            0.46745303796886984   ,    0.17136931456510717E-1 ,
     *            0.44816057788302606   ,    0.21417949011113340E-1 ,
     *            0.42468380686628499   ,    0.25499029631188088E-1 ,
     *            0.39724189798397120   ,    0.29342046739267774E-1 ,
     *            0.36609105937014484   ,    0.32911111388180923E-1 ,
     *            0.33152213346510760   ,    0.36172897054424253E-1 ,
     *            0.29385787862038116   ,    0.39096947893535153E-1 ,
     *            0.25344995446611470   ,    0.41655962113473378E-1 ,
     *            0.21067563806531767   ,    0.43826046502201906E-1 ,
     *            0.16593430114106382   ,    0.45586939347881942E-1 ,
     *            0.11964368112606854   ,    0.46922199540402283E-1 ,
     *            0.72235980791398250E-1   ,    0.47819360039637430E-1 ,
     *            0.24153832843869158E-1   ,    0.48270044257363900E-1 /
     $,NQ/2,4,8,10,12,16,24,32/,
     $ NS/1,3,7,15,25,37,53,77/ 
      IF(N .EQ. 1) GO TO 200
      DO 300 L=1,8
      IF(N.EQ.NQ(L)) GO TO 301
 300  CONTINUE 
      WRITE(6,905) N
 905  FORMAT('0 CALLING PARAMETER =',I5,'   INTEGRATION NOT POSSIBLE'//)
      RETURN
 200  AX=0.5*(XU+XL)
      CVAL=CF(AX)*CMPLX(XU-XL,0.0)
      RETURN
 301  CONTINUE
      NP=NS(L)
      NE=NP+N-1
      AX=0.5*(XU+XL)
      BX=(XU-XL)
      CVAL=(0.0,0.0)
      DO 350 J=NP,NE,2
      DX=QG(J)*BX
      CVAL=CVAL+CMPLX(QG(J+1),0.0)*(CF(AX+DX)+CF(AX-DX))
 350  CONTINUE 
      CVAL=CVAL*CMPLX(BX,0.0)
      RETURN 
      END 
