
C THIS PROGRAM IS 
C IS WRITTEN FOR CALCULATION OF DIELECTRIC CONSTANT OF A CIRCULAR
C DIELECTRIC POST IN A WAVEGUIDE SUPPORTING TE10.
C THE METHOD IS BASED ON EIGEN ANALYSIS THAT IS IT FINDS
C THE EIGEN VALUES AND EIGEN VECTORS OF THE INCOMPLETE IMPEDANCE MATRIX
C AND THEN INVERTED IMPEDANCE MATRIX IS OBTAINED IN A EXPLICIT FORM IN TERMS
C OF THE DIELECTRIC CONSTANT.                   
C THE PROGRAM MUST BE BOUND WITH THE MTRXINV.BIN AND EISPAK.BIN
C In Sun, use
C   f77 -cg89 -o rditw31 dit.f mtrxinv.f -L/usr/um/lib -leispack_d

        REAL K0,LAM,H1,E(2),EER(0:200),EEI(0:200)
        COMPLEX I,G1,G2,GAMMA,GAMMA0,GGAA(0:200)
        COMPLEX P(200),Z(200,200),V(200),SV1(200),SUM(200) 
        COMPLEX Q(200,200),QINV(200,200),W(200)
        INTEGER KPVT(200) 
       DOUBLE PRECISION AR(200,200),AI(200,200)
       DOUBLE PRECISION QR(200,200),QI(200,200),WR(200),WI(200),
     &                 SV2(200),SV3(200),SV4(200) 
C        COMPLEX DG(200,200),DUM1(200,200),DUM2(200,200)
        COMPLEX LL(200),RR(200),TLL(200)
	character mfile*11,pfile*11,rfile*11,ifile*11
        COMMON /BLOCK1/K0,A
        COMMON /BLOCK2/DX,DY  
        COMMON /BLOCK3/ RR,LL,W,GAMMA0,H1,ITOT
        DATA MDIM /200/
C
        I=(0.0,1.0)
        PI = 3.14159265358979324
C                  

	print*, 'Enter mag data file name'
	read(*,'(a11)') mfile

	print*, 'Enter phase data file name'
	read(*,'(a11)') pfile

        open(66,file=mfile)
        open(67,file=pfile)
        DO IIIP=0,200
        READ(66,*)RHO
        READ(67,*)PHI                        
        PHI=-1.0*PHI*PI/180.0
        RHO=RHO/20.0
        RHO=10.**RHO
        GGAA(IIIP)=CMPLX(RHO*COS(PHI),RHO*SIN(PHI))
        ENDDO
        E(1)=40.0
        E(2)=20.0
C
        DO IE1=0,200
        GAMMA0=GGAA(IE1) 
        
        LAM=3.0E8/(4.0E9+2.0E7*IE1)  
        K0=(2.*PI)/LAM
        Z0=377.
        A=1.866*0.0254  
C       {WAVEGUIDE WIDTH in meters  

	print*, 'Enter diameter of post (meters)'
	read(*,*) d
C       {DIAMETER OF CYLINDRICAL POST 
C       PRINT*,'INPUT THE MEASURED REF. COEFFICIENT MAG. AND PHASE'
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
2       CONTINUE
C        PRINT *,'FINISHED WITH THE IMPADANCE MATRIX CALCULATION' 
        DO IND=1,ITOT
        DO INDP=1,ITOT
        work1= REAL(Z(IND,INDP))
        AR(IND,INDP)=DBLE(work1)
        work1=AIMAG(Z(IND,INDP))
        AI(IND,INDP)=DBLE(work1)
        ENDDO
        ENDDO
C-----------------------------------------------------------------
        MATZ=1    
C       {RETURNS BOTH THE EIGENVALUES AND THE EIGEN VECTORS
        CALL CG(MDIM,ITOT,AR,AI,WR,WI,MATZ,QR,QI,SV2,SV3,SV4,IER)
C        PRINT *,'FINISHED WITH EIGEN-ANALYSIS'
C  EIGEN VALUES
        DO IND=1,ITOT
        W(IND)=CMPLX(SNGL(WR(IND)),SNGL(WI(IND)))
        ENDDO
C EIGEN MATRIX
        DO IND=1,ITOT
        DO INDP=1,ITOT
        Q(IND,INDP)=CMPLX(SNGL(QR(IND,INDP)),SNGL(QI(IND,INDP)))
        QINV(IND,INDP)=Q(IND,INDP)
        ENDDO
        ENDDO
        MIND=ITOT
C----------------------------------------------------------------
CCC   L-U DECOMPOSE MATRIX [Q] 
C          
        CALL CGECO(QINV,MDIM,MIND,KPVT,RC,SV1)
C        PRINT *,'FINISHED WITH THE LU DECOMPOSITION'
C        PRINT *,'CONDITION IS:',RC                 
C----------------------------------------------------------------
        CALL CGEDI(QINV,MDIM,MIND,KPVT,DET,SV1,1)
CCC   INVESION OF MATRIX [Q]
C        PRINT*,'FINISHED WITH INVERSION OF [Q]'
C
C COMPUTE EXCITATION VECTOR
        H1=SQRT(K0**2-(PI/A)**2)  
        DO  IND=1,ITOT 
        X=REAL(P(IND))
        Y=AIMAG(P(IND))
        V(IND)=SIN(PI*X/A)*CEXP(I*H1*Y)
        ENDDO   
C COMPUTE RR=[Q]^-1V
        DO IND=1,ITOT
        RR(IND)=(0.0,0.0)
        DO INDP=1,ITOT
        RR(IND)=RR(IND)+QINV(IND,INDP)*V(INDP)
        ENDDO
        ENDDO 
C COMPUTE LL=V^t[Q]
        DO IND=1,ITOT
        LL(IND)=(0.0,0.0)
        DO INDP=1,ITOT
        LL(IND)=LL(IND)+V(INDP)*Q(INDP,IND)
        ENDDO
        ENDDO
C BEGIN INVERSION PROCESS
        n=2
        ftol=1.e-6
        call frprmn(E,n,ftol,iter,fret) 
C        print*,iter,fret 
        EER(IE1)=E(1)
        EEI(IE1)=E(2)
C        PRINT*,E(1),E(2)
        ENDDO   
	print*, 'Enter real output data file name'
	read(*,'(a11)') rfile

	print*, 'Enter imaginary output data file name'
	read(*,'(a11)') ifile


        open(101,file=rfile)
        open(111,file=ifile)
        WRITE(101,101)(EER(IND1),IND1=0,200)
        WRITE(111,101)(EEI(IND1),IND1=0,200)   
101     FORMAT(1H ,F12.7)
        STOP
        END 
C*****************************************************************
C
        FUNCTION FUNC(X)                                                
C
C******************************************************************   
C THIS FUNCTION CALCULATES REFLECTION COEFFICIENT IN TERMS OF 
C REAL AND IMAGINARY PART OF DIELECTIC CONSTANT X(1)=Er  X(2)=Ei
        COMPLEX RR(200),LL(200),W(200),I,GAMMA,GAMMA0,E
        REAL K0,Y0,X(2),H1
        COMMON /BLOCK1/K0,A
        COMMON /BLOCK2/DX,DY  
        COMMON /BLOCK3/ RR,LL,W,GAMMA0,H1,ITOT
        Z0=377.
        Y0=1./377.
        I=(0.0,1.0)
        E=CMPLX(X(1),X(2))        
        GAMMA=(0.0,0.0)
        DO IND=1,ITOT
        GAMMA=GAMMA+RR(IND)*LL(IND)/(W(IND)+I/(K0*Y0*(E-1.)))
        ENDDO
        GAMMA=-GAMMA*K0*Z0*DX*DY/(A*H1)
        FUNC=1000*CABS(GAMMA-GAMMA0)**2
        RETURN
        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 K,M,N,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
c        open(20,file='pos')
        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
C**********************************************************************
C
        SUBROUTINE MTMULT(N,A,B,C)
C
C***********************************************************************
        COMPLEX A(200,200),B(200,200),C(200,200)
        DO IND=1,N
        DO INDP=1,N
        C(IND,INDP)=(0.0,0.0) 
          DO K=1,N
          C(IND,INDP)=C(IND,INDP)+A(IND,K)*B(K,INDP)
          ENDDO
        ENDDO
        ENDDO
        RETURN
        END
c****************************************************************************
c       
        subroutine dfunc(x,gf)
c
c****************************************************************************
c        
c this subroutine computes the gadient func
c
        real x(2),gf(2),xp(2)
        delta=1.e-2 
        f=func(x)
        do j=1,2
        do ind1=1,2
        xp(ind1)=x(ind1)
        enddo
        xp(j)=x(j)+delta
        f2=func(xp)
        gf(j)=(f2-f)/(delta)
        enddo     
        return
        end
c*******************************************************************************
c
      SUBROUTINE FRPRMN(P,N,FTOL,ITER,FRET) 
c
c*****************************************************************************
c
c  Conjugate Gradient Method taken from Numerical Recepies pp. 305
c
      PARAMETER (NMAX=50,ITMAX=200,EPS=1.E-9)
      DIMENSION P(N),G(NMAX),H(NMAX),XI(NMAX)
      FP=FUNC(P) 
      CALL DFUNC(P,XI)
      DO 11 J=1,N
        G(J)=-XI(J)
        H(J)=G(J)
        XI(J)=H(J)
11    CONTINUE
      DO 14 ITER=1,ITMAX
        CALL LINMIN(P,XI,N,FRET) 
        IF(2.*ABS(FRET-FP).LE.FTOL*(ABS(FRET)+ABS(FP)+EPS))RETURN
        FP=FUNC(P)
        CALL DFUNC(P,XI)
        GG=0.
        DGG=0.
        DO 12 J=1,N
          GG=GG+G(J)**2
C         DGG=DGG+XI(J)**2
          DGG=DGG+(XI(J)+G(J))*XI(J)
12      CONTINUE   
        IF(GG.EQ.0.)RETURN
        GAM=DGG/GG
        DO 13 J=1,N
          G(J)=-XI(J)
          H(J)=G(J)+GAM*H(J)
          XI(J)=H(J)
13      CONTINUE
14    CONTINUE
      PAUSE 'FRPR maximum iterations exceeded'
      RETURN
      END
c********************************************************************
      FUNCTION DF1DIM(X)
      PARAMETER (NMAX=50)
      COMMON /F1COM/ NCOM,PCOM(NMAX),XICOM(NMAX)
      DIMENSION XT(NMAX),DF(NMAX)
      DO 11 J=1,NCOM
        XT(J)=PCOM(J)+X*XICOM(J)
11    CONTINUE
      CALL DFUNC(XT,DF)
      DF1DIM=0.
      DO 12 J=1,NCOM
        DF1DIM=DF1DIM+DF(J)*XICOM(J)
12    CONTINUE
      RETURN
      END
c*************************************************************************        
      FUNCTION F1DIM(X)
      PARAMETER (NMAX=50)
      COMMON /F1COM/ NCOM,PCOM(NMAX),XICOM(NMAX)
      DIMENSION XT(NMAX)
      DO 11 J=1,NCOM
        XT(J)=PCOM(J)+X*XICOM(J)
11    CONTINUE
      F1DIM=FUNC(XT)
      RETURN
      END
c****************************************************************************

      SUBROUTINE LINMIN(P,XI,N,FRET)
      PARAMETER (NMAX=50,TOL=1.E-4)
      EXTERNAL F1DIM
      EXTERNAL DF1DIM
      DIMENSION P(N),XI(N)
      COMMON /F1COM/ NCOM,PCOM(NMAX),XICOM(NMAX)
      NCOM=N
      DO 11 J=1,N
        PCOM(J)=P(J)
        XICOM(J)=XI(J)
11    CONTINUE
      AX=0.
      XX=1.
      BX=2.
      CALL MNBRAK(AX,XX,BX,FA,FX,FB,F1DIM)
      FRET=DBRENT(AX,XX,BX,F1DIM,DF1DIM,TOL,XMIN)
      DO 12 J=1,N
        XI(J)=XMIN*XI(J)
        P(J)=P(J)+XI(J)
12    CONTINUE
      RETURN
      END
c************************************************************************** 
      FUNCTION DBRENT(AX,BX,CX,F,DF,TOL,XMIN)
      PARAMETER (ITMAX=100,ZEPS=1.0E-10)
      LOGICAL OK1,OK2
      A=MIN(AX,CX)
      B=MAX(AX,CX)
      V=BX
      W=V
      X=V
      E=0.
      FX=F(X)
      FV=FX
      FW=FX
      DX=DF(X)
      DV=DX
      DW=DX
      DO 11 ITER=1,ITMAX
        XM=0.5*(A+B)
        TOL1=TOL*ABS(X)+ZEPS
        TOL2=2.*TOL1
        IF(ABS(X-XM).LE.(TOL2-.5*(B-A))) GOTO 3
        IF(ABS(E).GT.TOL1) THEN
          D1=2.*(B-A)
          D2=D1
          IF(DW.NE.DX) D1=(W-X)*DX/(DX-DW)
          IF(DV.NE.DX) D2=(V-X)*DX/(DX-DV)
          U1=X+D1
          U2=X+D2
          OK1=((A-U1)*(U1-B).GT.0.).AND.(DX*D1.LE.0.)
          OK2=((A-U2)*(U2-B).GT.0.).AND.(DX*D2.LE.0.)
          OLDE=E
          E=D
          IF(.NOT.(OK1.OR.OK2))THEN
            GO TO 1
          ELSE IF (OK1.AND.OK2)THEN
            IF(ABS(D1).LT.ABS(D2))THEN
              D=D1
            ELSE
              D=D2
            ENDIF
          ELSE IF (OK1)THEN
            D=D1
          ELSE
            D=D2
          ENDIF
          IF(ABS(D).GT.ABS(0.5*OLDE))GO TO 1
          U=X+D
          IF(U-A.LT.TOL2 .OR. B-U.LT.TOL2) D=SIGN(TOL1,XM-X)
          GOTO 2
        ENDIF
1       IF(DX.GE.0.) THEN
          E=A-X
        ELSE
          E=B-X
        ENDIF
        D=0.5*E
2       IF(ABS(D).GE.TOL1) THEN
          U=X+D
          FU=F(U)
        ELSE
          U=X+SIGN(TOL1,D)
          FU=F(U)
          IF(FU.GT.FX)GO TO 3
        ENDIF
        DU=DF(U)
        IF(FU.LE.FX) THEN
          IF(U.GE.X) THEN
            A=X
          ELSE
            B=X
          ENDIF
          V=W
          FV=FW
          DV=DW
          W=X
          FW=FX
          DW=DX
          X=U
          FX=FU
          DX=DU
        ELSE
          IF(U.LT.X) THEN
            A=U
          ELSE
            B=U
          ENDIF
          IF(FU.LE.FW .OR. W.EQ.X) THEN
            V=W
            FV=FW
            DV=DW
            W=U
            FW=FU
            DW=DU
          ELSE IF(FU.LE.FV .OR. V.EQ.X .OR. V.EQ.W) THEN
            V=U
            FV=FU
            DV=DU
          ENDIF
        ENDIF
11    CONTINUE
      PAUSE 'DBRENT exceeded maximum iterations.'
3     XMIN=X
      DBRENT=FX
      RETURN
      END 
c**********************************************************************       

      SUBROUTINE MNBRAK(AX,BX,CX,FA,FB,FC,FUNCC)
      PARAMETER (GOLD=1.618034, GLIMIT=100., TINY=1.E-20)
      FA=FUNCC(AX)
      FB=FUNCC(BX)
      IF(FB.GT.FA)THEN
        DUM=AX
        AX=BX
        BX=DUM
        DUM=FB
        FB=FA
        FA=DUM
      ENDIF
      CX=BX+GOLD*(BX-AX)
      FC=FUNCC(CX)
1     IF(FB.GE.FC)THEN
        R=(BX-AX)*(FB-FC)
        Q=(BX-CX)*(FB-FA)
        U=BX-((BX-CX)*Q-(BX-AX)*R)/(2.*SIGN(MAX(ABS(Q-R),TINY),Q-R))
        ULIM=BX+GLIMIT*(CX-BX)
        IF((BX-U)*(U-CX).GT.0.)THEN
          FU=FUNCC(U)
          IF(FU.LT.FC)THEN
            AX=BX
            FA=FB
            BX=U
            FB=FU
            GO TO 1
          ELSE IF(FU.GT.FB)THEN
            CX=U
            FC=FU
            GO TO 1
          ENDIF
          U=CX+GOLD*(CX-BX)
          FU=FUNCC(U)
        ELSE IF((CX-U)*(U-ULIM).GT.0.)THEN
          FU=FUNCC(U)
          IF(FU.LT.FC)THEN
            BX=CX
            CX=U
            U=CX+GOLD*(CX-BX)
            FB=FC
            FC=FU
            FU=FUNCC(U)
          ENDIF
        ELSE IF((U-ULIM)*(ULIM-CX).GE.0.)THEN
          U=ULIM
          FU=FUNCC(U)
        ELSE
          U=CX+GOLD*(CX-BX)
          FU=FUNCC(U)
        ENDIF
        AX=BX
        BX=CX
        CX=U
        FA=FB
        FB=FC
        FC=FU
        GO TO 1
      ENDIF
      RETURN
      END






