C THIS PROGRAM IS WRITTEN FOR CALCULATION OF POLARIZATION CURRENT OF A
C DIELECTRIC POST IN A WAVEGUIDE SUPPORTING TE10.
C OCT, 1992.
C PROGRAM MODIFIED TO VARY THE DIELECTRIC CONSTANT
        REAL K0,LAM,H1
        COMPLEX E,I,G1,G2,GAMMA,TAU,ZP(200,200)
        COMPLEX P(200),Z(200,200),V(200),SV(200),SUM(200) 
        real aref(0:200),atrn(0:200),pref(0:200),ptrn(0:200)
        INTEGER KPVT(200) 
        COMMON /BLOCK1/K0,A
        COMMON /BLOCK2/DX,DY
        DATA MDIM /200/
C
        DO IE1=0,200
        I=(0.0,1.0)
        PI = 3.14159265358979324
C
        LAM=3.0E8/(4.0E9+1.0E7*IE1)  
        K0=(2.*PI)/LAM
        Y0=1./377.
        Z0=377.
        A=1.866*0.0254   
        E=(20.00,-7.00)        
        D=0.003556   
C COMPUTE POSITION VECTOR OF THE CELLS
        ITOT=112
        CALL POSITION(A,D,P)   
C
C COMPUTE THE IMPEDANCE MATRIX 

        DO 2 IND=1,ITOT
        DO 2 INDP=1,ITOT
         IF(IND .EQ. INDP) THEN  
          CALL SELF(P(IND),G1)
          Z(IND,INDP)=K0*Z0*4*G1/A
         ELSE
           y=aimag(p(ind))
           yp=aimag(p(indp))
          if(y .eq. yp)then
           call row(p(ind),p(indp),g2) 
          else 
           CALL IMP(P(IND),P(INDP),G2)
          endif 
           Z(IND,INDP)=4.*K0*Z0*G2/A
         ENDIF
          ZP(IND,INDP)=Z(IND,INDP)
2       CONTINUE
C       PRINT *,'FINISHED WITH THE IMPADANCE MATRIX CALCULATION' 
C COMPUTE INTEGRAL VECTOR  
        H1=SQRT(K0**2-(PI/A)**2)  
        DO  IND=1,ITOT 
        X=REAL(P(IND))
        Y=AIMAG(P(IND))
        SUM(IND)=SIN(PI*X/A)*CEXP(I*H1*Y)
        ENDDO   
C                      
        MIND=ITOT
C        E=CMPLX(2.0,0.0)
         DO IND=1,ITOT 
         DO INDP=1,ITOT
         IF(IND .EQ. INDP)THEN
         Z(IND,IND)=ZP(IND,IND)+I/(K0*Y0*(E-1.))
         ELSE
         Z(IND,INDP)=ZP(IND,INDP)
         ENDIF
         ENDDO
         ENDDO
C----------------------------------------------------------------
CCC   L-U DECOMPOSE MATRIX [Z] AND SOLVE FOR CURRENTS
C          
        CALL CGECO(Z,MDIM,MIND,KPVT,RC,SV)
C       PRINT *,'FINISHED WITH THE MATRIX INVERSION'
C       PRINT *,'CONDITION IS:',RC             
C COMPUTE EXCITATION VECTOR
        DO  IND=1,ITOT 
        X=REAL(P(IND))
        Y=AIMAG(P(IND))
        V(IND)=SIN(PI*X/A)*CEXP(I*H1*Y)
        ENDDO 
C
C COMPUTE THE INDUCED CURRENT
        CALL CGESL(Z,MDIM,MIND,KPVT,V,0)
C COMPUTE REFLECTION AND TRANSMISSION COEFFICIENT 
        GAMMA=(0.0,0.0)               
        TAU=(0.0,0.0)
        DO IND=1,ITOT
        GAMMA=GAMMA+V(IND)*SUM(IND)
        TAU=TAU+V(IND)*CONJG(SUM(IND))
        ENDDO
        GAMMA=-GAMMA*K0*Z0*DX*DY/(A*H1)
        TAU=1.0-TAU*K0*Z0*DX*DY/(A*H1)
C       PRINT*,CABS(GAMMA)**2,CABS(TAU)**2  
        aref(ie1)=20.*alog10(cabs(GAMMA))
        atrn(ie1)=20.*alog10(cabs(tau))
        pref(ie1)=phase(gamma) 
        ptrn(ie1)=phase(tau)
        ENDDO            
        open(10,file='bw10rm.dat')
        open(11,file='bw10rp.dat')
        open(12,file='bw10tm.dat')
        open(13,file='bw10tp.dat')
        WRITE(10,101)(AREF(IND1),IND1=0,200)
        WRITE(11,101)(PREF(IND1),IND1=0,200)
        WRITE(12,101)(ATRN(IND1),IND1=0,200)
        WRITE(13,101)(PTRN(IND1),IND1=0,200)
101   FORMAT(1H ,F7.2)   
        STOP
        END
C*****************************************************************C
C                                                                 C
        SUBROUTINE POSITION(A,D,P)                     
C                                                                 C
C*****************************************************************C
C THIS SUBROUTINE IS CALLED FROM MAIN PROGRAM TO CALCULATE THE
C POSITION VECTOR OF THE CELLS
C
        COMPLEX P(200)
        REAL A  
        INTEGER IR(6)
        COMMON /BLOCK2/DX,DY 
        DY=D/FLOAT(12)
        DX=DY
        IR(1)=4
        IR(2)=8
        IR(3)=10
        IR(4)=10
        IR(5)=12
        IR(6)=12
        K=0
        DO M=1,6
        DO N=1,IR(M)
        K=K+1
        P(K)=CMPLX((N-IR(M)/2.-1/2.)*DX, (7-M-1/2.)*DY )
        P(113-K)=-P(K)
        ENDDO
        ENDDO   
        do ind=1,112
        p(ind)=a/2.+p(ind)
c        write(20,*)100*real(p(ind)),100*aimag(p(ind)) 
        enddo
        RETURN
        END       
C           
C*****************************************************************C
C                                                                 C
        SUBROUTINE SELF(P,G1)                     
C                                                                 C
C*****************************************************************C
C THIS SUBROUNTINE IS CALLED FROM MAIN TO EVALUATED SELF CELL FOR
C IMPEDANCE MATRIX
        COMPLEX G1,I,HN,ONE,P
        REAL K0,DX,DY,KNX,A 
        COMMON /BLOCK1/K0,A
        COMMON /BLOCK2/DX,DY
        I=(0.0,1.0) 
        ONE=(1.0,0.0)
        PI = 3.14159265358979324
        G1=(0.0,0.0)
        X=REAL(P)
        DO N=1,50 
        KNX=PI*N/A
        HN=CSQRT(K0**2-ONE*KNX**2)
        G1=-I*(SIN(KNX*X))**2*(CEXP(I*HN*DY/2.)-1.)*SIN(KNX*DX/2.)/
     &     (KNX*HN**2)+G1   
        ENDDO  
        RETURN   
        END
C*****************************************************************C
C                                                                 C
        SUBROUTINE ROW(P,PP,G)                     
C                                                                 C
C*****************************************************************C
C THIS SUBROUNTINE IS CALLED FROM MAIN TO EVALUATED SELF CELL FOR
C IMPEDANCE MATRIX
        COMPLEX G,I,HN,ONE,P,PP
        REAL K0,DX,DY,KNX,A 
        COMMON /BLOCK1/K0,A
        COMMON /BLOCK2/DX,DY
        I=(0.0,1.0) 
        ONE=(1.0,0.0)
        PI = 3.14159265358979324
        G=(0.0,0.0)
        X=REAL(P)
        XP=REAL(PP)
        DO N=1,50 
        KNX=PI*N/A
        HN=CSQRT(K0**2-ONE*KNX**2)
        G=-I*SIN(KNX*X)*SIN(KNX*XP)*SIN(KNX*DX/2.)*
     &     (CEXP(I*HN*DY/2.)-1.)/(KNX*HN**2)+G   
        ENDDO  
        RETURN   
        END
C*****************************************************************C
C                                                                 C
        SUBROUTINE IMP(P,PP,G)                     
C                                                                 C
C*****************************************************************C
C THIS SUBROUNTINE IS CALLED FROM MAIN TO EVALUATED ENTRIES FOR
C IMPEDANCE MATRIX WHEN ROW OF OBSERVATION AND SOURCE ARE EQUAL.
        COMPLEX G,I,HN,ONE,P,PP
        REAL K0,DX,DY,KNX 
        COMMON /BLOCK1/K0,A
        COMMON /BLOCK2/DX,DY
        I=(0.0,1.0) 
        ONE=(1.0,0.0)
        PI = 3.14159265358979324
        G=(0.0,0.0)
        X=REAL(P)
        XP=REAL(PP)
        Y=AIMAG(P)
        YP=AIMAG(PP)
        DO N=1,25 
        KNX=(PI*N/A)
        HN=CSQRT(K0**2-ONE*KNX**2) 
        G=SIN(KNX*X)*SIN(KNX*XP)*SIN(KNX*DX/2.)*SIN(HN*DY/2.)*
     &     CEXP(I*HN*ABS(Y-YP))/(HN**2*KNX)+G  
c        PRINT *,sin(hn*dy/2),G
        ENDDO
        RETURN
        END
C*******************************************************************
C
        FUNCTION PHASE(Z)
C
C******************************************************************
        COMPLEX Z 
        PI=4.*ATAN(1.)
        X=REAL(Z)
        Y=AIMAG(Z)
        if(x.eq.0 .and. y.eq.0)then
        phase=0
        else
        PHASE=(180./PI)*ATAN2(Y,X)
        endif
        RETURN
        END