ms.f90

是用FDTD对波导进行研究的程序

PROGRAM WG_3D

IMPLICIT NONE
REAL(KIND=8),ALLOCATABLE::HX(:,:,:),HY(:,:,:),HZ(:,:,:),EX(:,:,:),EY(:,:,:),EZ(:,:,:),&
                          HXY(:,:,:),HXZ(:,:,:),HYX(:,:,:),HYZ(:,:,:),HZX(:,:,:),HZY(:,:,:),&
						  EXY(:,:,:),EXZ(:,:,:),EYX(:,:,:),EYZ(:,:,:),EZX(:,:,:),EZY(:,:,:)
REAL(KIND=8),ALLOCATABLE::C1_HXY(:,:,:),C2_HXY(:,:,:),C1_HXZ(:,:,:),C2_HXZ(:,:,:),C1_HYX(:,:,:),C2_HYX(:,:,:),&
                          C1_HYZ(:,:,:),C2_HYZ(:,:,:),C1_HZX(:,:,:),C2_HZX(:,:,:),C1_HZY(:,:,:),C2_HZY(:,:,:),&
						  C1_EXY(:,:,:),C2_EXY(:,:,:),C1_EXZ(:,:,:),C2_EXZ(:,:,:),C1_EYX(:,:,:),C2_EYX(:,:,:),&
						  C1_EYZ(:,:,:),C2_EYZ(:,:,:),C1_EZX(:,:,:),C2_EZX(:,:,:),C1_EZY(:,:,:),C2_EZY(:,:,:)
REAL(KIND=8)::B,A,L
INTEGER::M,N,S,I,J,K,KL,K_PML1,K_PML2,K_SOURCE,NT,N0
REAL(KIND=8)::SIGMA_MAX,SIGMA_X,SIGMA_Y,SIGMA_Z
REAL(KIND=8)::D,DT,E0,U0,V0,DX,DY,DZ,REF,FREQ,MAG_S21,PHA_S21,PI
REAL(KIND=8)::C1,C2
COMPLEX::CJ=CMPLX(0.0,1.0),C_DFT
INTEGER::K_PORT1,K_PORT2,BUFFER_LAYERS
COMPLEX::VOLTAGE1(201),VOLTAGE2(201)
INTEGER::NF,N_FREQ
REAL(KIND=8)::FMIN,FMAX
INTEGER,PARAMETER::PML_LAYERS=20
REAL,EXTERNAL::EXIN
LOGICAL ALIVE
COMMON FMIN,FMAX
PI=3.14159265
E0=1.0E-9/(36.0*PI)           !真空中的介电常数
U0=4.0E-7*PI                     !真空中的磁导率
V0=1.0/SQRT(E0*U0)                    !光速


INQUIRE(FILE="INPUTFILE.DAT",EXIST=ALIVE)
IF(.NOT.ALIVE)THEN
  WRITE(*,*)"THE INPUTFILE DOES NOT EXIST."
END IF
OPEN(10,FILE="INPUTFILE.DAT")
  READ(10,*) DX,DY,DZ                                   !输入的每个小Yee元胞的长度
  READ(10,*) NT                                          !不出意外是时间的次数
  READ(10,*) A,B,L                                      !波导的实际长度和宽度
  READ(10,*) FMIN,FMAX,N_FREQ                             !FREQUENCY IN GHZ
CLOSE(10)
	DX=DX/1000.0;DY=DY/1000.0;DZ=DZ/1000.0
	A=A/1000.0;B=B/1000.0;L=L/1000.0
    FMIN=FMIN*1.0E9;FMAX=FMAX*1.0E9                     !这边仅是为了使单位统一一下

!确定边界
BUFFER_LAYERS=10
IF(B/DX-FLOAT(INT(B/DX))<0.5)THEN
    M=INT(B/DX)+1
  ELSE
    M=INT(B/DX)+1+1
END IF
DX=B/FLOAT(M)

IF(A/DY-FLOAT(INT(A/DY))<0.5)THEN
    N=1+INT(A/DY)
  ELSE
    N=1+INT(A/DY)+1
END IF
DY=A/FLOAT(N-1)

K_PML1=PML_LAYERS+1
K_SOURCE=K_PML1+BUFFER_LAYERS
K_PORT1=K_SOURCE+BUFFER_LAYERS
K_PORT2=K_PORT1+INT(L/DZ)
K_PML2=K_PORT2+BUFFER_LAYERS
DZ=L/FLOAT(K_PORT2-K_PORT1)
S=K_PML2+(K_PML1-1)

D=MIN(DX,DY,DZ)
DT=D/(2.0*V0)                        !
REF=1.0E-10
SIGMA_MAX=-(3.0/2.0*E0*V0*LOG(REF))/FLOAT(PML_LAYERS-1)/D
C1=DT/E0
C2=DT/U0

ALLOCATE(HX(M,N-1,S-1),HY(M-1,N,S-1),HZ(M-1,N-1,S),EX(M-1,N,S),EY(M,N-1,S),EZ(M,N,S-1),&
    HXY(M,N-1,S-1),HXZ(M,N-1,S-1),HYX(M-1,N,S-1),HYZ(M-1,N,S-1),HZX(M-1,N-1,S),HZY(M-1,N-1,S),&
	EXY(M-1,N,S),EXZ(M-1,N,S),EYX(M,N-1,S),EYZ(M,N-1,S),EZX(M,N,S-1),EZY(M,N,S-1))
ALLOCATE(C1_HXY(M,N-1,S-1),C2_HXY(M,N-1,S-1),C1_HXZ(M,N-1,S-1),C2_HXZ(M,N-1,S-1),C1_HYX(M-1,N,S-1),C2_HYX(M-1,N,S-1),&
         C1_HYZ(M-1,N,S-1),C2_HYZ(M-1,N,S-1),C1_HZX(M-1,N-1,S),C2_HZX(M-1,N-1,S),C1_HZY(M-1,N-1,S),C2_HZY(M-1,N-1,S),&
		 C1_EXY(M-1,N,S),C2_EXY(M-1,N,S),C1_EXZ(M-1,N,S),C2_EXZ(M-1,N,S),C1_EYX(M,N-1,S),C2_EYX(M,N-1,S),C1_EYZ(M,N-1,S),C2_EYZ(M,N-1,S),&
		 C1_EZX(M,N,S-1),C2_EZX(M,N,S-1),C1_EZY(M,N,S-1),C2_EZY(M,N,S-1))

! MAGNETIC COEFFICIENT'S CACULATION IN PML
DO K=1,S-1
  DO J=1,N-1
    DO I=2,M-1
	    SIGMA_Y=0.0
	    C1_HXY(I,J,K)=(2.0*E0-SIGMA_Y*DT)/(2.0*E0+SIGMA_Y*DT)
		C2_HXY(I,J,K)=2.0*E0*DT/U0/(2.0*E0+SIGMA_Y*DT)
	END DO
  END DO
END DO

DO K=1,S-1
  DO J=1,N-1
    DO I=2,M-1
	  IF(K<K_PML1) THEN
        SIGMA_Z=SIGMA_MAX*((FLOAT(K_PML1-K)-0.5)/FLOAT(PML_LAYERS))**2
	  ELSE IF(K>=K_PML2) THEN
	    SIGMA_Z=SIGMA_MAX*((FLOAT(K-K_PML2)+0.5)/FLOAT(PML_LAYERS))**2
	  ELSE
	    SIGMA_Z=0.0
	  END IF
	    C1_HXZ(I,J,K)=(2.0*E0-SIGMA_Z*DT)/(2.0*E0+SIGMA_Z*DT)
		C2_HXZ(I,J,K)=2.0*E0*DT/U0/(2.0*E0+SIGMA_Z*DT)
	END DO
  END DO
END DO

DO K=1,S-1
  DO J=2,N-1
    DO I=1,M-1
	    SIGMA_X=0.0
	    C1_HYX(I,J,K)=(2.0*E0-SIGMA_X*DT)/(2.0*E0+SIGMA_X*DT)
		C2_HYX(I,J,K)=2.0*E0*DT/U0/(2.0*E0+SIGMA_X*DT)
	END DO
  END DO
END DO

DO K=1,S-1
  DO J=2,N-1
    DO I=1,M-1
	  IF(K<K_PML1) THEN
        SIGMA_Z=SIGMA_MAX*((FLOAT(K_PML1-K)-0.5)/FLOAT(PML_LAYERS))**2
	  ELSE IF(K>=K_PML2) THEN
	    SIGMA_Z=SIGMA_MAX*((FLOAT(K-K_PML2)+0.5)/FLOAT(PML_LAYERS))**2
	  ELSE
	    SIGMA_Z=0.0
	  END IF
	    C1_HYZ(I,J,K)=(2.0*E0-SIGMA_Z*DT)/(2.0*E0+SIGMA_Z*DT)
		C2_HYZ(I,J,K)=2.0*E0*DT/U0/(2.0*E0+SIGMA_Z*DT)
	END DO
  END DO
END DO

DO K=2,S-1
  DO J=1,N-1
    DO I=1,M-1
	    SIGMA_X=0.0
	    C1_HZX(I,J,K)=(2.0*E0-SIGMA_X*DT)/(2.0*E0+SIGMA_X*DT)
		C2_HZX(I,J,K)=2.0*E0*DT/U0/(2.0*E0+SIGMA_X*DT)
	END DO
  END DO
END DO

DO K=2,S-1
  DO J=1,N-1
    DO I=1,M-1
	    SIGMA_Y=0.0
	    C1_HZY(I,J,K)=(2.0*E0-SIGMA_Y*DT)/(2.0*E0+SIGMA_Y*DT)
		C2_HZY(I,J,K)=2.0*E0*DT/U0/(2.0*E0+SIGMA_Y*DT)
	END DO
  END DO
END DO

!  ELECTRICAL COEFFICIENT'S CACULATION IN PML
DO K=2,S-1
  DO J=2,N-1
    DO I=1,M-1
	    SIGMA_Y=0.0
	    C1_EXY(I,J,K)=(2.0*E0-SIGMA_Y*DT)/(2.0*E0+SIGMA_Y*DT)
		C2_EXY(I,J,K)=2.0*DT/(2.0*E0+SIGMA_Y*DT)
	END DO
  END DO
END DO

DO K=2,S-1
  DO J=2,N-1
    DO I=1,M-1
	  IF(K<=K_PML1) THEN
        SIGMA_Z=SIGMA_MAX*(FLOAT(K_PML1-K)/FLOAT(PML_LAYERS))**2
	  ELSE IF(K>=K_PML2) THEN
	    SIGMA_Z=SIGMA_MAX*(FLOAT(K-K_PML2)/FLOAT(PML_LAYERS))**2
	  ELSE
	    SIGMA_Z=0.0
	  END IF
	    C1_EXZ(I,J,K)=(2.0*E0-SIGMA_Z*DT)/(2.0*E0+SIGMA_Z*DT)
		C2_EXZ(I,J,K)=2.0*DT/(2.0*E0+SIGMA_Z*DT)
	END DO
  END DO
END DO

DO K=2,S-1
  DO J=1,N-1
    DO I=2,M-1
	    SIGMA_X=0.0
	    C1_EYX(I,J,K)=(2.0*E0-SIGMA_X*DT)/(2.0*E0+SIGMA_X*DT)
		C2_EYX(I,J,K)=2.0*DT/(2.0*E0+SIGMA_X*DT)
	END DO
  END DO
END DO

DO K=2,S-1
  DO J=1,N-1
    DO I=2,M-1
	  IF(K<=K_PML1) THEN
        SIGMA_Z=SIGMA_MAX*(FLOAT(K_PML1-K)/FLOAT(PML_LAYERS))**2
	  ELSE IF(K>=K_PML2) THEN
	    SIGMA_Z=SIGMA_MAX*(FLOAT(K-K_PML2)/FLOAT(PML_LAYERS))**2
	  ELSE
	    SIGMA_Z=0.0
	  END IF
	    C1_EYZ(I,J,K)=(2.0*E0-SIGMA_Z*DT)/(2.0*E0+SIGMA_Z*DT)
		C2_EYZ(I,J,K)=2.0*DT/(2.0*E0+SIGMA_Z*DT)
	END DO
  END DO
END DO

DO K=1,S-1
  DO J=2,N-1
    DO I=2,M-1
	    SIGMA_X=0.0
	    C1_EZX(I,J,K)=(2.0*E0-SIGMA_X*DT)/(2.0*E0+SIGMA_X*DT)
		C2_EZX(I,J,K)=2.0*DT/(2.0*E0+SIGMA_X*DT)
	END DO
  END DO
END DO

DO K=1,S-1
  DO J=2,N-1
    DO I=2,M-1
	    SIGMA_Y=0.0
	    C1_EZY(I,J,K)=(2.0*E0-SIGMA_Y*DT)/(2.0*E0+SIGMA_Y*DT)
		C2_EZY(I,J,K)=2.0*DT/(2.0*E0+SIGMA_Y*DT)
	END DO
  END DO
END DO

HX=0.0;HY=0.0;HZ=0.0;EX=0.0;EY=0.0;EZ=0.0
HXY=0.0;HXZ=0.0;HYX=0.0;HYZ=0.0;HZX=0.0;HZY=0.0;EXY=0.0;EXZ=0.0;EYX=0.0;EYZ=0.0;EZX=0.0;EZY=0.0
VOLTAGE1=0.0;VOLTAGE2=0.0

OPEN(1,FILE='PORT_VOLT.DAT')

DO N0=1,NT
   IF(MOD(N0,10)==0) WRITE(*,*) "N0=",N0,EX(1,(1+N)/2,K_PORT1),EX(1,(1+N)/2,K_PORT2)     !不明白为什么

!CACULATE THE MAGNETIC FIELD
!HX FIELD
   DO K=K_PML1,K_PML2-1
      DO J=1,N-1
	     DO I=2,M-1
		    IF(K==K_PML1)THEN
			   EY(I,J,K)=EYX(I,J,K)+EYZ(I,J,K)
			ELSE IF(K==K_PML2-1)THEN
			   EY(I,J,K+1)=EYX(I,J,K+1)+EYZ(I,J,K+1)
			END IF
			HX(I,J,K)=HX(I,J,K)-C2*((EZ(I,J+1,K)-EZ(I,J,K))/DY-(EY(I,J,K+1)-EY(I,J,K))/DZ)
		  END DO
	   END DO
   END DO

!HY FIELD
   DO K=K_PML1,K_PML2-1
      DO J=2,N-1
	     DO I=1,M-1
		    IF(K==K_PML1)THEN
			   EX(I,J,K)=EXY(I,J,K)+EXZ(I,J,K)
			ELSE IF(K==K_PML2-1)THEN
			   EX(I,J,K+1)=EXY(I,J,K+1)+EXZ(I,J,K+1)
			END IF
			HY(I,J,K)=HY(I,J,K)-C2*((EX(I,J,K+1)-EX(I,J,K))/DZ-(EZ(I+1,J,K)-EZ(I,J,K))/DX)
		  END DO
	   END DO
   END DO

!HZ FIELD
   DO K=K_PML1+1,K_PML2-1
      DO J=1,N-1
	     DO I=1,M-1
			HZ(I,J,K)=HZ(I,J,K)-C2*((EY(I+1,J,K)-EY(I,J,K))/DX-(EX(I,J+1,K)-EX(I,J,K))/DY)
		  END DO
	   END DO
	END DO

!CALCULATE THE MAGNETIC FIELDS IN THE PML
!HX IN PML PARALELL TO Z=CONSTANT PLANE
   DO K=1,K_PML1-1
      KL=K+K_PML2-1
      DO J=1,N-1
         DO I=2,M-1
            HXY(I,J,K)=C1_HXY(I,J,K)*HXY(I,J,K)-C2_HXY(I,J,K)*(EZX(I,J+1,K)-EZX(I,J,K)+EZY(I,J+1,K)-EZY(I,J,K))/DY
            HXZ(I,J,K)=C1_HXZ(I,J,K)*HXZ(I,J,K)+C2_HXZ(I,J,K)*(EYX(I,J,K+1)-EYX(I,J,K)+EYZ(I,J,K+1)-EYZ(I,J,K))/DZ
            HXY(I,J,KL)=C1_HXY(I,J,KL)*HXY(I,J,KL)-C2_HXY(I,J,KL)*(EZX(I,J+1,KL)-EZX(I,J,KL)+EZY(I,J+1,KL)-EZY(I,J,KL))/DY
            HXZ(I,J,KL)=C1_HXZ(I,J,KL)*HXZ(I,J,KL)+C2_HXZ(I,J,KL)*(EYX(I,J,KL+1)-EYX(I,J,KL)+EYZ(I,J,KL+1)-EYZ(I,J,KL))/DZ
         END DO
      END DO
   END DO

!HY IN PML PARALELL TO Z=CONSTANT PLANE
   DO K=1,K_PML1-1
      KL=K+K_PML2-1
      DO J=2,N-1
         DO I=1,M-1
		    HYX(I,J,K)=C1_HYX(I,J,K)*HYX(I,J,K)+C2_HYX(I,J,K)*(EZX(I+1,J,K)-EZX(I,J,K)+EZY(I+1,J,K)-EZY(I,J,K))/DX
            HYZ(I,J,K)=C1_HYZ(I,J,K)*HYZ(I,J,K)-C2_HYZ(I,J,K)*(EXY(I,J,K+1)-EXY(I,J,K)+EXZ(I,J,K+1)-EXZ(I,J,K))/DZ
		    HYX(I,J,KL)=C1_HYX(I,J,KL)*HYX(I,J,KL)+C2_HYX(I,J,KL)*(EZX(I+1,J,KL)-EZX(I,J,KL)+EZY(I+1,J,KL)-EZY(I,J,KL))/DX
            HYZ(I,J,KL)=C1_HYZ(I,J,KL)*HYZ(I,J,KL)-C2_HYZ(I,J,KL)*(EXY(I,J,KL+1)-EXY(I,J,KL)+EXZ(I,J,KL+1)-EXZ(I,J,KL))/DZ
         END DO
      END DO
   END DO

!HZ IN PML PARALELL TO Z=CONSTANT PLANE
   DO K=2,K_PML1
      KL=K+K_PML2-2
      DO J=1,N-1
         DO I=1,M-1
            HZX(I,J,K)=C1_HZX(I,J,K)*HZX(I,J,K)-C2_HZX(I,J,K)*(EYX(I+1,J,K)-EYX(I,J,K)+EYZ(I+1,J,K)-EYZ(I,J,K))/DX
            HZY(I,J,K)=C1_HZY(I,J,K)*HZY(I,J,K)+C2_HZY(I,J,K)*(EXY(I,J+1,K)-EXY(I,J,K)+EXZ(I,J+1,K)-EXZ(I,J,K))/DY
            HZX(I,J,KL)=C1_HZX(I,J,KL)*HZX(I,J,KL)-C2_HZX(I,J,KL)*(EYX(I+1,J,KL)-EYX(I,J,KL)+EYZ(I+1,J,KL)-EYZ(I,J,KL))/DX
            HZY(I,J,KL)=C1_HZY(I,J,KL)*HZY(I,J,KL)+C2_HZY(I,J,KL)*(EXY(I,J+1,KL)-EXY(I,J,KL)+EXZ(I,J+1,KL)-EXZ(I,J,KL))/DY
         END DO
      END DO
   END DO

!CALCULATE THE ELECTRIC FIELD 

! EX FIELD
  DO K=K_PML1+1,K_PML2-1
    DO J=2,N-1
	   DO I=1,M-1
		     EX(I,J,K)=EX(I,J,K)+C1*((HZ(I,J,K)-HZ(I,J-1,K))/DY-(HY(I,J,K)-HY(I,J,K-1))/DZ)
          IF(K==K_SOURCE) THEN
             EX(I,J,K)=EX(I,J,K)+EXIN(N0,DT,A,DY,J)
          END IF
	   END DO
	END DO
  END DO

!EY FIELD
  DO K=K_PML1+1,K_PML2-1
    DO J=1,N-1
	   DO I=2,M-1
		     EY(I,J,K)=EY(I,J,K)+C1*((HX(I,J,K)-HX(I,J,K-1))/DZ-(HZ(I,J,K)-HZ(I-1,J,K))/DX)
	   END DO
	END DO
  END DO

!EZ FIELD
  DO K=K_PML1,K_PML2-1
    DO J=2,N-1
	   DO I=2,M-1
		     EZ(I,J,K)=EZ(I,J,K)+C1*((HY(I,J,K)-HY(I-1,J,K))/DX-(HX(I,J,K)-HX(I,J-1,K))/DY)
	   END DO
	END DO
  END DO
 
 ! ELECTRIC FIELDS CALCULATION  IN PML
 !EX IN PML PARALELL TO Z=CONSTANT PLANE
  DO K=2,K_PML1
    KL=K+K_PML2-2
    DO J=2,N-1
	   DO I=1,M-1
	      EXY(I,J,K)=C1_EXY(I,J,K)*EXY(I,J,K)+C2_EXY(I,J,K)*(HZX(I,J,K)-HZX(I,J-1,K)+HZY(I,J,K)-HZY(I,J-1,K))/DY
	      EXY(I,J,KL)=C1_EXY(I,J,KL)*EXY(I,J,KL)+C2_EXY(I,J,KL)*(HZX(I,J,KL)-HZX(I,J-1,KL)+HZY(I,J,KL)-HZY(I,J-1,KL))/DY
		  IF(K==K_PML1)THEN
		    EXZ(I,J,K)=C1_EXZ(I,J,K)*EXZ(I,J,K)-C2_EXZ(I,J,K)*(HY(I,J,K)-HYX(I,J,K-1)-HYZ(I,J,K-1))/DZ
		  ELSE
		    EXZ(I,J,K)=C1_EXZ(I,J,K)*EXZ(I,J,K)-C2_EXZ(I,J,K)*(HYX(I,J,K)-HYX(I,J,K-1)+HYZ(I,J,K)-HYZ(I,J,K-1))/DZ
		  END IF

		  IF(KL==K_PML2)THEN
            EXZ(I,J,KL)=C1_EXZ(I,J,KL)*EXZ(I,J,KL)-C2_EXZ(I,J,KL)*(HYX(I,J,KL)+HYZ(I,J,KL)-HY(I,J,KL-1))/DZ
		  ELSE
		    EXZ(I,J,KL)=C1_EXZ(I,J,KL)*EXZ(I,J,KL)-C2_EXZ(I,J,KL)*(HYX(I,J,KL)-HYX(I,J,KL-1)+HYZ(I,J,KL)-HYZ(I,J,KL-1))/DZ
		  END IF
	   END DO
	 END DO
   END DO

 !EY IN PML PARALELL TO Z=CONSTANT PLANE
  DO K=2,K_PML1
     KL=K+K_PML2-2
     DO J=1,N-1
	   DO I=2,M-1
	        EYX(I,J,K)=C1_EYX(I,J,K)*EYX(I,J,K)-C2_EYX(I,J,K)*(HZX(I,J,K)-HZX(I-1,J,K)+HZY(I,J,K)-HZY(I-1,J,K))/DX
	        EYX(I,J,KL)=C1_EYX(I,J,KL)*EYX(I,J,KL)-C2_EYX(I,J,KL)*(HZX(I,J,KL)-HZX(I-1,J,KL)+HZY(I,J,KL)-HZY(I-1,J,KL))/DX
		    IF(K==K_PML1)THEN
		      EYZ(I,J,K)=C1_EYZ(I,J,K)*EYZ(I,J,K)+C2_EYZ(I,J,K)*(HX(I,J,K)-HXY(I,J,K-1)-HXZ(I,J,K-1))/DZ
		    ELSE
		      EYZ(I,J,K)=C1_EYZ(I,J,K)*EYZ(I,J,K)+C2_EYZ(I,J,K)*(HXY(I,J,K)-HXY(I,J,K-1)+HXZ(I,J,K)-HXZ(I,J,K-1))/DZ
		    END IF

		    IF(KL==K_PML2)THEN
		      EYZ(I,J,KL)=C1_EYZ(I,J,KL)*EYZ(I,J,KL)+C2_EYZ(I,J,KL)*(HXY(I,J,KL)+HXZ(I,J,KL)-HX(I,J,KL-1))/DZ
		    ELSE
		      EYZ(I,J,KL)=C1_EYZ(I,J,KL)*EYZ(I,J,KL)+C2_EYZ(I,J,KL)*(HXY(I,J,KL)-HXY(I,J,KL-1)+HXZ(I,J,KL)-HXZ(I,J,KL-1))/DZ
		    END IF
       END DO
	 END DO
  END DO
 
 !EZ IN PML PARALELL TO Z=CONSTANT PLANE
  DO K=1,K_PML1-1
    KL=K+K_PML2-1
    DO J=2,N-1
	   DO I=2,M-1
            EZX(I,J,K)=C1_EZX(I,J,K)*EZX(I,J,K)+C2_EZX(I,J,K)*(HYX(I,J,K)-HYX(I-1,J,K)+HYZ(I,J,K)-HYZ(I-1,J,K))/DX
		    EZY(I,J,K)=C1_EZY(I,J,K)*EZY(I,J,K)-C2_EZY(I,J,K)*(HXY(I,J,K)-HXY(I,J-1,K)+HXZ(I,J,K)-HXZ(I,J-1,K))/DY
            EZX(I,J,KL)=C1_EZX(I,J,KL)*EZX(I,J,KL)+C2_EZX(I,J,KL)*(HYX(I,J,KL)-HYX(I-1,J,KL)+HYZ(I,J,KL)-HYZ(I-1,J,KL))/DX
		    EZY(I,J,KL)=C1_EZY(I,J,KL)*EZY(I,J,KL)-C2_EZY(I,J,KL)*(HXY(I,J,KL)-HXY(I,J-1,KL)+HXZ(I,J,KL)-HXZ(I,J-1,KL))/DY
	   END DO
	END DO
  END DO

  DO NF=1,N_FREQ
    FREQ=FMIN+FLOAT(NF-1)*(FMAX-FMIN)/FLOAT(N_FREQ-1)
    C_DFT=2.0/FLOAT(NT)*EXP(-CJ*2.0*PI*FREQ*FLOAT(N0-1)*DT)
    DO I=1,M-1		   
       VOLTAGE1(NF)=VOLTAGE1(NF)+C_DFT*EX(I,(1+N)/2,K_PORT1)
       VOLTAGE2(NF)=VOLTAGE2(NF)+C_DFT*EX(I,(1+N)/2,K_PORT2)
    END DO
  END DO
  WRITE(1,*) N0,EX(1,(1+N)/2,K_PORT1),EX(1,(1+N)/2,K_PORT2)
END DO                    !END OF TIME ITREATION.

CLOSE(1)

OPEN(UNIT=21,FILE="S2P.DAT",STATUS="UNKNOWN")

DO NF=1,N_FREQ
   FREQ=FMIN+FLOAT(NF-1)*(FMAX-FMIN)/FLOAT(N_FREQ-1)
   MAG_S21=CABS(VOLTAGE2(NF)/VOLTAGE1(NF))
   PHA_S21=ATAN2(AIMAG(VOLTAGE2(NF)/VOLTAGE1(NF)),REAL(VOLTAGE2(NF)/VOLTAGE1(NF)))
   WRITE(21,*) FREQ/1.0E9, MAG_S21,PHA_S21*180.0/PI
END DO
CLOSE(21)
STOP
END                      ! END OF ALL

!激励函数
FUNCTION EXIN(N0,DT,A,DY,J)
IMPLICIT NONE
COMMON FMIN,FMAX
REAL(KIND=8)::DT,DY,A,T_WIDTH,T_CENTER,FMIN,FMAX,PI
INTEGER::N0,J
REAL(KIND=8)::EXIN
PI=3.14159265
T_WIDTH=1.0/(PI*FMAX)
T_CENTER=2.0*T_WIDTH
EXIN=SIN(PI*(J-1)*DY/A)*SIN(PI*(FMIN+FMAX)*(N0*DT))*EXP(-(FLOAT(N0-1)*DT-T_CENTER)**2/T_WIDTH**2)
RETURN
END