📄 dadi-fdtd.f90
字号:
X1=FLOAT(ICmin); X2=FLOAT(ICmax)
DO J=JCmin,JCmax; Y=FLOAT(J)
DO K=KCmin,KCmax-1; Z=FLOAT(K)+0.5
T1=X1*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi; T1=Time1-(T1-Source)/S
T2=X2*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi; T2=Time1-(T2-Source)/S
CALL IncWave(S1,T1); S1=SS*S1*Cin(3)
CALL IncWave(S2,T2); S2=SS*S2*Cin(3)
Hyx(icmin-1,j,k)=Hyx(icmin-1,j,k)-S1
Hyx(icmax,j,k)=Hyx(icmax,j,k)+S2
ENDDO
ENDDO
!****** 通过求解线性方程组来完成Exz的场量迭代 *****
TmpE(ICmin:ICmax-1,JCmin+1:JCmax-1,KCmin+1:KCmax-1)=Exz(ICmin:ICmax-1,JCmin+1:JCmax-1,KCmin+1:KCmax-1)
DO I=Imin,Imax-1
X=FLOAT(I)+0.5; Condition(1)=I.GE.ICmin.AND.I.LT.ICmax
DO J=Jmin+1,Jmax-1
Y=FLOAT(J); Condition(2)=J.GT.JCmin.AND.J.LT.JCmax
DO K=Kmin+1,Kmax-1; Condition(3)=K.GT.KCmin.AND.K.LT.KCmax
KK=K-Kmin
n1=Object(i,j-1,k-1,2); n2=Object(i,j-1,k,2)
n3=Object(i,j,k-1,2); n4=Object(i,j,k,2)
Perm1=0.25*(CMedia(n1,1)+CMedia(n2,1)+CMedia(n3,1)+CMedia(n4,1))
Perm2=0.25*(CMedia(n1,2)+CMedia(n2,2)+CMedia(n3,2)+CMedia(n4,2))
Sigme=0.25*(CMedia(n1,3)+CMedia(n2,3)+CMedia(n3,3)+CMedia(n4,3))
IF(Condition(1).AND.Condition(2).AND.Condition(3)) THEN
Ccc=1.+Sigme*Dt/4.+(Dt*(Perm1-1.)+4.*Tau*(Perm2-1.))/(Dt+4.*Tau)
Cee=1.-Sigme*Dt/4.-(Dt*(Perm1-1.)-4.*Tau*(Perm2-1.))/(Dt+4.*Tau); Cee=Cee/Ccc
Cep=2.*Dt/(Eps0*(Dt+4.*Tau)); Cep=Cep/Ccc
Ceh=SS; Ceh=Ceh/Ccc
Che=SS
Az(KK)=-Ceh*Che; Cz(KK)=-Ceh*Che
Bz(KK)=1.+Ceh*S
THy1=(Hyz(i,j,k)-Che*(Exy(i,j,k+1)-Exy(i,j,k)))+Hyx(i,j,k)
THy2=(Hyz(i,j,k-1)-Che*(Exy(i,j,k)-Exy(i,j,k-1)))+Hyx(i,j,k-1)
Rz(KK)=Cee*Exz(i,j,k)+Cep*Pxz(i,j,k)-Ceh*(THy1-THy2)
ELSE
Ccc=1.+Sigme*Dt/4.
Cee=1.-Sigme*Dt/4.; Cee=Cee/Ccc
Ceh=SS; Ceh=Ceh/Ccc
n=Object(i,j,k,9); Sigmh=CMedia(n,3)
Ccc=1.+Sigmh*Dt/4.
Chh1=1.-Sigmh*Dt/4.; Chh1=Chh1/Ccc
Che1=SS; Che1=Che1/Ccc
n=Object(i,j,k-1,9); Sigmh=CMedia(n,3)
Ccc=1.+Sigmh*Dt/4.
Chh2=1.-Sigmh*Dt/4.; Chh2=Chh2/Ccc
Che2=SS; Che2=Che2/Ccc
IF(K.NE.Kmin+1) Az(KK)=-Ceh*Che2
IF(K.NE.Kmax-1) Cz(KK)=-Ceh*Che1
Bz(KK)=1.+Ceh*(Che1+Che2)
THy1=(Chh1*Hyz(i,j,k)-Che1*(Exy(i,j,k+1)-Exy(i,j,k)))+Hyx(i,j,k)
THy2=(Chh2*Hyz(i,j,k-1)-Che2*(Exy(i,j,k)-Exy(i,j,k-1)))+Hyx(i,j,k-1)
Rz(KK)=Cee*Exz(i,j,k)-Ceh*(THy1-THy2)
IF(Condition(1).AND.Condition(2)) THEN
SELECTCASE(K)
CASE(KCmin-1)
Z=FLOAT(K+1); T1=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T1=Time2-(T1-Source)/S; CALL IncWave(S1,T1); S1=S1*Cin(1)
Rz(KK)=Rz(KK)-SS*SS*S1
CASE(KCmin)
Z=FLOAT(K); T1=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T1=Time2-(T1-Source)/S; CALL IncWave(S1,T1); S1=S1*Cin(1)
Z=FLOAT(K)-0.5; T2=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T2=Time2-(T2-Source)/S; CALL IncWave(S2,T2); S2=S2*Cin(5)
Rz(KK)=Rz(KK)+SS*(SS*S1+S2)
CASE(KCmax)
Z=FLOAT(K); T1=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T1=Time2-(T1-Source)/S; CALL IncWave(S1,T1); S1=S1*Cin(1)
Z=FLOAT(K)+0.5; T2=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T2=Time2-(T2-Source)/S; CALL IncWave(S2,T2); S2=S2*Cin(5)
Rz(KK)=Rz(KK)+SS*(SS*S1-S2)
CASE(KCmax+1)
Z=FLOAT(K-1); T1=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T1=Time2-(T1-Source)/S; CALL IncWave(S1,T1); S1=S1*Cin(1)
Rz(KK)=Rz(KK)-SS*SS*S1
ENDSELECT
ENDIF
ENDIF
ENDDO
!CALL MyLSLTR(Kmax-Kmin-1,Az,Bz,Cz,Rz)
Exz(i,j,kmin+1:kmax-1)=Rz(:)
ENDDO
ENDDO
!*********** 直接完成 Pxz 场量迭代 ****************
DO I=ICmin,ICmax-1
DO J=JCmin+1,JCmax-1
DO K=KCmin+1,KCmax-1
n1=Object(i,j-1,k-1,2); n2=Object(i,j-1,k,2)
n3=Object(i,j,k-1,2); n4=Object(i,j,k,2)
Perm1=0.25*(CMedia(n1,1)+CMedia(n2,1)+CMedia(n3,1)+CMedia(n4,1))
Perm2=0.25*(CMedia(n1,2)+CMedia(n2,2)+CMedia(n3,2)+CMedia(n4,2))
Ccc=Dt+4.*Tau
Cpe1=Eps0*(Dt*(Perm1-1.)+4.*Tau*(Perm2-1.)); Cpe1=Cpe1/Ccc
Cpe2=Eps0*(Dt*(Perm1-1.)-4.*Tau*(Perm2-1.)); Cpe2=Cpe2/Ccc
Cpp=Dt-4.*Tau; Cpp=Cpp/Ccc
Pxz(i,j,k)=Cpe1*Exz(i,j,k)+Cpe2*TmpE(i,j,k)-Cpp*Pxz(i,j,k)
ENDDO
ENDDO
ENDDO
!*********** 直接完成 Hyz 场量迭代 ****************
DO J=Jmin+1,Jmax-1
DO I=Imin,Imax-1
DO K=Kmin,Kmax-1
n=Object(i,j,k,9); Sigm=CMedia(n,3)
Ccc=1.+Sigm*Dt/4.
Chh=1.-Sigm*Dt/4.; Chh=Chh/Ccc
Che=SS; Che=Che/Ccc
THyz=Che*((Exy(i,j,k+1)+Exz(i,j,k+1))-(Exy(i,j,k)+Exz(i,j,k)))
Hyz(i,j,k)=Chh*Hyz(i,j,k)-THyz
ENDDO
ENDDO
ENDDO
Z1=FLOAT(KCmin); Z2=FLOAT(KCmax)
DO J=JCmin,JCmax; Y=FLOAT(J)
DO I=ICmin,ICmax-1; X=FLOAT(I)+0.5
T1=X*Sxi*Cpi+Y*Sxi*Spi+Z1*Cxi; T1=Time2-(T1-Source)/S
T2=X*Sxi*Cpi+Y*Sxi*Spi+Z2*Cxi; T2=Time2-(T2-Source)/S
CALL IncWave(S1,T1); S1=SS*S1*Cin(1)
CALL IncWave(S2,T2); S2=SS*S2*Cin(1)
Hyz(i,j,kcmin-1)=Hyz(i,j,kcmin-1)+S1
Hyz(i,j,kcmax)=Hyz(i,j,kcmax)-S2
ENDDO
ENDDO
!****** 通过求解线性方程组来完成Eyx的场量迭代 *****
TmpE(ICmin+1:ICmax-1,JCmin:JCmax-1,KCmin+1:KCmax-1)=Eyx(ICmin+1:ICmax-1,JCmin:JCmax-1,KCmin+1:KCmax-1)
DO J=Jmin,Jmax-1
Y=FLOAT(J)+0.5; Condition(1)=J.GE.JCmin.AND.J.LT.JCmax
DO K=Kmin+1,Kmax-1
Z=FLOAT(K); Condition(2)=K.GT.KCmin.AND.K.LT.KCmax
DO I=Imin+1,Imax-1; Condition(3)=I.GT.ICmin.AND.I.LT.ICmax
II=I-Imin
n1=Object(i-1,j,k-1,4); n2=Object(i,j,k-1,4)
n3=Object(i-1,j,k,4); n4=Object(i,j,k,4)
Perm1=0.25*(CMedia(n1,1)+CMedia(n2,1)+CMedia(n3,1)+CMedia(n4,1))
Perm2=0.25*(CMedia(n1,2)+CMedia(n2,2)+CMedia(n3,2)+CMedia(n4,2))
Sigme=0.25*(CMedia(n1,3)+CMedia(n2,3)+CMedia(n3,3)+CMedia(n4,3))
IF(Condition(1).AND.Condition(2).AND.Condition(3)) THEN
Ccc=1.+Sigme*Dt/4.+(Dt*(Perm1-1.)+4.*Tau*(Perm2-1.))/(Dt+4.*Tau)
Cee=1.-Sigme*Dt/4.-(Dt*(Perm1-1.)-4.*Tau*(Perm2-1.))/(Dt+4.*Tau); Cee=Cee/Ccc
Cep=2.*Dt/(Eps0*(Dt+4.*Tau)); Cep=Cep/Ccc
Ceh=SS; Ceh=Ceh/Ccc
Che=SS
Ax(II)=-Ceh*Che; Cx(II)=-Ceh*Che
Bx(II)=1.+Ceh*S
THz1=(Hzx(i,j,k)-Che*(Eyz(i+1,j,k)-Eyz(i,j,k)))+Hzy(i,j,k)
THz2=(Hzx(i-1,j,k)-Che*(Eyz(i,j,k)-Eyz(i-1,j,k)))+Hzy(i-1,j,k)
Rx(II)=Cee*Eyx(i,j,k)+Cep*Pyx(i,j,k)-Ceh*(THz1-THz2)
ELSE
Ccc=1.+Sigme*Dt/4.
Cee=1.-Sigme*Dt/4.; Cee=Cee/Ccc
Ceh=SS; Ceh=Ceh/Ccc
n=Object(i,j,k,11); Sigmh=CMedia(n,3)
Ccc=1.+Sigmh*Dt/4.
Chh1=1.-Sigmh*Dt/4.; Chh1=Chh1/Ccc
Che1=SS; Che1=Che1/Ccc
n=Object(i-1,j,k,11); Sigmh=CMedia(n,3)
Ccc=1.+Sigmh*Dt/4.
Chh2=1.-Sigmh*Dt/4.; Chh2=Chh2/Ccc
Che2=SS; Che2=Che2/Ccc
IF(I.NE.Imin+1) Ax(II)=-Ceh*Che2
IF(I.NE.Imax-1) Cx(II)=-Ceh*Che1
Bx(II)=1.+Ceh*(Che1+Che2)
THz1=(Chh1*Hzx(i,j,k)-Che1*(Eyz(i+1,j,k)-Eyz(i,j,k)))+Hzy(i,j,k)
THz2=(Chh2*Hzx(i-1,j,k)-Che2*(Eyz(i,j,k)-Eyz(i-1,j,k)))+Hzy(i-1,j,k)
Rx(II)=Cee*Eyx(i,j,k)-Ceh*(THz1-THz2)
IF(Condition(1).AND.Condition(2)) THEN
SELECTCASE(I)
CASE(ICmin-1)
X=FLOAT(I+1); T1=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T1=Time2-(T1-Source)/S; CALL IncWave(S1,T1); S1=S1*Cin(2)
Rx(II)=Rx(II)-SS*SS*S1
CASE(ICmin)
X=FLOAT(I); T1=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T1=Time2-(T1-Source)/S; CALL IncWave(S1,T1); S1=S1*Cin(2)
X=FLOAT(I)-0.5; T2=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T2=Time2-(T2-Source)/S; CALL IncWave(S2,T2); S2=S2*Cin(6)
Rx(II)=Rx(II)+SS*(SS*S1+S2)
CASE(ICmax)
X=FLOAT(I); T1=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T1=Time2-(T1-Source)/S; CALL IncWave(S1,T1); S1=S1*Cin(2)
X=FLOAT(I)+0.5; T2=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T2=Time2-(T2-Source)/S; CALL IncWave(S2,T2); S2=S2*Cin(6)
Rx(II)=Rx(II)+SS*(SS*S1-S2)
CASE(ICmax+1)
X=FLOAT(I-1); T1=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T1=Time2-(T1-Source)/S; CALL IncWave(S1,T1); S1=S1*Cin(2)
Rx(II)=Rx(II)-SS*SS*S1
ENDSELECT
ENDIF
ENDIF
ENDDO
!CALL MyLSLTR(Imax-Imin-1,Ax,Bx,Cx,Rx)
Eyx(imin+1:imax-1,j,k)=Rx(:)
ENDDO
ENDDO
!*********** 直接完成 Pyx 场量迭代 ****************
DO J=JCmin,JCmax-1
DO K=KCmin+1,KCmax-1
DO I=ICmin+1,ICmax-1
n1=Object(i-1,j,k-1,4); n2=Object(i,j,k-1,4)
n3=Object(i-1,j,k,4); n4=Object(i,j,k,4)
Perm1=0.25*(CMedia(n1,1)+CMedia(n2,1)+CMedia(n3,1)+CMedia(n4,1))
Perm2=0.25*(CMedia(n1,2)+CMedia(n2,2)+CMedia(n3,2)+CMedia(n4,2))
Ccc=Dt+4.*Tau
Cpe1=Eps0*(Dt*(Perm1-1.)+4.*Tau*(Perm2-1.)); Cpe1=Cpe1/Ccc
Cpe2=Eps0*(Dt*(Perm1-1.)-4.*Tau*(Perm2-1.)); Cpe2=Cpe2/Ccc
Cpp=Dt-4.*Tau; Cpp=Cpp/Ccc
Pyx(i,j,k)=Cpe1*Eyx(i,j,k)+Cpe2*TmpE(i,j,k)-Cpp*Pyx(i,j,k)
ENDDO
ENDDO
ENDDO
!*********** 直接完成 Hzx 场量迭代 ****************
DO K=Kmin+1,Kmax-1
DO J=Jmin,Jmax-1
DO I=Imin,Imax-1
n=Object(i,j,k,11); Sigm=CMedia(n,3)
Ccc=1.+Sigm*Dt/4.
Chh=1.-Sigm*Dt/4.; Chh=Chh/Ccc
Che=SS; Che=Che/Ccc
THzx=Che*((Eyz(i+1,j,k)+Eyx(i+1,j,k))-(Eyz(i,j,k)+Eyx(i,j,k)))
Hzx(i,j,k)=Chh*Hzx(i,j,k)-THzx
ENDDO
ENDDO
ENDDO
X1=FLOAT(ICmin); X2=FLOAT(ICmax)
DO K=KCmin,KCmax; Z=FLOAT(K)
DO J=JCmin,JCmax-1; Y=FLOAT(J)+0.5
T1=X1*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi; T1=Time2-(T1-Source)/S
T2=X2*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi; T2=Time2-(T2-Source)/S
CALL IncWave(S1,T1); S1=SS*S1*Cin(2)
CALL IncWave(S2,T2); S2=SS*S2*Cin(2)
Hzx(icmin-1,j,k)=Hzx(icmin-1,j,k)+S1
Hzx(icmax,j,k)=Hzx(icmax,j,k)-S2
ENDDO
ENDDO
!****** 通过求解线性方程组来完成Ezy的场量迭代 *****
TmpE(ICmin+1:ICmax-1,JCmin+1:JCmax-1,KCmin:KCmax-1)=Ezy(ICmin+1:ICmax-1,JCmin+1:JCmax-1,KCmin:KCmax-1)
DO K=Kmin,Kmax-1
Z=FLOAT(K)+0.5; Condition(1)=K.GE.KCmin.AND.K.LT.KCmax
DO I=Imin+1,Imax-1
X=FLOAT(I); Condition(2)=I.GT.ICmin.AND.I.LT.ICmax
DO J=Jmin+1,Jmax-1; Condition(3)=J.GT.JCmin.AND.J.LT.JCmax
JJ=J-Jmin
n1=Object(i-1,j-1,k,6); n2=Object(i-1,j,k,6)
n3=Object(i,j-1,k,6); n4=Object(i,j,k,6)
Perm1=0.25*(CMedia(n1,1)+CMedia(n2,1)+CMedia(n3,1)+CMedia(n4,1))
Perm2=0.25*(CMedia(n1,2)+CMedia(n2,2)+CMedia(n3,2)+CMedia(n4,2))
Sigme=0.25*(CMedia(n1,3)+CMedia(n2,3)+CMedia(n3,3)+CMedia(n4,3))
IF(Condition(1).AND.Condition(2).AND.Condition(3)) THEN
Ccc=1.+Sigme*Dt/4.+(Dt*(Perm1-1.)+4.*Tau*(Perm2-1.))/(Dt+4.*Tau)
Cee=1.-Sigme*Dt/4.-(Dt*(Perm1-1.)-4.*Tau*(Perm2-1.))/(Dt+4.*Tau); Cee=Cee/Ccc
Cep=2.*Dt/(Eps0*(Dt+4.*Tau)); Cep=Cep/Ccc
Ceh=SS; Ceh=Ceh/Ccc
Che=SS
Ay(JJ)=-Ceh*Che; Cy(JJ)=-Ceh*Che
By(JJ)=1.+Ceh*S
THx1=(Hxy(i,j,k)-Che*(Ezx(i,j+1,k)-Ezx(i,j,k)))+Hxz(i,j,k)
THx2=(Hxy(i,j-1,k)-Che*(Ezx(i,j,k)-Ezx(i,j-1,k)))+Hxz(i,j-1,k)
Ry(JJ)=Cee*Ezy(i,j,k)+Cep*Pzy(i,j,k)-Ceh*(THx1-THx2)
ELSE
Ccc=1.+Sigme*Dt/4.
Cee=1.-Sigme*Dt/4.; Cee=Cee/Ccc
Ceh=SS; Ceh=Ceh/Ccc
n=Object(i,j,k,7); Sigmh=CMedia(n,3)
Ccc=1.+Sigmh*Dt/4.
Chh1=1.-Sigmh*Dt/4.; Chh1=Chh1/Ccc
Che1=SS; Che1=Che1/Ccc
n=Object(i,j-1,k,7); Sigmh=CMedia(n,3)
Ccc=1.+Sigmh*Dt/4.
Chh2=1.-Sigmh*Dt/4.; Chh2=Chh2/Ccc
Che2=SS; Che2=Che2/Ccc
IF(J.NE.Jmin+1) Ay(JJ)=-Ceh*Che2
IF(J.NE.Jmax-1) Cy(JJ)=-Ceh*Che1
By(JJ)=1.+Ceh*(Che1+Che2)
THx1=(Chh1*Hxy(i,j,k)-Che1*(Ezx(i,j+1,k)-Ezx(i,j,k)))+Hxz(i,j,k)
THx2=(Chh2*Hxy(i,j-1,k)-Che2*(Ezx(i,j,k)-Ezx(i,j-1,k)))+Hxz(i,j-1,k)
Ry(JJ)=Cee*Ezy(i,j,k)-Ceh*(THx1-THx2)
IF(Condition(1).AND.Condition(2)) THEN
SELECTCASE(J)
CASE(JCmin-1)
Y=FLOAT(J+1); T1=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T1=Time2-(T1-Source)/S; CALL IncWave(S1,T1); S1=S1*Cin(3)
Ry(JJ)=Ry(JJ)-SS*SS*S1
CASE(JCmin)
Y=FLOAT(J); T1=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T1=Time2-(T1-Source)/S; CALL IncWave(S1,T1); S1=S1*Cin(3)
Y=FLOAT(J)-0.5; T2=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T2=Time2-(T2-Source)/S; CALL IncWave(S2,T2); S2=S2*Cin(4)
Ry(JJ)=Ry(JJ)+SS*(SS*S1+S2)
CASE(JCmax)
Y=FLOAT(J); T1=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T1=Time2-(T1-Source)/S; CALL IncWave(S1,T1); S1=S1*Cin(3)
Y=FLOAT(J)+0.5; T2=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T2=Time2-(T2-Source)/S; CALL IncWave(S2,T2); S2=S2*Cin(4)
Ry(JJ)=Ry(JJ)+SS*(SS*S1-S2)
CASE(JCmax+1)
Y=FLOAT(J-1); T1=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T1=Time2-(T1-Source)/S; CALL IncWave(S1,T1); S1=S1*Cin(3)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -