	pROGRAM BPAT 
c this is a program for calculating two apertures on an infinite
c extention conducting plane backed waveguides and network.  
	IMPLICIT NONE  
	INTEGER N1,N2  
	REAL UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,
     1	R0,WL,ANG0,FAI
	COMPLEX, ARRAY(2,2) :: YL       
	INTEGER I,J,M,N 
	COMMON/T0/ R0,ANG0   
        COMMON/T2/ UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,WL,FAI  
	PRINT*,'We use implse based mode function expansion in the 
     1  moment method. Please input the expansion parameters.'   
	PRINT*,'I and J are the number of pulse base mode functions,'  
	PRINT*,'please input I and J.'  
	READ(*,*) I,J 
	PRINT*,'I=', I,'J=',J  
	PRINT*,'M and N are the number of pulse base functions, please
     1	input M and N.'  
	READ(*,*) M,N  
	PRINT*,'M=',M,'N=',N  
	PRINT*,'N1 is the sum of I and J.'  
	N1=I+J  
	PRINT*,'N1=',N1 
	PRINT*,'For getting scattering field strength pattern, we 
     1  divide the space by eaqual angals. The total measurement 
     2  point number is N2. Please input N2.'   
	READ(*,*) N2  
	PRINT*,'N2=',N2  
	PRINT*,'The following are the parameters of the two apertures  
     1  and the material medium. Please input the values one by one'  
	PRINT*,'aperture widths, input a, then input b(in meter) '  
	READ(*,*) HA, HB  
	PRINT*,'a=',HA, 'b=',HB  
	PRINT*,'input aperture separation d(in meter)'  
	READ(*,*) HD  
	PRINT*,'d=',HD  
	PRINT*,'input lengths of the two parallel plate waveguide(in   
     1  meter)'  
	READ(*,*) AL1,AL2  
	PRINT*,'L1=',AL1,'L2=',AL2  
	PRINT*,'input material medium parameters according to this order
     1  (1) space permittivity and (2) permeability. Then (3) wavequide
     2  permittivity and (4)  permeability.'    
	READ(*,*) EM1,UM1,EM2,UM2  
	PRINT*,'EM1=',EM1,'UM1=',UM1,'EM2=',EM2,'UM2=',UM2  
	PRINT*,'Input the electromagnetic wave incident angle(in
     1  radian) and wave-length(in meter)' 
	READ(*,*) FAI,WL   
	PRINT*,'Fi=',FAI,'WL=',WL    
	PRINT*,'Input the distance value from measurement points  
     1  to the middle of the two apertures(all point share the same 
     2  radiate distance). Then input the begin angal.(meter & radian)'  
	READ(*,*) R0,ANG0  
	PRINT*,'r=',R0,'The begin point angal=',ANG0  
	PRINT*,'Input the admittances of the microwave network.'  
	PRINT*,'input Y11,Y12,Y21,Y22 one by one(in 1/ohm).'  
	READ(*,*) YL(1,1),YL(1,2),YL(2,1),YL(2,2)    
	PRINT*,'Y11=',YL(1,1),'Y12=',YL(1,2),  
     1  'Y21=',YL(2,1),'Y22=',YL(2,2)   
	PRINT*,'Ok! You finish inputing oraginal data, afterwards you
     1  can pick up the results in the files: "bpouts1.dat" and "
     2  "bpouts52.dat".' 
	CALL MAIN(N1,N2,I,J,M,N,YL)    
	STOP  
	END  
c 
	SUBROUTINE MAIN(N1,N2,IS,JS,MS,NS,YL)  
	IMPLICIT NONE  
	INTEGER N1,N2  
	COMPLEX YA(N1,N1)  
	COMPLEX, ARRAY(N2,N1) :: SMAG,FMAG  
	COMPLEX, ARRAY(N2,N1) :: SHM,FHM
	COMPLEX, ARRAY(N2) :: HFIELD,FARH,HSCRR,
     1  AG,HMAG,HSC,HSCN,HTAT,HSCT
	COMPLEX, ARRAY(N1) :: AIJ,COEF,VM  
	COMPLEX, ARRAY(N1,1) :: AIJS   
	COMPLEX, ARRAY(N1,N1) :: CMTT,YB,YIJ,YTOT,YII,YJI,YJJ   
	REAL, ARRAY(N1,N1) :: YBE   
	REAL, ARRAY(N1*N1) :: YE   
	COMPLEX, ARRAY(2,2) :: YL,CMT       
        COMPLEX YMI,YMJ,ZJ  
	REAL, ARRAY(N1) :: AIR,AII,AIA,AVR,AVA,AVI
        REAL, ARRAY(N2) :: HMGR,HMGI,HMGA,HTTR,HTTI
     1  ,HSCR,HSCI,X0,Y0,AGA,AA1,AA2,BB1,BB2,HSNR,HTTA,HSCA,ANG
	REAL UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,TIME,
     1	R0,WL,CSQ,DANG,ANG0,FAI
	INTEGER, ARRAY(N2) :: LA
	INTEGER, ARRAY(N1) :: IJA
	INTEGER JJ,II,III,IERF,IERS,IER,ITIM,mpTimerElapsed          
	INTEGER NIJ,IJ,I,J,M,N,L,IS,JS,MS,NS,IN,ID    
	COMMON/T0/ R0,ANG0   
	COMMON/T1/ L,NIJ,I,J,M,N 
        COMMON/T2/ UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,WL,FAI
C OPEN em9out.dat ON 11 FOR OUTPUT  
	CALL mpTimerStart()       !start timer 0 
	OPEN (11,FILE='em9out1.dat',STATUS='UNKNOWN')
	OPEN (49,FILE='em9out2.dat',STATUS='UNKNOWN')
C READ IN ORIGINAL DATA.
	I=IS  
	J=JS  
	M=MS  
	N=NS   
	L=N2  
	NIJ=N1  
	YA=(0.,0.)  
	SHM=(0.,0.)  
	FHM=(0.,0.)  
	HFIELD=(0.,0.)  
	FARH=(0.,0.)  
	HSCRR=(0.,0.)  
	AG=(0.,0.)  
	HMAG=(0.,0.)  
	HSC=(0.,0.)  
	HSCN=(0.,0.)  
	HTAT=(0.,0.)  
	HSCT=(0.,0.)  
	AIJ=(0.,0.)  
	COEF=(0.,0.)  
	VM=(0.,0.)  
	CMTT=(0.,0.)  
	YB=(0.,0.) 
	YII=(0.,0.)  
	YJI=(0.,0.)  
	YJJ=(0.,0.) 
	YIJ=(0.,0.)  
	YTOT=(0.,0.)  
	AIR=0.  
	AII=0.  
	AIA=0.  
	AVR=0.  
	AVA=0.  
	AVI=0.  
	HMGR=0.  
	HMGI=0.  
	HMGA=0.  
	HTTR=0.  
	HTTI=0.  
	HSCR=0.  
	HSCI=0.  
	X0=0.  
	Y0=0.  
	AGA=0.  
	AA1=0.  
	AA2=0.  
	BB1=0.  
	BB2=0.  
	HSNR=0.  
	HTTA=0.  
	HSCA=0.  
	ANG=0. 
	IJ=I+J
	LA=[1:L]
	IJA=[1:NIJ]
c calculate the equivalent microwave network parameters
	ZJ=(0.,1.) 
	IN=IJ*IJ  
	ID=IJ+1  
	YE=0.   
	YMI=1/SQRT(UM1/EM1)/377*ZJ
	YMJ=1/SQRT(UM2/EM2)/377*ZJ
	FORALL(II=1:IN:ID) YE(II)=1. 
	YBE=RESHAPE(SOURCE=YE,SHAPE=[IJ,IJ])   
	FORALL(II=2:I,JJ=2:I) YB(II,JJ)=YMI/SQRT(1.-(( 
     1   HA/((II-1)*WL)*2.)**2)*(UM1*EM1))  
	FORALL(II=I+2:IJ,JJ=I+2:IJ) YB(II,JJ)=YMJ/SQRT(1.-
     1  ((HB/((II-I-1)*WL)*2.)**2)*(UM2*EM2))
	YB=YB*YBE  
	CALL CNETWORK(YL,CMT)
	CMTT(1:I+1:I,1:I+1:I)=CMT   
	YB=YB+CMTT 
	WRITE(11,*) 'Yb is the admittance of the network'   
	WRITE(11,*) YB(1:IJ,1:IJ)  
	CALL APULSE(I,I,M,M,HA,HA,YA,HD,HD)  
	YII(1:I,1:I)=YA(1:I,1:I)  
	CALL APULSE(I,J,M,N,HA,HB,YA,HD,-HD)
	YIJ(1:I,I+1:IJ)=YA(1:I,1:J)  
	CALL APULSE(J,I,N,M,HB,HA,YA,-HD,HD)  
	YJI(I+1:IJ,1:I)=YA(1:J,1:I)  
	CALL APULSE(J,J,N,N,HB,HB,YA,-HD,-HD)
	YJJ(I+1:IJ,I+1:IJ)=YA(1:J,1:J)  
	YTOT=YII+YJI+YJJ+YIJ+YB  
	WRITE(11,*) 'this is total admittances'  
	WRITE(11,*) YTOT(1:IJ,1:IJ)    
	CALL CURRENT(AIJ)   
	AIA=ABS(AIJ)
	WRITE(11,*) 'the following are current values'
	WRITE(11,*) AIJ(1:IJ)      
	CALL GAUSS(IJ,YTOT(1:IJ,1:IJ),AIJ(1:IJ),VM)  
	WRITE(11,*) 'the following are volages values'
	WRITE(11,*) VM(1:IJ)   
	AVA=ABS(VM)
	AA1=HD
	BB1=HD+HA
	AA2=-HD-HB
	BB2=-HD
	COEF(1:I)=VM(1:I)*SQRT(2./HA)
	COEF(I+1:IJ)=VM(I+1:IJ)*SQRT(2./HB) 
	CSQ=1./SQRT(2.) 
	WHERE((IJA.EQ.1).OR.(IJA.EQ.I+1)) 
	COEF=COEF*CSQ 
	END WHERE  
	DANG=3.141592654/(L-1)   !when L=901 DANG=3.49065852E-3   
	ANG=(LA-1)*DANG+ANG0
	X0=R0*COS(ANG)
	Y0=R0*SIN(ANG)
	CALL HF(I,WL,X0,Y0,AA1,BB1,HA,HD,SMAG)
	SHM(:,1:I)=SMAG(:,1:I)  
	CALL HF(J,WL,X0,Y0,AA2,BB2,HB,-HB-HD,SMAG)
	SHM(:,I+1:IJ)=SMAG(:,1:J)   
	HFIELD=MATMUL(SHM(:,1:IJ),COEF(1:IJ))   
	CALL AFHF(I,WL,R0,ANG,AA1,BB1,HA,HD,FMAG)
	FHM(:,1:I)=FMAG(:,1:I)
	CALL AFHF(J,WL,R0,ANG,AA2,BB2,HB,-HB-HD,FMAG)
	FHM(:,I+1:IJ)=FMAG(:,1:J)
	FARH=MATMUL(FHM(:,1:IJ),COEF(1:IJ))   
	HMAG=2.*3.141592654*3.E8/WL*8.841941283E-12*HFIELD
	AG=2.*3.141592654*3.E8/WL*8.841941283E-12*FARH
	HSCRR=-ZJ*2.*3.141592654/WL*(X0*COS(FAI)+Y0*SIN(FAI))
	HSC=EXP(HSCRR)
	HSCN=EXP(-1.*HSCRR)
	HSNR=ABS(HSCN)
	HTAT=HMAG+HSC+HSCN
	HSCT=HSCN+HMAG
	HMGR=REAL(HMAG)
	HMGI=AIMAG(HMAG)
	HMGA=ABS(HMAG)
	HTTR=REAL(HTAT)
	HTTI=AIMAG(HTAT)
	HTTA=ABS(HTAT)
	HSCR=REAL(HSCT)
	HSCI=AIMAG(HSCT)
	HSCA=ABS(HSCT)
	AGA=ABS(AG)  
	WRITE(49,33) HMGA(1:L) 
33	FORMAT(1x,F12.6)   
	ITIM=mpTimerElapsed()
	print*,ITIM,'milliseconds'  
	TIME=ITIM/1000.   
	PRINT*,'time=',TIME,'secondes'  
	RETURN   
	END    
C 
C
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 N,NC
        COMPLEX*8 UL(N,N),PIVOT,EM
      REAL*4 SCL(N)
        INTEGER IPS(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,ULF,BF,XF)
        IMPLICIT NONE
        INTEGER N,NC
      COMPLEX*8 UL(N,N),B(N),X(N),SUM
	COMPLEX, ARRAY(N) :: BF,XF  
	COMPLEX, ARRAY(N,N) :: ULF  
      INTEGER IPS(N)
        INTEGER IP,I,J,IM1,IBAOK,IP1,NP1
	ULF=ULF  
	BF=BF 
	XF=XF 
	DO 102 I=1,N  
	DO 101 J=1,N  
	UL(I,J)=ULF(I,J)  
101	CONTINUE  
102	CONTINUE  
	DO 103 I=1,N  
	B(I)=BF(I)  
103	CONTINUE  
      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)
	DO 104 I=1,N  
	XF(I)=X(I)  
104	CONTINUE 
	XF=XF   
      RETURN
      END
C 
C 
C
C
	SUBROUTINE CURRENT(AIJ)
C this is a subroutine for calculating the excitation vector
	IMPLICIT NONE 
	INTEGER NIJ,II,IJ,I,J,M,N,L 
	COMMON/T1/ L,NIJ,I,J,M,N 
        COMMON/T2/ UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,WL,FAI
	COMPLEX, ARRAY(NIJ) :: AIJ    
	COMPLEX, ARRAY(I) :: AII
	COMPLEX, ARRAY(J) :: BII  
	COMPLEX, ARRAY(I,M) :: AID,AX
	COMPLEX, ARRAY(J,N) :: BID,BX
	COMPLEX ZJ
	INTEGER, ARRAY(I) :: KI
	INTEGER, ARRAY(J) :: KJ
	INTEGER, ARRAY(M) :: KM
	INTEGER, ARRAY(N) :: KN
	INTEGER, ARRAY(I,M) :: KMI,KIM
	INTEGER, ARRAY(J,N) :: KNJ,KJN
	REAL, ARRAY(I,M) :: RXM,CXKI,AIM,AAX
	REAL, ARRAY(J,N) :: RXN,CXKJ,BIM,BBX
	REAL UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,WL,FAI,
     1  DELT1,DELT2,HK,AKM,AKN,TES,EFC1,EFC2
	ZJ=(0.,1.)  
	AIJ=(0.,0.)  
	IJ=I+J   
	AID=(0.,0.)  
	BID=(0.,0.)  
	AX=(0.,0.)  
	BX=(0.,0.)  
	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-8) THEN  
	EFC1=2*DELT1
	EFC2=2*DELT2
	ELSE  
	EFC1=2.*SIN(AKM*COS(FAI))/(HK*COS(FAI))
	EFC2=2.*SIN(AKN*COS(FAI))/(HK*COS(FAI))
	END IF 
	KI=[1:I]     ! KI=[1:I] is faster than FORALL(II=1:I) KI(II)=II
	KJ=[1:J]
	KM=[1:M]
	KN=[1:N] 
	KIM=SPREAD(KI,DIM=2,NCOPIES=M)
	KMI=SPREAD(KM,DIM=1,NCOPIES=I)
	KJN=SPREAD(KJ,DIM=2,NCOPIES=N)
	KNJ=SPREAD(KN,DIM=1,NCOPIES=J)
	RXM=(KMI+0.5)*DELT1
	CXKI=KIM*RXM*3.141592654/HA  
	WHERE(KIM.EQ.1) 
	AIM=1./SQRT(HA)
	ELSEWHERE
	AIM=SQRT(2./HA)*COS(CXKI)
	END WHERE
	AAX=HK*(RXM+HD)*COS(FAI)
	AX=-ZJ*AAX
	WHERE (AAX.LE.-78) 
	AID=0
	ELSEWHERE
	AID=AIM*EXP(AX)*EFC1
	END WHERE
	AII=SUM(AID,DIM=2)
	RXN=-(KNJ+0.5)*DELT2
        CXKJ=KJN*RXN*3.141592654/HB
        WHERE(KJN.EQ.1) 
	BIM=1./SQRT(HB)
        ELSEWHERE
        BIM=SQRT(2./HB)*COS(CXKJ)
        END WHERE
        BBX=HK*(RXN-HD)*COS(FAI)
        BX=ZJ*BBX
        WHERE (BBX.LE.-78) 
	BID=0
        ELSEWHERE
        BID=BIM*EXP(BX)*EFC2
        END WHERE
        BII=SUM(BID,DIM=2)  
	AIJ(1:I)=AII(1:I)  
	AIJ(I+1:IJ)=BII(1:J)  
	RETURN 
	END
C
C
	SUBROUTINE APULSE(KI,KJ,KM,KN,HW1,HW2,YA,START1,START2)
C this is a subroutine for calculating the  eguivalent admittance of 'A' 
C region by using the pulse basis mode function expansion.  
	IMPLICIT NONE
	INTEGER L,II,NIJ,JJ,I,J,M,N,KI,KJ,KN,KM,IER  
	COMMON/T1/ L,NIJ,I,J,M,N
        COMMON/T2/ UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,WL,FAI
	INTEGER, ARRAY(KI) :: KII
	INTEGER, ARRAY(KJ) :: KJJ
	INTEGER, ARRAY(KM) :: KMM
	INTEGER, ARRAY(KN) :: KNN
	INTEGER, ARRAY(KI,KM) :: KIA
	INTEGER, ARRAY(KN,KJ) :: KJA  
	INTEGER, ARRAY(NIJ,KI) :: ATI  
	COMPLEX, ARRAY(NIJ,KJ) :: BTI  
	REAL, ARRAY(KM) :: A,B,XM
	REAL, ARRAY(KN) :: C,D,XN  
	COMPLEX, ARRAY(KM,KJ) :: YMJ 
	COMPLEX, ARRAY(KI,KM) :: YIZ  
	COMPLEX, ARRAY(KN,KJ) :: YJZ   
	REAL, ARRAY(KM,KN) :: YS,XF,YF,ZF 
	REAL, ARRAY(KI,KM) :: YI,XMA
	REAL, ARRAY(KN,KJ) :: YJ,XNA
	COMPLEX, ARRAY(KM,KN) :: SOUT,SOUT1
	COMPLEX, ARRAY(KI,KJ) :: YA1
	COMPLEX, ARRAY(NIJ,NIJ) :: YA 
	COMPLEX UT   
	REAL YT,HW1,HW2,GMN,UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,
     1  WL,FAI,SGN,START1,START2
	EXTERNAL SGN  
	A=0. 
	B=0.  
	C=0.  
	D=0.  
	YIZ=(0.,0.)  
	YJZ=(0.,0.)   
	BTI=(0.,0.) 
	ATI=0  
	YI=0. 
	XMA=0.  
	XNA=0.  
	YJ=0.  
	XN=0.  
	XM=0.   
	YS=0.  
	XF=0.  
	YF=0.  
	ZF=0. 
	SOUT=(0.,0.)  
	SOUT1=(0.,0.)  
	YMJ=(0.,0.)  
	YA1=(0.,0.)   
	YA=(0.,0.) 
	YMJ=YMJ  
	SOUT=SOUT  
	SOUT1=SOUT1  
	YA1=YA1  
	YA=YA   
	KII=[1:KI] 
	KJJ=[1:KJ]
	KMM=[1:KM]
	KNN=[1:KN]
	KIA=SPREAD(KII,DIM=2,NCOPIES=KM)
	KJA=SPREAD(KJJ,DIM=1,NCOPIES=KN)
	XM=(KMM+0.5)*HW1/KM*SGN(START1)
	XN=(KNN+0.5)*HW2/KN*SGN(START2)
	XMA=SPREAD(XM,DIM=1,NCOPIES=KI)
	XNA=SPREAD(XN,DIM=2,NCOPIES=KJ)
	WHERE(KIA.EQ.1)
	YI=1./SQRT(HW1)
	ELSEWHERE
	YI=SQRT(2./HW1)*COS(3.141592654*(KIA-1)*XMA/HW1)
	END WHERE
	WHERE(KJA.EQ.1)
	YJ=1./SQRT(HW2)
	ELSEWHERE
	YJ=SQRT(2./HW2)*COS(3.141592654*(KJA-1)*XNA/HW2) 
	END WHERE
	A=START1+(KMM-1)*HW1/KM*SGN(START1)
	B=START1+(KMM)*HW1/KM*SGN(START1)
	C=START2+(KNN-1)*HW2/KN*SGN(START2)
	D=START2+(KNN)*HW2/KN*SGN(START2) 
	GMN=SQRT(UM1*EM1)*2.*3.141592654/WL*HW1/KN*1.63804957E-1
	YT=6.366197725E-1*ALOG(GMN)  
	UT=CMPLX(1.,YT)*HW1*HW2/KM/KN  
	CALL CGAUS(KM,KN,A,B,C,D,WL,UT,SOUT)
	SOUT1=SOUT*SGN(START1)*SGN(START2)  
	YIZ=CMPLX(YI,0.)  
	YJZ=CMPLX(YJ,0.)  
	YMJ=MATMUL(SOUT1,YJZ)   
	YA1=MATMUL(YIZ,YMJ)   
	YA(1:KI,1:KJ)=YA1(1:KI,1:KJ)/WL*1.66663E-2   
	RETURN
	END
C
C
      SUBROUTINE CGAUS(M,N,A,B,C,D,WL,TR,SOUT)
      IMPLICIT NONE      
	INTEGER I,J,M,N
	COMPLEX, ARRAY(M,N) :: SS,ZZ,SOU,SOUT
	REAL, ARRAY(M) :: A,B
	REAL, ARRAY(N) :: C,D
	REAL, ARRAY(M,N) :: DELTAX,DELTAZ,AL,BL,EL,A4  
	REAL, ARRAY(M,N) :: CL,DL,XM,XR,ZM,ZR,A1,A2,A3
	REAL, ARRAY(5) :: X,W
	REAL WL
	COMPLEX, ARRAY(M,N) :: GN,GN1,GN2,GN3,GN4  
	COMPLEX TR   
      DATA X/1.488743389E-1,4.333953941E-1,
     1          6.794095682E-1,8.60633666E-1,9.739065285E-1/
      DATA W/2.955242247E-1,2.692667193E-1,
     2          2.190863625E-1,1.494513491E-1,6.66713443E-2/
	AL=0.  
	BL=0. 
	EL=0.   
	CL=0.  
	DL=0.  
	XM=0.  
	XR=0.  
	ZM=0.  
	ZR=0.  
	DELTAX=0.  
	DELTAZ=0. 
	A=A  
	B=B  
	C=C  
	D=D  
	AL=SPREAD(A,DIM=2,NCOPIES=N)
	BL=SPREAD(B,DIM=2,NCOPIES=N)
	CL=SPREAD(C,DIM=1,NCOPIES=M)
	DL=SPREAD(D,DIM=1,NCOPIES=M)
	EL=ABS(AL-CL) 
	XM=0.5*(AL+BL)
	XR=0.5*(BL-AL)	
	ZM=0.5*(CL+DL)
	ZR=0.5*(DL-CL)  
	SOUT=(0.,0.) 
	SOU=(0.,0.)    
        SS=CMPLX(0.0,0.0) 
	GN=(0.,0.)  
	GN1=(0.,0.)  
	GN2=(0.,0.)  
	GN3=(0.,0.)  
	GN4=(0.,0.)  
	A1=0.
	A2=0.
	A3=0.
	A4=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)
        CALL H20(A1,M,N,GN) 
        GN1=GN    
        CALL H20(A2,M,N,GN) 
        GN2=GN 
        CALL H20(A3,M,N,GN)  
        GN3=GN  
        CALL H20(A4,M,N,GN)  
        GN4=GN 
        ZZ=ZZ+W(I)*(GN1+GN2+GN3+GN4)  
12      CONTINUE
	SS=SS+W(J)*ZZ*XR
11    CONTINUE
      SOU=ZR*SS
	WHERE(EL.LE.1.E-8) 
	SOUT=TR  
	ELSEWHERE 
	SOUT=SOU   
	END WHERE  
      RETURN
      END
C
C
        SUBROUTINE H20(A,IM,IN,HMN)
	IMPLICIT NONE    
	INTEGER IM,IN
	REAL, ARRAY(IM,IN) :: BJ,BY,F0,T0,A 
	COMPLEX, ARRAY(IM,IN) :: ZO  
	COMPLEX, ARRAY(IM,IN) :: HMN  
	ZO=(0.,0.) 
	BJ=0.  
	BY=0.  
	F0=0.  
	T0=0.  
	HMN=(0.,0.) 
	A=A  
	ZO=ZO   
	WHERE(A.LE.1.E-8)  
	HMN=ZO  
	END WHERE 
	WHERE((A.GT.0.).and.(A.LE.3.)) 	 
        BJ=A*A/9.  
        BJ=1.0+BJ*(-2.2499997+BJ*(1.2656208+BJ*
     1  (-3.163866E-1+BJ*(4.44479E-2+BJ*(-3.9444E-3+BJ*
     2  2.1E-4)))))
        BY=A*A/9.  
        BY=2.0/3.14159265*ALOG(A/2.0)*BJ+3.674669E-1
     1  +BY*(6.0559366E-1+BY*(-7.4350384E-1+BY*(2.5300117E-1+BY*(
     2  -4.261214E-2+BY*(4.27916E-3-BY*2.4846E-4)))))  
	HMN=CMPLX(BJ,-BY)  
	END WHERE 
	WHERE(A.GT.3.) 
        BJ=3.0/A
        F0=7.9788456E-1+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-7.8539816E-1+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
        HMN=CMPLX(BJ,-BY)
	END WHERE  
        RETURN
        END
C
C
	REAL FUNCTION SGN(X)
	IMPLICIT NONE 
	REAL 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,CMT)
	IMPLICIT NONE
	COMPLEX, ARRAY(2,2) :: YL,CMAT,CMT  
	COMPLEX SINH,COSH,YL11,YL12,YL21,YL22,Z
     1  ,GA1,GAMA1,GA2,GAMA2,BAMT1,BAMT2,BAMT3,BAMT4,DELT1,DELT2,
     2  BMT1,BMT2,BMT3,BMT4,DMAT11,DMAT12,DMAT13,DMAT14,DMAT1,DMBT11,
     3  DMBT12,DMBT13,DMBT14,DMBT1,DMAT21,DMAT22,DMAT23,DMAT24,DMAT2,
     4  DMBT21,DMBT22,DMBT23,DMBT24,DMBT2,ZJ
	REAL UM1,UM2,EM1,EM2,HA,HB,HD,AL1,AL2,WL,FAI,HK1,
     1  HK2,ETA11,ETA22
	INTEGER L,IJ,I,J,M,N,NN  
	COMMON/T1/ L,IJ,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))
	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
	ZJ=(0.,1.)
	NN=0  ! This is (transmission mode number minuse 1) 
 	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 
	CMT=CMT   
	CMT=CMAT  
	RETURN
	END
C
	SUBROUTINE HF(K,WL,XX,YY,AA,BB,HW,HD1,SMAG)
C THIS IS FOR CALCULATE SCATTERING FIELD IN SPACE.
	IMPLICIT NONE
	INTEGER L,II,NIJ,K,I,J,M,N
	COMMON/T1/ L,NIJ,I,J,M,N
	COMPLEX, ARRAY(L,K) :: CVAL
	COMPLEX, ARRAY(L,NIJ) :: SMAG
	REAL, ARRAY(L) :: AA,BB,XX,YY   ! 901=L 
	REAL WL,HW,HD1,WLT,HHW,HHD,AB
	INTEGER LK,KK
	COMMON/T13/ LK,KK
	COMMON/T3/ WLT,HHW,HHD,AB
	WLT=WL
	HHW=HW
	HHD=HD1
	AB=HD1  
	KK=K
	LK=L  
	AA=AA
	BB=BB  
	XX=XX  
	YY=YY  
	SMAG=(0.,0.) 
	CVAL=(0.,0.)   
	CALL CGQ1(AA,BB,XX,YY,10,CVAL)
	SMAG(:,1:K)=CVAL(:,1:K)  
	RETURN
	END
C
	SUBROUTINE FIELD(X,Y,Z,KKA,CF)
	IMPLICIT NONE
	INTEGER LK,KK
	COMMON/T13/ LK,KK
	REAL WLT,HHW,HHD,AB,SGN  
	REAL, ARRAY(LK,KK) :: RR,X1,X,Y,Z 
	COMPLEX, ARRAY(LK,KK) :: HMN,CF    
	COMMON/T3/ WLT,HHW,HHD,AB
	INTEGER, ARRAY(LK,KK) :: KKA     
	EXTERNAL SGN 
	RR=0.  
	X1=0.  
	HMN=(0.,0.)  
	CF=(0.,0.)   
	RR=2.*3.141592654/WLT*SQRT((X-Y)*(X-Y)+Z*Z)    
	X1=(KKA-1)*3.141592654*(X-HHD*SGN(AB))/HHW 
	CALL H20(RR,LK,KK,HMN)   
	CF=HMN*COS(X1)  
 	RETURN
	END
C
	SUBROUTINE AFHF(K,WL,RR,ANG,AA,BB,HW,HD1,FMAG)
	IMPLICIT NONE
	INTEGER L,II,NIJ,I,J,M,N
	COMMON/T1/ L,NIJ,I,J,M,N
	REAL WL,RR,HD1,WLT,HW,HWT,HDT
	INTEGER K,LK,KK
	COMMON/T5/ WLT,HWT,HDT   
	COMMON/T51/ LK,KK
	REAL, ARRAY(L) :: AA,BB,ANG 
	COMPLEX, ARRAY(L) :: ZW,ZEFF,ZK,ZEFFE
	COMPLEX, ARRAY(L,K) :: ZF,ZFF,CVAL,FMAG1
	COMPLEX, ARRAY(L,NIJ) :: FMAG
	COMPLEX ZJ 
	WLT=WL
	HWT=HW
	HDT=HD1
	KK=K
	LK=L  
	AA=AA  
	BB=BB  
	ANG=ANG  
	ZJ=(0.,1.)
	ZW=ZJ*WL/RR
	ZEFF=SQRT(ZW)/3.141592654
	ZF=SPREAD(ZEFF,DIM=2,NCOPIES=KK)
	ZK=-ZJ*2.*3.141592654/WL*RR
	ZEFFE=EXP(ZK)
	ZFF=SPREAD(ZEFFE,DIM=2,NCOPIES=KK)
	CALL CGQ2(AA,BB,RR,ANG,10,CVAL)
	FMAG1=CVAL*ZF*ZFF
	FMAG=(0.,0.)  
	FMAG(:,1:K)=FMAG1(:,1:K)   
	RETURN
	END
C
	SUBROUTINE APROX(X,Y,Z,KKA,APX)
	IMPLICIT NONE
	INTEGER LK,KK  
	COMMON/T5/ WLT,HWT,HDT  
	COMMON/T51/ LK,KK  
	COMPLEX, ARRAY(:,:) :: APX  
	COMPLEX, ARRAY(LK,KK) :: ZX,ZZX  
	COMPLEX ZJ
	REAL  WLT,HWT,HDT,Y    
	REAL, ARRAY(LK,KK) :: ZM,ZZM,Z,X
	INTEGER, ARRAY(:,:) ::  KKA
        ZJ=(0.,1.)  
	Y=Y  
        ZX=ZJ*2.*3.141592654*X/WLT*COS(Z)
        ZZX=EXP(ZX)
        ZM=(KKA-1)*3.141592654*(X-HDT)/HWT
        ZZM=COS(ZM)
        APX=ZZX*ZZM
	RETURN
	END
C
C
C----------SUBROUTINE--CGQ1-------------------------------------------
C                                                                     
      SUBROUTINE CGQ1(XXL,XXU,X0,Y0,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
	INTEGER LK,KK,N,L,NP,NE,J,NS(8),NQ(8)
	COMMON/T13/ LK,KK
	COMPLEX, ARRAY(LK,KK) :: CF,CF1,CF2,CF3,CVAL
	REAL, ARRAY(LK) :: X0,Y0,XXL,XXU
	REAL, ARRAY(LK,KK) :: XX0,YY0,AX,BX,DX,XL,XU
	INTEGER, ARRAY(KK) :: KA
	INTEGER, ARRAY(LK,KK) :: KKA  
        REAL Q1(52),Q2(24),Q3(32),QG(108)
        EQUIVALENCE (Q1(1),QG(1)),(Q2(1),QG(53)),(Q3(1),QG(77))
      DATA Q1
     $               /2.88675134594812882E-1,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/ 
	KA=[1:KK]
	KKA=0.  
	XL=0.  
	XU=0.  
	XX0=0.  
	YY0=0.  
	CF=(0.,0.)  
	CF1=(0.,0.)  
	CF2=(0.,0.)  
	CF3=(0.,0.)  
	CVAL=(0.,0.)  
	AX=0.  
	BX=0. 
	DX=0.   
	KKA=SPREAD(KA,DIM=1,NCOPIES=LK) 
	XL=SPREAD(XXL,DIM=2,NCOPIES=KK)
	XU=SPREAD(XXU,DIM=2,NCOPIES=KK) 
	XX0=SPREAD(X0,DIM=2,NCOPIES=KK) 
	YY0=SPREAD(Y0,DIM=2,NCOPIES=KK)
      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)  
      CALL FIELD(AX,XX0,YY0,KKA,CF)  
      CF1=CF  
      CVAL=CF1*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  
      CALL FIELD(AX+DX,XX0,YY0,KKA,CF) 
      CF2=CF  
      CALL FIELD(AX-DX,XX0,YY0,KKA,CF)  
      CF3=CF  
      CVAL=CVAL+CMPLX(QG(J+1),0.0)*(CF2+CF3)  
 350  CONTINUE 
      CVAL=CVAL*CMPLX(BX,0.0)
      RETURN 
      END 
C----------SUBROUTINE--CGQ2-------------------------------------------
c
C                                                                     
      	SUBROUTINE CGQ2(XXL,XXU,X0,Y0,N,CVAL)                                
        IMPLICIT NONE
	INTEGER LK,KK,N,L,NP,NE,J,NS(8),NQ(8)
	COMMON/T13/ LK,KK
	COMPLEX, ARRAY(LK,KK) :: APX,CF1,CF2,CF3,CVAL  
	REAL X0   
	REAL, ARRAY(LK) :: Y0,XXL,XXU
	REAL, ARRAY(LK,KK) :: YY0,AX,BX,DX,XL,XU
	INTEGER, ARRAY(KK) :: KA
	INTEGER, ARRAY(LK,KK) :: KKA  
        REAL Q1(52),Q2(24),Q3(32),QG(108)
        EQUIVALENCE (Q1(1),QG(1)),(Q2(1),QG(53)),(Q3(1),QG(77))
      DATA Q1
     $               /2.88675134594812882E-1,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/ 
	KA=[1:KK]
	KKA=SPREAD(KA,DIM=1,NCOPIES=LK) 
	XL=SPREAD(XXL,DIM=2,NCOPIES=KK)
	XU=SPREAD(XXU,DIM=2,NCOPIES=KK) 
	YY0=SPREAD(Y0,DIM=2,NCOPIES=KK)
      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)  
      CALL APROX(AX,X0,YY0,KKA,APX)  
      CF1=APX  
      CVAL=CF1*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  
      CALL APROX(AX+DX,X0,YY0,KKA,APX) 
      CF2=APX  
      CALL APROX(AX-DX,X0,YY0,KKA,APX)  
      CF3=APX  
      CVAL=CVAL+CMPLX(QG(J+1),0.0)*(CF2+CF3)  
 350  CONTINUE 
      CVAL=CVAL*CMPLX(BX,0.0)
      RETURN 
      END 
