C*********************************************************************
C
C  2DFDTD-BASIC
C
C  Two-Dimensional Time-Domain Finete-Difference Analysis Program
C
C  Written by M.Fujii, July  6, 1998
C  Completed,          July  9, 1998
C  Modified for FDTD,  July 31, 2000
C  Corrected,          May  20, 2004
C
C  Features:
C  -- TE wave propagating on xz plane 
C     (Structures are uniform in y-direction.)
C  -- Rectangular uniform mesh
C  -- Homogeneous lossy dielectric material
C  -- No absorbing boundary conditions (to be implemented later)
C  -- 2D FDTD only
C
C  PARAMETERS
C
C     IN,KN : MAXIMUM CALCULATION REGION IN X AND Z DIRECTIONS
C     NP    : NUMBER OF PORTS FOR TIME SIGNAL EXTRACTION
C     NN    : MAXIMUM NUMBER OF TIME SERIES DATA
C     NF    : NUMBER OF FIELDS (EY,HX,HZ)
C     NS    : NUMBER OF FIELD DISTRIBUTION SNAPSHOTS (MAX 99)
C
      PARAMETER (IN=300, KN=300)
      PARAMETER (NP=3)
      PARAMETER (NN=10000)
      PARAMETER (NF=3)
      PARAMETER (NS=10)
C
C  CALCULATION REGION PARAMETERS
C     DEFINITION :
C     |----+----+----+----+----+- - ---+----+----+----+----|  X(I)
C     IL  ILP  ILP2 ...  KLP4 ...     IUM4 ...  IUM2  IUM  IU  
C     <--------->
C      YEE'S CELL
C
      INTEGER
     +    IL, ILP, ILP2, IU, IUM, IUM2, ID,
     +    KL, KLP, KLP2, KU, KUM, KUM2, KD
C
C  MAIN DIMENSION
C     V    : FIELD COMPONENTS
C     SJ   : CURRENT SORCE FOR HAAR BASIS 
C
      REAL 
     +    V (IN,KN),   SJ(IN,KN),  PC(IN,KN)
C
C  UPDATE EQUATION COEFFICIENTS
C
      REAL
     +    C0, CE, CMX, CMZ
C
C  CALCULATION CONDITIONS
C
      REAL
     +    DT, UNIT, CVACUUM
      INTEGER
     +    NMAX, NSTORE, NSTORE2, NPRNT, NMOD, NPD
C
C  OUTPUT DATA
C
      REAL
     +    FLDOUT(0:IN+1,0:KN+1,NF,NS), TMOUT(0:NN,0:NP)
C
C  FIELD DISTRIBUTION FILE NAME
C
      CHARACTER*9
     +    FL(NF,NS)
C
C  PORT POSITIONS
C
      INTEGER
     +    IPORTL(NP), IPORTU(NP), KPORTL(NP), KPORTU(NP)
C
C  FOR PRINTING FIELD
C
      REAL
     +    FLDPRNT(IN,KN)
C
C  INITIAL DATA
C
      DATA IL,KL       /3,3/
      DATA PI, CVACUUM /3.14159265, 2.998E8/
C
      WRITE (6,6000)
 6000 FORMAT (
     + /' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
     + /' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
     + /' !!!!                                               !!!!'
     + /' !!!!                    2D-FDTD                    !!!!'
     + /' !!!!                 BASIC PROGRAM                 !!!!'
     + /' !!!!                                               !!!!'
     + /' !!!!   2D Time-Domain Finite-Difference Analysis   !!!!'
     + /' !!!!                                               !!!!'
     + /' !!!!               Written by M.Fujii              !!!!'
     + /' !!!!                  July 31, 2000                !!!!'
     + /' !!!!                                               !!!!'
     + /' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
     + /' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
C
C  ********************************
C  INPUT THE CALCULATION CONDITIONS
C  ********************************
C
C  -------------------------------------------------
C  (1) SET UNIT IN LENGTH AND OTHER BASIC PARAMETERS
C  -------------------------------------------------
C
      UNIT = 1.0E-3                
      GHZ  = 1.0E9*UNIT/CVACUUM
      Z0   = 120*PI
C
C  --------------------------
C  (2) SET CALCULATION REGION
C  -------------------------- 
C
      XMAX0 = 100.0
      ZMAX0 = 100.0
C
C  ------------------------------------------------------------------
C  (3) SET NUMBERS OF GRIDS (EVEN NUMBERS, 2 X NUMBER OF YEE'S CELLS)
C  --------------------------============--==========----------------
C
      ID = 100
      KD = 100
C
C  GRID SPACINGS (1/2 OF YEE'S CELL)
C
      DXH = XMAX0/FLOAT(ID)
      DZH = ZMAX0/FLOAT(KD)
C
C  CALCULATION REGION PARAMETERS
C
      ILP  = IL+1
      ILP2 = IL+2
      ILP4 = IL+4
      KLP  = KL+1
      KLP2 = KL+2
      KLP4 = KL+4
      IU   = IL+ID
      IUM  = IU-1
      IUM2 = IU-2
      IUM4 = IU-4
      KU   = KL+KD
      KUM  = KU-1
      KUM2 = KU-2
      KUM4 = KU-4
C
C  ----------------------------
C  (4) SET DIELECTRIC CONSTANTS
C  ----------------------------
C
      ER = 1.0
C
C  ----------------------------------------------------
C  (5) SET CONDUCTIVITY 
C      SG SHOULD BE NORMALIZED BY MULTIPLYING (Z0*UNIT)
C  ----------------------------------------------------
C
C     SG = 1.0E-2*Z0*UNIT
      SG = 0.0
C
C  --------------------------------------------------------
C  (6) SET PEC/PMC BOUNDARY POSITIONS AND DIRECTION
C      THE NUMBERS SHOULD SPECIFY ELECTRIC FIELD (EY) NODES
C  --------------------------------------------------------
C
      DO 100 K=KL,KU
      DO 100 I=IL,IU
        PC(I,K) = 1.0
 100  CONTINUE

      DO 110 K=KL+30,KL+70
      DO 110 I=IL+30,IL+70
        PC(I,K) = 0.0
 110  CONTINUE
C
C  ---------------------------------------------------------
C  (10) SET OUTPUT PORT POSITIONS (MAXIMUM NP(=3) POINTS)
C       THE NUMBERS SHOULD SPECIFY ELECTRIC FIELD (EY) NODES
C  ---------------------------------------------------------
C
      IPORTL(1) = IL
      IPORTU(1) = IU
      KPORTL(1) = 10
      KPORTU(1) = 10
      IPORTL(2) = IL
      IPORTU(2) = IU
      KPORTL(2) = KU-10
      KPORTU(2) = KU-10
      IPORTL(3) = IL
      IPORTU(3) = IU
      KPORTL(3) = 10
      KPORTU(3) = 10
C
C  ---------------------------------------------
C  (11) SET EXCITATION POSITION (ON ACTUAL NODE)
C  ---------------------------------------------
C
      IEXCT1 = IL
      IEXCT2 = IU
      KEXCT1 = KL
      KEXCT2 = KL+10
C
C  SMOOTH EXCITATION DISTRIBUTION (RAISED-COSINE DISTRIBUTION)
C
      DO 500 K=KEXCT1,KEXCT2,2
      DO 500 I=IEXCT1,IEXCT2,2
        SJ(I,K) = SIN(PI*FLOAT(I-IEXCT1)/FLOAT(IEXCT2-IEXCT1))**2 
     +           *SIN(PI*FLOAT(K-KEXCT1)/FLOAT(KEXCT2-KEXCT1))**2 
 500  CONTINUE
C
C  ---------------------------------------------------
C  (12) SET TIME INTERVAL
C        DTMRGN : MARGIN FOR DT TO MAXIMUM LIMIT OF DT
C        DT     : NORMALIZED TIME DISCRETIZATION
C        DTREAL : UNNORMALIZED VALUE OF DT IN SECOND
C  ---------------------------------------------------
C
      DTMRGN = 0.98
C
      DD = 1.0/SQRT(1.0/(2.0*DXH)**2+1.0/(2.0*DZH)**2)
      DT     = DTMRGN*DD/SQRT(ER)
      DTREAL = DT*UNIT/CVACUUM
C
C  ---------------------------------------------------------------------
C  (13) SET EXCITING PULSE DURATION (SINE-MODULATED RAISED-COSINE PULSE)
C        FTGT : TARGET FREQUENCY
C        NMOD : NUMBER OF SINE MODULATION IN UNIT OF A HALF PERIOD
C        TT   : PULSE DURATION TIME
C        NPD  : NUMBER OF PULSE DURATION
C  ---------------------------------------------------------------------
C
      FTGT = 10.0*GHZ
      NMOD = 2
C
      TT   = 1.0/FTGT
      NPD  = INT(TT/DT)*NMOD
C
C  ----------------------------------------------------------
C  (14) SET TIME STEP
C        NMAX    : NUMBER OF TIME STEPS ( NMAX/NSTORE2 < NN )
C        NSTORE2 : TIME SERIES DATA STORE INTERVAL
C        NSTORE  : FIELD DISTRIBUTION DATA STORE INTERVAL       
C        NPRNT   : FIELD DATA PRINT INTERVAL
C        FNQ     : NYQUIST FREQUENCY (FOR REFERENCE)
C  ----------------------------------------------------------
C
      NMAX    = INT(4.42E-10/DTREAL)
      NSTORE2 = 1
C
      NSTORE  = NMAX/NS
      NPRNT   = NMAX/20
      FNQ     = 1.0/(2.0*DT*FLOAT(NSTORE2))
      IF (NMAX/NSTORE2.GT.NN) STOP 'DIMENSION OVER! NMAX/NSTORE2 > NN'
C
C  DISPLAY CALCULATION CONDITIONS
C
      WRITE (6,6010) ID/2,KD/2,2.0*DXH,2.0*DZH,UNIT
      WRITE (6,6020) DT,DTREAL,DTMRGN,FTGT/GHZ,FTGT,FNQ,
     +               TT,NMOD,NPD,
     +               NMAX,NSTORE,NPRNT,NSTORE2
C
 6010 FORMAT (/'[2D TD-MRA CALCULATION CONDITONS]'
     +      /,' A.MODEL GEOMETRIES'
     +      /,'   YEE''S GRID COUNTS FOR' 
     +      /,'     X,Z-DIRECTION    = ',2I6
     +      /,'   YEE''S GRID DIMENSIONS FOR'
     +      /,'     X-DIRECTION                      = ',G15.5
     +      /,'     Z-DIRECTION                      = ',G15.5
     +      /,'   LENGTH UNIT IN M                   = ',G15.5)
 6020 FORMAT (' B.TIME SPACE CONDITIONS'
     +      /,'   TIME INTERVAL, DT                  = ',G15.5
     +      /,'   TIME INTERVAL /SEC, DTREAL         = ',G15.5
     +      /,'   TIME STEP MARGIN, DTMRGN           = ',G15.5
     +      /,'   TERGET FREQUENCY /GHZ, FTGT        = ',G15.5
     +      /,'   TERGET FREQUENCY, FTGT             = ',G15.5
     +      /,'   NYQIST FREQUENCY, FNQ              = ',G15.5
     +      /,'   PULSE DURATION TIME ,TT            = ',G15.5
     +      /,'   NUMBER OF HALF-SINE MUDULATION     = ',I8
     +      /,'   PULSE DURATION /IT, NPD            = ',I8
     +      /,'   MAXIMUM NUMBER OF TIME STEPS, NMAX = ',I8
     +      /,'   FIELD DATA STORE INTERVAL, NSTORE  = ',I8
     +      /,'   FIELD DATA PRINT INTERVAL, NPRNT   = ',I8
     +      /,'   TIME DATA STORE INTERVAL, NSTORE2  = ',I8)
C
C  ***************************
C  PREPROCESSING OF INPUT DATA
C  ***************************
C
C  UPDATE-EQUATION COEFFICIENTS
C
      C0 = (2.0*ER-SG*DT)/(2.0*ER+SG*DT)
      CE = 2.0*DT/(2.0*ER+SG*DT)/(4.0*DXH*DZH)
      CMX = -DT*DXH/DZH
      CMZ = -DT*DZH/DXH
C
C  *****************************
C  *****************************
C  START TIME MARCHING ALGORITHM
C  *****************************
C  *****************************
C
      PRINT *,'SOLVING UPDATE EQUATIONS ...'
C
      DO 10000 IT=1,NMAX
C
        IF (IT.LT.NPD) THEN
          ST = SIN(PI*FLOAT(IT-1)/FLOAT(NPD))**2
     +        *SIN(FLOAT(NMOD)*PI*FLOAT(IT-1)/FLOAT(NPD))
        ELSE
          ST = 0.0
        END IF
C
C  *************************************************
C  SOLVE UPDATE EQUATIONS FOR MAXWELL CURL EQUATIONS
C  *************************************************
C
C  --------------------------------
C  EY NODE UPDATE (INTERNAL POINTS)
C  --------------------------------
C
      DO 1000 K=KLP2,KUM2,2
      DO 1000 I=ILP2,IUM2,2
        V(I,K) = PC(I,K)*(C0*V(I,K)
     +                   +CE*(V(I,K+1)-V(I,K-1)-V(I+1,K)+V(I-1,K)
     +                   -ST*SJ(I,K)))
 1000 CONTINUE
C
C  ---------------------------------------
C  HX AND HZ NODE UPDATE (INTERNAL POINTS)
C  ---------------------------------------
C
C  HX NODE
C
        DO 1620 K=KLP,KUM,2
        DO 1620 I=IL,IU,2
          V(I,K) = PC(I,K)*(V(I,K)+CMX*(-V(I,K+1)+V(I,K-1)))
 1620   CONTINUE
C
C  HZ NODE
C
        DO 1640 K=KL,KU,2
        DO 1640 I=ILP,IUM,2
          V(I,K) = PC(I,K)*(V(I,K)+CMZ*( V(I+1,K)-V(I-1,K)))
 1640   CONTINUE
C
C  ****************
C  END OF EQUATIONS
C  ****************
C
C  ************************************
C  FIELD DISTRIBUTION SNAPSHOTS STORAGE
C  ************************************
C
        IF (IT/NSTORE*NSTORE.EQ.IT) THEN
C
          IS = IT/NSTORE
C
C  --------
C  EY FIELD
C  --------
C
          DO 3000 K=KL,KU,2
          DO 3000 I=IL,IU,2
            FLDOUT(I,K,1,IS) = FLDOUT(I,K,1,IS)+V(I,K)
 3000    CONTINUE
C
C  --------
C  HX FIELD
C  --------
C
          DO 3200 K=KLP,KUM,2
          DO 3200 I=IL,IU,2
            FLDOUT(I,K,2,IS) = FLDOUT(I,K,2,IS)+V(I,K)
 3200     CONTINUE
C
C  --------
C  HZ FIELD
C  --------
C
          DO 3400 K=KL,KU,2
          DO 3400 I=ILP,IUM,2
            FLDOUT(I,K,3,IS) = FLDOUT(I,K,3,IS)+V(I,K)
 3400     CONTINUE
C
        ENDIF
C
C  **************************************
C  EY FIELD DISTRIBUTION PRINT ON DISPLAY
C  **************************************
C
        IF (IT/NPRNT*NPRNT.EQ.IT) THEN
C
          DO 4000 K=KL,KU
          DO 4000 I=IL,IU
            FLDPRNT(I,K) = 0.0
 4000     CONTINUE
C
C         EY FIELD
C
          DO 4050 K=KL,KU,2
          DO 4050 I=IL,IU,2
              FLDPRNT(I,K) = FLDPRNT(I,K)+V(I,K)
 4050     CONTINUE
C
          VMAX = 0.0
          DO 4100 K=KL,KU
          DO 4100 I=IL,IU
            VMAX = AMAX1(VMAX,ABS(FLDPRNT(I,K)))
 4100     CONTINUE
          IF (VMAX.EQ.0.0) THEN
            VMAX = 1.0
          ENDIF
          WRITE (6,6050) IT,VMAX
 6050     FORMAT (//' IT = ',I5,/' MAX EY = ',G15.7/)
C
          DO 4200 I=IU+1,IL,-1
            WRITE (6,6100) (INT(10.0*ABS(FLDPRNT(I,K)/VMAX)),K=KL,KU)
 4200     CONTINUE
 6100     FORMAT (1H ,500I1)
C
        ENDIF
C
C  **************************************************
C  TIME SERIES DATA STORAGE (WITH DETECTING EY VALUE)
C  **************************************************
C
        IF (IT/NSTORE2*NSTORE2.EQ.IT) THEN
C
          ITS = IT/NSTORE2
C
          DO 5000 IP=1,NP
            DO 5100 I=IL+IPORTL(IP),IL+IPORTU(IP),2
            DO 5100 K=KL+KPORTL(IP),KL+KPORTU(IP),2
              TMOUT(ITS,IP) = TMOUT(ITS,IP)-V(I,K)
 5100       CONTINUE
            TMOUT(ITS,IP) = TMOUT(ITS,IP)
     +                   /FLOAT(ABS(IPORTU(IP)-IPORTL(IP))/2+1)
     +                   /FLOAT(ABS(KPORTU(IP)-KPORTL(IP))/2+1)
 5000     CONTINUE 
C
C  ----------------------
C  STORE EXCITATION PULSE
C  ----------------------
C
          TMOUT(ITS,0) = ST
C
        ENDIF
C
C  ******************************
C  ******************************
C  END OF TIME MARCHING ALGORITHM
C  ******************************
C  ******************************
C
10000 CONTINUE
C
C  *****************************
C  POSTPROCESSING OF OUTPUT DATA
C  *****************************
C  ------------------------------------------------------------
C  FIELD DISTRIBUTION NORMALIZATION AND DATA STORAGE INTO FILES
C  ------------------------------------------------------------
C
      PRINT *,'POST-PROCESSING ...'
C
C  CREATE FILE NAME
C
      DO 7600 IS=1,NS
        NDIGIT1 = 48+(IS-IS/10*10)
        NDIGIT2 = 48+IS/10
        FL(1,IS) = 'EY'//CHAR(NDIGIT2)//CHAR(NDIGIT1)
        FL(2,IS) = 'HX'//CHAR(NDIGIT2)//CHAR(NDIGIT1)
        FL(3,IS) = 'HZ'//CHAR(NDIGIT2)//CHAR(NDIGIT1)
 7600 CONTINUE
C
C  ------------------
C  FIELD DATA STORAGE
C  ------------------
C
      DO 7700 IS=1,NS  !NUMBER OF SNAPSHOTS
C
C       EY
C
        OPEN (11,FILE=FL(1,IS),FORM='FORMATTED')
        WRITE (11,6300) FL(1,IS)
        DO 7820 K=KL,KU,2
          Z = DZH*FLOAT(K)
          DO 7920 I=IL,IU,2
            X = DXH*FLOAT(I)
            WRITE (11,6400) Z,X,FLDOUT(I,K,1,IS)
 7920     CONTINUE
          WRITE (11,*)
 7820   CONTINUE
        CLOSE(11)
C
C       HX
C
        OPEN (11,FILE=FL(2,IS),FORM='FORMATTED')
        WRITE (11,6300) FL(2,IS)
        DO 7840 K=KLP,KUM,2
          Z = DZH*FLOAT(K)
          DO 7940 I=IL,IU,2
            X = DXH*FLOAT(I)
            WRITE (11,6400) Z,X,FLDOUT(I,K,1,IS)
 7940     CONTINUE
          WRITE (11,*)
 7840   CONTINUE
        CLOSE(11)
C
C       HZ
C
        OPEN (11,FILE=FL(3,IS),FORM='FORMATTED')
        WRITE (11,6300) FL(3,IS)
        DO 7860 K=KL,KU,2
          Z = DZH*FLOAT(K)
          DO 7960 I=ILP,IUM,2
            X = DXH*FLOAT(I)
            WRITE (11,6400) Z,X,FLDOUT(I,K,1,IS)
 7960     CONTINUE
          WRITE (11,*)
 7860   CONTINUE
        CLOSE(11)
C
 7700 CONTINUE
C
 6300 FORMAT ('# ',A3,' FIELD ON ZX PLANE'
     +       /'# Z',14X,'X',14X,'FIELD VALUE')
 6400 FORMAT (' ',3(G16.8,3X))
C
C  ----------------------
C  STORE TIME SERIES DATA
C  ----------------------
C
      NDATA    = NMAX/NSTORE2
      TIMESTEP = DT*FLOAT(NSTORE2)*UNIT/CVACUUM
C
      OPEN (10,FILE='TIMEDATA',FORM='FORMATTED')
      WRITE(10,6500)
 6500 FORMAT('# TD-MRA TIME SERIES DATA (IN MKSA UNIT)'
     +   /'# FOR S-PARAMETER EXTRACTION'
     +   /'     TIME(SEC)',5X,'EXCITATION',6X,'V1',14X,'V2',14X,'V3')
      DO 8000 N=0,NDATA
	TIME = TIMESTEP*FLOAT(N)
        WRITE(10,6510) TIME,(TMOUT(N,I),I=0,NP)
 8000 CONTINUE
 6510 FORMAT(1H ,5G16.8)
      CLOSE (10)
C
      STOP
      END



C*********************************************************************
      SUBROUTINE PRNT (N,M,A)
C*********************************************************************
C
C     Supplemental subroutine for checking values in an array A with 
C     dimensions N x M.
C
      DIMENSION A(N,M)
C
      VMAX = 0.0
      DO 1100 K=1,M
      DO 1100 I=1,N
        VMAX = AMAX1(VMAX,A(I,K))
 1100 CONTINUE
      IF (VMAX.EQ.0.0) THEN
        VMAX = 1.0
      ENDIF
      WRITE (6,6050) N,M,VMAX
 6050 FORMAT (//' ',I5,' X ',I5,' MATRIX, WITH MAX VALUE = ',G15.7/)
C
      DO 1200 I=1,N
        WRITE (6,6100) (INT(10.0*A(I,K)/VMAX),K=1,M)
 1200 CONTINUE
 6100 FORMAT (1H ,1X,200I1)
C
      RETURN
      END


