C**********************************************************************
      PROGRAM FDTD_1D_ver_1_00
C**********************************************************************
C
C   Aug. 27, 1998  Written by M.Fujii
C
C   This program calculates 1-D plane wave propagating in z-direction.
C
C
C   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C   !!!  IMPEDANCE BOUNDARY TEST !!!
C   !!!       AUG.27, 1998       !!!
C   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
C
C   [VARIABLES AND PARAMETERS]
C     KDMAX: MAXIMUM DIMENSION IN Z-DIRECTION
C     V:   E- AND H-FIELD
C     X:   EXICITATION POSITION
C     ER:  RELATIVE PERMITTIVITY
C     MR:  RELATIVE PERMEABILITY
C     SG:  ELECTRIC CONDUCTIBITY
C     KD:  DIMENSION (NUMBER OF YEE'S GRIDS * 2) IN Z-DIRECTION
C     TE:  TIME SERIES DATA OF E-FIELDS
C     TX:  TIME SERIES DATA OF EXCITATION ST
C     FC:  CENTER FREQUENCY
C     FS:  FREQUENCY SPAN
C     FR:  FREQUENCY RESOLUTION
C     TT:  PULSE DURATION
C     FTGT: TARGET FREQUENCY FOR EXCITATION
C     FNQ: NYQUIST FREQUENCY
C
C ------------------------------------------------------------------
      PARAMETER (KDMAX=300)
      PARAMETER (NNMAX=16384)
      PARAMETER (NIMP0=2)
      COMMON /BASICS/ UNIT,CVACUUM,PI
      COMMON /FIELDS/ V(0:KDMAX), X(0:KDMAX)
      COMMON /DOMAIN/ KL,KLP,KLP2,KU,KUM,KUM2,KD
      COMMON /TMCOND/ DT,NMAX,NPULSE,TT,NPRNT,NSTORE,NSTORE2,DTREAL
      COMMON /SPCOND/ DZ
      COMMON /TMDATA/ TE(NNMAX),TX(NNMAX)
      COMMON /DFTCND/ FC,FS,FR,PSD(NNMAX)
C
      REAL        V0(0:KDMAX)   ! FIELD VALUE AT PREVIOUS TIME STEP
      INTEGER     KIMP(NIMP0)   ! IMPEDANCE BOUNDARIES POSITIONS
      CHARACTER*1 DIMP(NIMP0)   ! DIRECTION OF IMPEDANCE BOUNDAIES
      REAL        ZIMP(NIMP0)   ! NORMALIZED IMPEDANCE VALUE
C
C     CHARACTER*3 SWIMP
      REAL MR
      DATA KL,KLP,KLP2 /0,1,2/
      DATA Q /0.25/
C
      WRITE(6,'(''================================================'')')
      WRITE(6,'(''               1-D FDTD PROGRAM                 '')')
      WRITE(6,'(''                    VER. 1                      '')')
      WRITE(6,'(''        WRITTEN BY M.FUJII. AUG.27, 1998        '')')
      WRITE(6,'(''================================================'')')
C
      PI = 2.0*ASIN(1.0)
      CVACUUM = 2.998E8
C
C----------------------------------------------------------------------
C   SET CALCULATION CONDITIONS
C----------------------------------------------------------------------
C
C   CALCULATION GEOMETRIES
C
      KD   = 2*32         !Space division with respect to 2*Yee's cells
      KU   = KL+KD
      KUM  = KU-1
      KUM2 = KU-2
C
C   LENGTH OF STRUCTURE
C
      X1 = 0.5
C
C   VARIABLE CONVERSION
C   (NORMALIZED PERMITTIVITY AND PERMEABILITY OF FREE SPACE: E0=M0=1)
C
      ER = 1.0                     !Rerative permittivity
      MR = 1.0                     !Rerative permeavility
      SG = 0.0                     !Normalized conducttivity
C
C   EXCITATION POSITION (DEFAULT VALUE = 0.0)
C
      X(KL+KD/2) = 1.0             !Excitation occure where X = 1.0.
C
C   TIME SERIES E-FILED DETECTION POSITION (SPECIFY E-NODE)
C
      KDTCT = KL+KD/2
C
C   IMPEDANCE BOUNDARY POSITIONS (SPECIFY E-NODE) 
C   NUMBER OF IMPEDANCE BOUNDARIES
C
      NIMP = 2
C
      KIMP(1) = KU
      DIMP(1) = 'U'
*     ZIMP(1) = SQRT(MR/ER)    !NORMALIZED IMPEDANCE (1.0)
      ZIMP(1) = 0.2
C
      KIMP(2) = KL
      DIMP(2) = 'L'
*     ZIMP(2) = SQRT(MR/ER)    !NORMALIZED IMPEDANCE (1.0)
      ZIMP(2) = 0.2
C
C   TIME AND SPACE CONDITIONS
C
      UNIT = 1.0                   !Unit length
C
      DZ      = X1/(KD/2)          !Length of Yee's cell for FDTD
      DT      = 0.8*DZ             !Time discretization for FDTD
      DTREAL  = DT*UNIT/CVACUUM
      NMAX    = INT(2.0/DT)        !Maximum iterations
      NPRNT   = (NMAX-1)/10+1      !Interval of printing fields
      NSTORE  = 1                  !Interval of storing field distributions
      NSTORE2 = 1                  !Interval of storing time signal
      IF (NMAX/NSTORE2.GT.NNMAX) STOP "!! NMAX/NSTORE2 > NNMA !!"
C
C   TARGET FREQUENCY
C
      FTGT = 1.0                   !Center freq. of excitation pulse
      FNQ  = 1.0/(2.0*DT)          !Nyquist freq.
C
C   EXCITATION PULSE DURATION
C
      TT = 1.0/FTGT                !Time period
      NPULSE = INT(0.5*TT/DT)      !Excitation pulse duration
C
C   DISCRETE FOURIER TRANSFORMATION 
C
      FC = 1.0                     !Center frequency in DFT
      FS = 0.1                     !Frequency span in DFT
      FR = 0.001                   !Frequency resolution in DFT
C
C----------------------------------------------------------------------
C   This is the end of user definition.
C----------------------------------------------------------------------
C
C   --------------------------------
C   PRINT OUT CALCULATION CONDITIONS
C   --------------------------------
C
      WRITE (6,6100) KD
      WRITE (6,6110) ER,MR,SG
      WRITE (6,6120) UNIT,DZ,DT,NMAX,NPRNT,NSTORE,NSTORE2
      WRITE (6,6130) FTGT,TT,NPULSE,FNQ
      WRITE (6,6140) FC,FS,FR
 6100 FORMAT ('[CALCULATION CONDITIONS]'
     +       /'  GEOMETRIC CONDITIONS'
     +       /'                       KD =',I5)
 6110 FORMAT ('  MATERIALS'
     +       /'                       ER =',G15.7,
     +       /'                       MR =',G15.7,
     +       /'                       SG =',G15.7)
 6120 FORMAT ('  TIME CONDISIONS'
     +       /'                     UNIT =',G15.7,
     +       /'                       DZ =',G15.7,
     +       /'                       DT =',G15.7,
     +       /'                     NMAX =',I8,
     +       /'                    NPRNT =',I8,
     +       /'                   NSTORE =',I8,
     +       /'                  NSTORE2 =',I8)
 6130 FORMAT ('  EXCITATION CONDITIONS'
     +       /'                     FTGT =',G15.7,
     +       /'                       TT =',G15.7,
     +       /'                   NPULSE =',I8,
     +       /'      NYQUIST FREQ.   FNQ =',G15.7)
 6140 FORMAT ('  FOURIER TRANSFORMATION'
     +       /'      CENTER FREQ.     FC =',G15.7,
     +       /'      FREQ. SPAN       FS =',G15.7,
     +       /'      FREQ. RESOLUTION FR =',G15.7)
C
C   -----------------------------------------------------
C
C   SET EXCITATION POSITION
C
      WRITE (6,6020) 'EXCITATION POSITION AND DISTRIBUTION',
     +                (INT(ABS(10.0*X(K))),K=KL,KU)
 6020 FORMAT (1H ,A,/,1H ,240I1)
C
C   SET COEFFICIENTS
C
      DNM = 2.0*ER+SG*DT
      C0  =(2.0*ER-SG*DT)/DNM
      CE  = 2.0*DT/DNM/DZ
      CS  = 2.0*DT/DNM
      CM  = DT/MR/DZ
C
C   ===============================
C   START TIME MARCHING ARGOLITHM
C   ===============================
C
      PRINT *,'START TIME MARCHING ARGOLITHM'
      DO 10000 N=1,NMAX
C
C     RAISED COSINE PULSE EXCITATION
C
      IF (N.LE.NPULSE) THEN
        ST = SIN(PI*FLOAT(N)/FLOAT(NPULSE))**2
      ELSE
        ST = 0.0
      ENDIF
C
C   ================   
C   CALCULATE FIELDS
C   ================   
C
C   -----------------------
C   E-FIELD INTERNAL POINTS
C   -----------------------
C
      DO 100 K=KLP2,KUM2,2
        V(K) = C0*V(K)-CE*(V(K+1)-V(K-1))-CS*ST*X(K)
 100  CONTINUE
C
C   ----------------
C   E-FIELD BOUNDARY 
C   ----------------
C
C   PERFECT ELECTRIC CONDUCTOR (PEC) BOUNDARY AT E-NODES
C
      V(KL) = 0.0
      V(KU) = 0.0
C
C   -----------------------
C   H-FIELD INTERNAL POINTS
C   -----------------------
C
      DO 200 K=KLP,KUM,2
        V(K) = V(K)-CM*(V(K+1)-V(K-1))
 200  CONTINUE
C
C   ----------------
C   H-FIELD BOUNDARY
C   ----------------
C
C   IMPEDANCE BOUNDARY
C
      DO 210 M=1,NIMP
        K = KIMP(M)
        IF (DIMP(M).EQ.'U') THEN
          ANM = 4.0*MR*DZ-3.0*ZIMP(M)*DT
          DNM = 4.0*MR*DZ+3.0*ZIMP(M)*DT
          CI0 = ANM/DNM
          CI1 = DT/DNM
          V(K-1) = CI0*V0(K-1)
     +            +CI1*(4.0*V(K-2)+ZIMP(M)*(V(K-3)+V0(K-3)))
        ELSEIF (DIMP(M).EQ.'L') THEN
          ANM = 4.0*MR*DZ-3.0*ZIMP(M)*DT
          DNM = 4.0*MR*DZ+3.0*ZIMP(M)*DT
          CI0 = ANM/DNM
          CI1 = DT/DNM
          V(K+1) = CI0*V0(K+1)
     +            -CI1*(4.0*V(K+2)-ZIMP(M)*(V(K+3)+V0(K+3)))
        ENDIF
 210  CONTINUE
C
C   ========================
C   END OF FIELD CALCULATION
C   ========================
C
C   ------------
C   PRINT FIELDS
C   ------------
C     
      IF (N/NPRNT*NPRNT.EQ.N.OR.N/NSTORE*NSTORE.EQ.N) THEN
        CALL PRNT(N)
      END IF
C
      IF (N/NSTORE2*NSTORE2.EQ.N) THEN
        CALL VIDATA(N,ST,KDTCT)
      END IF
C
C   ----------------------------------------
C   STORE FIELD VALUES AT PREVIOUS TIME STEP
C   ----------------------------------------
C
      DO 80 K=KL,KU
        V0(K) = V(K)
 80   CONTINUE
C
C   ==============================
C   END OF TIME MARCHING ARGORITHM
C   ==============================
C
10000 CONTINUE
C
C   POSTPROCESSING DATA
C
      CALL PSTPRCS
C
      STOP
      END

C**********************************************************************
      SUBROUTINE PRNT(N)
C**********************************************************************
      PARAMETER (KDMAX=300)
      COMMON /BASICS/ UNIT,CVACUUM,PI
      COMMON /FIELDS/ V(0:KDMAX), X(0:KDMAX)
      COMMON /DOMAIN/ KL,KLP,KLP2,KU,KUM,KUM2,KD
      COMMON /TMCOND/ DT,NMAX,NPULSE,TT,NPRNT,NSTORE,NSTORE2,DTREAL
      COMMON /SPCOND/ DZ
      INTEGER         IE(0:KDMAX), IH(0:KDMAX)
      CHARACTER*1     B (0:KDMAX)
C
      IF (N/NPRNT*NPRNT.EQ.N) THEN
C
C   INITIALIZATION
C
        DO 10 K=KL,KU
          IE(K) = 0
          IH(K) = 0
          B(K) = '-'
 10     CONTINUE
        DO 20 K=KL,KU,5
          B(K) = '+'
 20     CONTINUE
        DO 30 K=KL,KU,10
          B(K) = '0'
 30     CONTINUE
C
C   SEARCH MAXIMUM FIELD VALUE
C
        EMAX = 0.0
        HMAX = 0.0
        DO 200 K=KL,KU,2
          EMAX = AMAX1(EMAX,ABS(V(K)))
 200    CONTINUE
        DO 210 K=KLP,KUM,2
          HMAX = AMAX1(HMAX,ABS(V(K)))
 210    CONTINUE
        IF (EMAX.EQ.0.0) THEN
          EMAX = 1.0
          PRINT *,'NO E-FIELDS ARE DETECTED.'
        ENDIF
C
C   CONVERT FIELD VALUES INTO INTEGER NUMBER
C
        DO 300 K=KL,KU,2
          IE(K) = INT(10.0*ABS(V(K))/EMAX)
 300    CONTINUE
        DO 310 K=KLP,KUM,2
          IH(K) = INT(10.0*ABS(V(K))/HMAX)
 310    CONTINUE
C
C   PRINT FIELD DISTRIBUTIONS ON DISPLAY
C
        WRITE (6,6000) 'TIME ITTERATIONS N=',N,
     +                 'E-FIELD (MAX:',EMAX,')'
        WRITE (6,6010) (IE(K),K=KL,KU)
        WRITE (6,6020) 'H-FIELD (MAX:',HMAX,')'
        WRITE (6,6010) (IH(K),K=KL,KU)
        WRITE (6,6030) (B(K),K=KL,KU)
 6000   FORMAT (A,I5,/,1H ,A,G10.3,A)
 6010   FORMAT (1H ,240I1)
 6020   FORMAT (1H ,A,G10.3,A)
 6030   FORMAT (1H ,240A)
C
      ENDIF
C
C   STORE TIME VARING FIELDS 
C
      IF (N/NSTORE*NSTORE.EQ.N) THEN
C
        OPEN (10,FILE='TIMEEFLD')
        OPEN (11,FILE='TIMEHFLD')
C
C   HEADER
C
        IF (N.EQ.NSTORE) WRITE (10,6100)
 6100   FORMAT ('# Time varing electric field distribution. ',
     +         /'#   TIME',14X,'Z',14X,'E')
        IF (N.EQ.NSTORE) WRITE (11,6110)
 6110   FORMAT ('# Time varing magnetic field distribution. ',
     +         /'#   TIME',14X,'Z',14X,'H')
C
C   STORE FIELD DISTRIBUTION DATA
C
        DO 400 K=KL,KU,2
          WRITE (10,6120) DT*FLOAT(N),0.5*DZ*FLOAT(K),V(K)
 400    CONTINUE
        WRITE (10,*)
        DO 410 K=KLP,KUM,2
          WRITE (11,6120) DT*FLOAT(N),0.5*DZ*FLOAT(K),V(K)
 410    CONTINUE
        WRITE (11,*)
 6120   FORMAT (1H ,G15.7,3X,G15.7,3X,G15.7)
C
      ENDIF
C
      RETURN
      END

C**********************************************************************
      SUBROUTINE VIDATA(N,ST,KDTCT)
C**********************************************************************
      PARAMETER (KDMAX=300)
      PARAMETER (NNMAX=16384)
      COMMON /BASICS/ UNIT,CVACUUM,PI
      COMMON /FIELDS/ V(0:KDMAX),X(0:KDMAX)
      COMMON /DOMAIN/ KL,KLP,KLP2,KU,KUM,KUM2,KD
      COMMON /TMCOND/ DT,NMAX,NPULSE,TT,NPRNT,NSTORE,NSTORE2,DTREAL
      COMMON /SPCOND/ DZ
      COMMON /TMDATA/ TE(NNMAX),TX(NNMAX)
C
      NS = N/NSTORE2
C
C   STORE EXCITATION 
C
      TX(NS) = ST
C
C   STORE E-FILED
C
      TE(NS) = V(KDTCT)
C
      RETURN
      END

C**********************************************************************
      SUBROUTINE PSTPRCS
C**********************************************************************
      PARAMETER (KDMAX=300)
      PARAMETER (NNMAX=16384)
      COMMON /BASICS/ UNIT,CVACUUM,PI
      COMMON /FIELDS/ V(0:KDMAX),XJ(0:KDMAX)
      COMMON /DOMAIN/ KL,KLP,KLP2,KU,KUM,KUM2,KD
      COMMON /TMCOND/ DT,NMAX,NPULSE,TT,NPRNT,NSTORE,NSTORE2,DTREAL
      COMMON /SPCOND/ DZ
      COMMON /TMDATA/ TE(NNMAX),TX(NNMAX)
      COMMON /DFTCND/ FC,FS,FR,PSD(NNMAX)
C
C   ----------------------
C   STORE TIME SERIES DATA
C   ----------------------
C 
      OPEN (12,FILE='TIMEDATA')
C
      WRITE(12,6000)
      DO 100 N=1,NMAX/NSTORE2
C       WRITE (12,6010) DTREAL*FLOAT(N*NSTORE2),TX(N),TE(N)
        WRITE (12,6010) DT*FLOAT(N*NSTORE2),TX(N),TE(N)
 100  CONTINUE
 6000 FORMAT ('# TIME SERIES DATA',
     +       /'# ',
     +       /'# TIME',10X,'EXCITE',10X,'E-FIELD')
 6010 FORMAT (1H ,G15.7,1X,G15.7,1X,G15.7)
C
C   ------------------------
C   SPECTRUM ANALYSIS BY DFT
C   ------------------------
C
      PRINT *,'PROCESSING SPECTRUM ANARYSIS BY DFT ...'
      N1 = INT(NPULSE/NSTORE2)
      N2 = INT(NMAX/NSTORE2)
      OMG1 = 2.0*PI*(FC-0.5*FS)
      OMGR = 2.0*PI*FR
      DTS  = DT*FLOAT(NSTORE2)
      DO 200 M=1,INT(FS/FR)+1
        OMG = OMG1+OMGR*FLOAT(M-1)
        R = 0.0
        X = 0.0
        DO 210 N=N1,N2
C   DFT
          R = R+TE(N)*COS(OMG*FLOAT(N)*DTS)*DTS
          X = X-TE(N)*SIN(OMG*FLOAT(N)*DTS)*DTS
 210    CONTINUE
      PSD(M) = R*R+X*X
 200  CONTINUE
C
C   SEARCH MAXMUM
C
      PSDMAX = 0.0
      F1 = FC-0.5*FS
      DO 300 M=1,INT(FS/FR)+1
        IF (PSD(M).GT.PSDMAX) THEN
          PSDMAX = PSD(M)
          F0 = F1+FR*FLOAT(M-1)
        ENDIF
 300  CONTINUE
C
      WRITE (6,6050) 'FREQUENCY RESOLUTION =',FR,
     +               'RESONANT FREQUENCY   =',F0
 6050 FORMAT (1H ,A,G15.7,/1H ,A,G15.7)
C
C   STORE POWER SPECTRUM DENSITY
C
      OPEN (13,FILE='FREQDATA')
      WRITE (13,6100)
      F1 = FC-0.5*FS
      DO 400 M=1,INT(FS/FR)+1
        FREQ = F1+FR*FLOAT(M-1)
        WRITE (13,6110) FREQ,10.0*LOG10(PSD(M))
 400  CONTINUE
 6100 FORMAT ('# FREQENCY DOMAIN DATA',
     +       /'# ',
     +       /'# FREQ.',9X,'PSD')
 6110 FORMAT (1H ,G15.7,1X,G15.7)
C
      RETURN
      END
















