      PROGRAM DTE
C
C     THIS PROGRAM COMPUTES NEAR FIELDS FOR THE
C     TE POLARIZATION FROM PEC AND LOSSY CYLINDERS.
C     THE PROGRAM PROVIDES FIELD STRENGTH IN dB V/M.
C     THE INCIDENCE ANGLE IS ASSUMED TO BE 180.0 DEGREES.
C
C     WRITTEN BY:  KEITH LYSIAK
C     PORTIONS TAKEN FROM JOHN BEGGS' TEBIST PROGRAM
C     DATE:  15 APRIL 1996
C
C     PROGRAM INPUTS:  dte.inp  (DATA FILE OF INPUT PARAMETERS)
C     DATA FILE FORMAT:  FREQ
C                        XWIDTH,XINCR,YWIDTH,YINCR
C                        RADIUS,MODE
C                        EPSR,SIGMA
C                        XMUR,SIGMAC
C     
C     FREQ =   FREQUENCY OF INTEREST
C     XWIDTH = CALCULATION WIDTH (IN METERS)  \
C     XINCR =  WIDTH INCREMENT (IN METERS)     | CENTERED AROUND
C     XWIDTH = CALCULATION WIDTH (IN METERS)   | CYLINDER
C     XINCR =  WIDTH INCREMENT (IN METERS)    /
C     RADIUS = RADIUS OF THE CYLINDER IN METERS
C     MODE =   DETERMINES TYPE OF CYLINDER (MODE=0 FOR PEC CYLINDER,
C              MODE=1 FOR DIELECTRIC OR MAGNETIC)
C     EPSR =   RELATIVE DIELECTRIC CONSTANT
C     SIGMA =  CONDUCTIVITY (IN SIEMENS PER METER)
C     XMUR =   RELATIVE PERMEABILITY
C     SIGMAC = MAGNETIC CONDUCTIVITY
C
C     PROGRAM OUTPUTS:  dte.dat (DATA FILE OF ELECTRIC FIELDS
C                       VERSUS X/Y LOCATIONS)
C
C     DATA FILE FORMAT:  X,Y,Ex,Ey,Erho,Ephi
C
C     PROGRAM CALLS:  FUNCTIONS BESSJ, BESSY, HANK, BESSJP, HANKP
C                     BESSJ0, BESSY0, BESSJ1, BESSY1
C
C     FOR EVALUATING BESSEL FUNCTIONS OF COMPLEX ARGUMENTS
C     AND THEIR DERIVATIVES
C
      IMPLICIT NONE

      EXTERNAL RHO,PHI

      REAL*8 RHO,PHI,XWIDTH,YWIDTH,XINCR,YINCR,X,Y,RADIUS
      REAL*8 W,FREQ,LAMBDA,EPS0,XMU0,PI
      REAL*8 EPSR,SIGMA,XMUR,SIGMAC
      REAL*8 EXMAG,EYMAG,Rmag,Pmag
      
      COMPLEX*16 J,EX,EY
      COMPLEX*16 K0,K1,EPSRC,XMURC,P,R
      COMPLEX*16 AN0,BN0,AN(115),BN(115)
      
      INTEGER MODE

      PARAMETER (EPS0=8.841941E-12,XMU0=1.256637E-6)
      PARAMETER (J=(0.0,1.0))
      PARAMETER (PI=3.1415926536)

      COMMON /PAR/ RADIUS,K0,K1,W,EPSRC,XMURC,MODE
  
      OPEN (UNIT=20, FILE='dte.dat',STATUS='unknown')
      
      CALL INPUT(FREQ,XWIDTH,XINCR,YWIDTH,YINCR,RADIUS,MODE,EPSR,
     $ SIGMA,XMUR,SIGMAC)
      CALL SETUP(W,FREQ,LAMBDA,K0,K1,EPSR,SIGMA,EPSRC,XMUR,SIGMAC,
     $ XMURC,MODE,RADIUS,AN0,BN0,AN,BN)
      
C     BEGIN AREA CALCULATIONS OVER RANGE OF X AND Y

      DO 10 X=-XWIDTH/2.0,XWIDTH/2.0,XINCR
	DO 20 Y=-YWIDTH/2.0,YWIDTH/2.0,YINCR

C     CALCULATE EX AND EY

	  CALL EPHI(X,Y,P,AN0,BN0,AN,BN)
	  CALL ERHO(X,Y,R,AN,BN)
      
	  EX = (R)*DCOS(PHI(X,Y))-(P)*DSIN(PHI(X,Y))
	  EY = (R)*DSIN(PHI(X,Y))+(P)*DCOS(PHI(X,Y))
      
	  EXMAG=CDABS(EX)
	  EYMAG=CDABS(EY)
	  
	  Rmag=cdabs(R)
	  Pmag=cdabs(P)

C     WRITE X,Y,EX, AND EY TO exact.dat

	  WRITE (20,200) X,Y,EXMAG,EYMAG,Rmag,Pmag
 20     CONTINUE
 10   CONTINUE

 200  FORMAT (2(F6.2,1X),4(e10.4,1X))

      CLOSE (UNIT=20)
      STOP
      END

*--------------------------------------------------------------------*
      SUBROUTINE INPUT(FREQ,XW,XINCR,YW,YINCR,R,MODE,EPSR,SIGMA,XMUR,
     $ SIGMAC)
*--------------------------------------------------------------------*
      IMPLICIT NONE
      
      REAL*8 XW,XINCR,XSAMP,YW,YINCR,YSAMP,R
      REAL*8 FREQ,EPSR,SIGMA,XMUR,SIGMAC
      
      INTEGER MODE

      OPEN (UNIT=10,FILE='dte.inp',STATUS='OLD')
      
      READ (10,*) FREQ
      READ (10,*) XW,XSAMP,YW,YSAMP
      READ (10,*) R,MODE
      READ (10,*) EPSR,SIGMA
      READ (10,*) XMUR,SIGMAC
      
      XINCR = XW/XSAMP
      YINCR = YW/YSAMP

      CLOSE (UNIT=10)

      RETURN
      END
*--------------------------------------------------------------------*
      SUBROUTINE SETUP(W,FREQ,LAMBDA,K0,K1,EPSR,SIGMA,EPSRC,XMUR,
     & SIGMAC,XMURC,MODE,A,AN0,BN0,AN,BN)
*--------------------------------------------------------------------*
      IMPLICIT NONE
      EXTERNAL BESSJ,BESSJP,HANK,HANKP
      
      INTEGER MODE,N
      REAL*8 PI,W,FREQ,LAMBDA,EPS0,EPSR,SIGMA,XMU0,XMUR,SIGMAC
      REAL*8 A
      COMPLEX*16 J,EPSRC,XMURC,K0,K1
      COMPLEX*16 AN0,BN0,BESSJ,BESSJP,HANK,HANKP,AN(115),BN(115)
      PARAMETER (PI=3.1415926536)
      PARAMETER (EPS0=8.841941E-12,XMU0=1.256637E-6)
      PARAMETER (J=(0.0,1.0))
      

C  INITIALIZE VARIABLES

      SIGMAC=SIGMAC*XMU0/EPS0
      W=2.0*PI*FREQ
      LAMBDA=3.0E8/FREQ
      K0=DCMPLX(W/3.0E8,0.0)
      EPSRC=EPSR-J*SIGMA/(W*EPS0)
      XMURC=XMUR-J*SIGMAC/(W*XMU0)
      K1=W*CDSQRT(EPSRC*XMURC)/3.0E8
      
C  INITIALIZE COEFFICIENTS

C  BEGIN WITH ZERO TERMS FOR EPHI
      
      IF (MODE.EQ.0) THEN
        AN0=BESSJP(0,K0*A)/HANKP(2,0,K0*A)
      ELSE IF (MODE.EQ.1) THEN
        AN0=(CDSQRT(XMURC/EPSRC)*BESSJ(0,K0*A)*BESSJP(0,K1*A)-
     $  BESSJP(0,K0*A)*BESSJ(0,K1*A))/(CDSQRT(XMURC/EPSRC)*
     $  BESSJP(0,K1*A)*HANK(2,0,K0*A)-BESSJ(0,K1*A)*HANKP(2,0,K0*A))
      ENDIF
      
      DO 10 N=1,115
      AN(N) = (0.0,0.0)
 10   CONTINUE

      BN0=(BESSJP(0,K0*A)*HANK(2,0,K0*A)-
     $  BESSJ(0,K0*A)*HANKP(2,0,K0*A))
     $  /(CDSQRT(XMURC/EPSRC)*BESSJP(0,K1*A)*HANK(2,0,K0*A)-
     $  BESSJ(0,K1*A)*HANKP(2,0,K0*A))
      
      DO 20 N=1,115
      BN(N) = (0.0,0.0)
 20   CONTINUE
     
      RETURN
      END
*--------------------------------------------------------------------*
      SUBROUTINE EPHI(X,Y,P,AN0,BN0,AN,BN)
*--------------------------------------------------------------------*
      IMPLICIT NONE
      
      EXTERNAL BESSJ,BESSJP,HANK,HANKP,RHO,PHI
      
      REAL*8 EPS0,H0,W,X,Y,RADIUS,RHO,PHI,A
      
      COMPLEX*16 TERM,BESSJ,BESSJP,HANK,HANKP,K0,K1
      COMPLEX*16 AN0,AN(115),BN0,BN(115)
      COMPLEX*16 EPSRC,XMURC,P,J
      
      INTEGER N,NTERM,MODE
      
      PARAMETER (EPS0=8.841941E-12)
      PARAMETER (J=(0.0,1.0))
      
      COMMON /PAR/ RADIUS,K0,K1,W,EPSRC,XMURC,MODE

      A=RADIUS
      P=(1.0e-4,1.0e-4)
      H0=1.0/377.0
      
C  DO FIELDS INSIDE REAL CYLINDER   
      
      IF (RHO(X,Y).LT.RADIUS) THEN
      
C  IF PEC (MODE=0) SKIP CALCULATING FIELDS INSIDE   
      
	IF (MODE.EQ.0) GO TO 50

	IF (RHO(X,Y).LT.1.0E-3) GO TO 50
	
	  P=BN0*BESSJP(0,K1*RHO(X,Y))

C  BEGIN SERIES SUMMATION   
	  
	NTERM=0
	DO 10 N=1,115
	  IF (CDABS(BN(N)).GT.0.0) GOTO 888
c	    write(*,*) 'computing bn for term ',n
	    BN(N) = ((BESSJP(N,K0*A)*HANK(2,N,K0*A)-
     $        BESSJ(N,K0*A)*HANKP(2,N,K0*A))
     $        /(CDSQRT(XMURC/EPSRC)*BESSJP(N,K1*A)*HANK(2,N,K0*A)-
     $        BESSJ(N,K1*A)*HANKP(2,N,K0*A)))
c      write(*,*) 'bn= ',bn(n)
888       CONTINUE
		
       TERM=((-J)**N)*2.0*BN(N)*BESSJP(N,K1*RHO(X,Y))*DCOS(N*PHI(X,Y))

C   DETERMINE IF SUM IS CONVERGING  ****

	  IF (CDABS(TERM/(P+TERM)).GT.(1E-3)) THEN

C  NOT YET, ADD IN NEW VALUE  ****

	    P=P+TERM
	    NTERM=0
	  ELSE
	    NTERM=NTERM+1

C    IF CONVERGED FOR 3 CONSECUTIVE ORDERS, THEN YOU'RE FINISHED

	    IF(NTERM.EQ.3) GO TO 20
	    P=P+TERM
	  ENDIF

 10     CONTINUE
 
	WRITE (6,*) 'DID NOT CONVERGE, RHO= ',RHO(X,Y),' PHI= ',PHI(X,Y)

C  MULTIPLY FINAL TERM BY CONSTANTS   
      
 20     CONTINUE
        P = ((-K0*H0)/(J*W*EPS0))*CDSQRT(XMURC/EPSRC)*P

        ELSE	  

C  DO FIELDS OUTSIDE CYLINDER   
	
	P=BESSJP(0,K0*RHO(X,Y))-(AN0*HANKP(2,0,K0*RHO(X,Y)))
	
	NTERM=0
	DO 30 N=1,115
	  IF (CDABS(AN(N)).GT.0.0) GOTO 999
c	    write(*,*) 'computing AN for term ',n
      
            IF (MODE.EQ.0) THEN
	    AN(n)=BESSJP(N,K0*A)/HANKP(2,N,K0*A)
	    ELSE IF (MODE.EQ.1) THEN
	    AN(n)=(CDSQRT(XMURC/EPSRC)*BESSJ(N,K0*A)*BESSJP(N,K1*A)-
     $        BESSJP(N,K0*A)*BESSJ(N,K1*A))/(CDSQRT(XMURC/EPSRC)*
     $        BESSJP(N,K1*A)*HANK(2,N,K0*A)-BESSJ(N,K1*A)*
     $        HANKP(2,N,K0*A))
	    ENDIF
c	    write(*,*) 'AN= ',AN(n)

999       CONTINUE	  
	  
	  TERM=((-J)**N)*2.0*
     $      (BESSJP(N,K0*RHO(X,Y))-(AN(N)*HANKP(2,N,K0*RHO(X,Y))))*
     $       DCOS(N*PHI(X,Y))

C  DETERMINE IF SUM IS CONVERGING  ****

	  IF (CDABS(TERM/(P+TERM)).GT.(1E-3)) THEN

C  NOT YET, ADD IN NEW VALUE  ****

	    P=P+TERM
	    NTERM=0
	  ELSE
	    NTERM=NTERM+1

C  IF CONVERGED FOR 3 CONSECUTIVE ORDERS, THEN YOU'RE FINISHED

	    IF(NTERM.EQ.3) GO TO 40
	      P=P+TERM
	  ENDIF

 30     CONTINUE
	WRITE (6,*) 'DID NOT CONVERGE, RHO= ',RHO(X,Y),' PHI= ',PHI(X,Y)

C  MULTIPLY FINAL TERM BY CONSTANTS   
      
 40     CONTINUE
	P = (-K0*H0/(J*W*EPS0))*P
      
      ENDIF
      
 50   CONTINUE

      RETURN
      END
*--------------------------------------------------------------------*
      SUBROUTINE ERHO(X,Y,R,AN,BN)
*--------------------------------------------------------------------*
      IMPLICIT NONE
      
      EXTERNAL BESSJ,BESSJP,HANK,HANKP,RHO,PHI
      
      REAL*8 EPS0,H0,W,X,Y,RADIUS,RHO,PHI,A
      
      COMPLEX*16 TERM,BESSJ,BESSJP,HANK,HANKP,K0,K1
      COMPLEX*16 AN(115),BN(115)
      COMPLEX*16 EPSRC,XMURC,R,J
      
      INTEGER N,NTERM,MODE
      
      PARAMETER (J=(0.0,1.0))
      PARAMETER (EPS0=8.841941E-12)
      
      COMMON /PAR/ RADIUS,K0,K1,W,EPSRC,XMURC,MODE

      A=RADIUS
      R=(1.0e-4,1.0e-4)
      H0=1.0/377.0
      
C  DO FIELDS INSIDE CYLINDER   
      
      IF (RHO(X,Y).LE.RADIUS) THEN

C  IF PEC (MODE=0) SKIP CALCULATING FIELDS INSIDE   
      
	IF (MODE.EQ.0) GO TO 50

	IF (RHO(X,Y).LT.1.0E-3) GO TO 50

	NTERM=0
      
	DO 10 N=1,115
	  IF (CDABS(BN(N)).GT.0.0) GOTO 888
c	    write(*,*) 'computing bn for term ',n
	    BN(N) = ((BESSJP(N,K0*A)*HANK(2,N,K0*A)-
     $        BESSJ(N,K0*A)*HANKP(2,N,K0*A))
     $        /(CDSQRT(XMURC/EPSRC)*BESSJP(N,K1*A)*HANK(2,N,K0*A)-
     $        BESSJ(N,K1*A)*HANKP(2,N,K0*A)))
c      write(*,*) 'bn= ',bn(n)
888       CONTINUE

	  TERM=((-J)**N)*N*BN(N)*BESSJ(N,K1*RHO(X,Y))*
     $         DSIN(N*PHI(X,Y))

C  DETERMINE IF SUM IS CONVERGING

	   IF (CDABS(TERM/(R+TERM)).GT.(1E-3)) THEN

C  NOT YET, ADD IN NEW VALUE
	     R=R+TERM
	     NTERM=0
	   ELSE
	     NTERM=NTERM+1

C  IF CONVERGED FOR 3 CONSECUTIVE ORDERS, THEN YOU'RE FINISHED

	     IF(NTERM.EQ.3) GO TO 20
	     R=R+TERM
	   ENDIF

 10   CONTINUE
      WRITE (6,*) 'DID NOT CONVERGE, RHO= ',RHO(X,Y),' PHI= ',PHI(X,Y)

C  MULTIPLY FINAL TERM BY CONSTANTS   
      
 20   CONTINUE
      R = ((-2.0*H0)/(J*W*EPS0*EPSRC*RHO(X,Y)))*(R)

      ELSE

C  DO TERMS OUTSIDE CYLINDER   

	NTERM=0
      
	DO 30 N=1,115
	  IF (CDABS(AN(N)).GT.0.0) GOTO 999
c	    write(*,*) 'computing an term for term ',n
            IF (MODE.EQ.0) THEN
	      AN(N)=BESSJP(N,K0*A)/HANKP(2,N,K0*A)
            ELSE IF (MODE.EQ.1) THEN
	      AN(N)=(CDSQRT(XMURC/EPSRC)*BESSJ(N,K0*A)*BESSJP(N,K1*A)-
     $        BESSJP(N,K0*A)*BESSJ(N,K1*A))/(CDSQRT(XMURC/EPSRC)*
     $        BESSJP(N,K1*A)*HANK(2,N,K0*A)-BESSJ(N,K1*A)*
     $        HANKP(2,N,K0*A))
            ENDIF
c	    write(*,*) 'an= ',an(n)
999       CONTINUE

	  TERM=((-J)**N)*2.0*N*(BESSJ(N,K0*RHO(X,Y))-
     $     (AN(N)*HANK(2,N,K0*RHO(X,Y))))*DSIN(N*PHI(X,Y))

C  DETERMINE IF SUM IS CONVERGING

	   IF (CDABS(TERM/(R+TERM)).GT.(1E-3)) THEN

C  NOT YET, ADD IN NEW VALUE

	     R=R+TERM
	     NTERM=0
	   ELSE
	     NTERM=NTERM+1

C  IF CONVERGED FOR 3 CONSECUTIVE ORDERS, THEN YOU'RE FINISHED

	     IF(NTERM.EQ.3) GO TO 40
	     R=R+TERM
	   ENDIF

 30   CONTINUE
      WRITE (6,*) 'DID NOT CONVERGE, RHO= ',RHO(X,Y),' PHI= ',PHI(X,Y)

C  MULTIPLY FINAL TERM BY CONSTANTS
      
 40   CONTINUE
      R = ((-H0)/(J*W*EPS0*RHO(X,Y)))*(R)
      
      ENDIF

 50   CONTINUE

      RETURN
      END      
*--------------------------------------------------------------------*
      REAL*8 FUNCTION RHO(X,Y)
*--------------------------------------------------------------------*
      IMPLICIT NONE
      REAL*8 X,Y

      RHO=DSQRT((X**2)+(Y**2))

      RETURN
      END
*--------------------------------------------------------------------*
      REAL*8 FUNCTION PHI(X,Y)
*--------------------------------------------------------------------*
      IMPLICIT NONE
      REAL*8 X,Y

      IF (X.EQ.0.0.AND.Y.EQ.0.0) THEN
	PHI=0.0
      ELSE
	PHI=DATAN2(Y,X)
      ENDIF

      RETURN
      END
