📄 dadi-fdtd.f90
字号:
ENDDO
!*********** 直接完成 Eyx 场量迭代 ****************
TmpE(Imin:Imax,Jmin:Jmax-1,Kmin:Kmax)=Eyx(Imin:Imax,Jmin:Jmax-1,Kmin:Kmax)
DO J=Jmin,Jmax-1; Condition(1)=J.GE.JCmin.AND.J.LT.JCmax
DO K=Kmin+1,Kmax-1; Condition(2)=K.GT.KCmin.AND.K.LT.KCmax
DO I=Imin+1,Imax-1; Condition(3)=I.GT.ICmin.AND.I.LT.ICmax
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))
Sigm=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.+Sigm*Dt/4.+(Dt*(Perm1-1.)+4.*Tau*(Perm2-1.))/(Dt+4.*Tau)
Cee=1.-Sigm*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
TEyx=Ceh*((Hzx(i,j,k)+Hzy(i,j,k))-(Hzx(i-1,j,k)+Hzy(i-1,j,k)))
Eyx(i,j,k)=Cee*Eyx(i,j,k)+Cep*Pyx(i,j,k)-TEyx
ELSE
Ccc=1.+Sigm*Dt/4.
Cee=1.-Sigm*Dt/4.; Cee=Cee/Ccc
Ceh=SS; Ceh=Ceh/Ccc
TEyx=Ceh*((Hzx(i,j,k)+Hzy(i,j,k))-(Hzx(i-1,j,k)+Hzy(i-1,j,k)))
Eyx(i,j,k)=Cee*Eyx(i,j,k)-TEyx
ENDIF
ENDDO
ENDDO
ENDDO
X1=FLOAT(ICmin)-0.5; X2=FLOAT(ICmax)+0.5
DO J=JCmin,JCmax-1; Y=FLOAT(J)+0.5
DO K=KCmin,KCmax; Z=FLOAT(K)
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(6)
CALL IncWave(S2,T2); S2=SS*S2*Cin(6)
Eyx(icmin,j,k)=Eyx(icmin,j,k)+S1
Eyx(icmax,j,k)=Eyx(icmax,j,k)-S2
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)+TmpE(i+1,j,k))-(Eyz(i,j,k)+TmpE(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(Imin:Imax,Jmin:Jmax,Kmin:Kmax-1)=Ezy(Imin:Imax,Jmin:Jmax,Kmin:Kmax-1)
DO K=Kmin,Kmax-1; Condition(1)=K.GE.KCmin.AND.K.LT.KCmax
DO I=Imin+1,Imax-1; Condition(2)=I.GT.ICmin.AND.I.LT.ICmax
DO J=Jmin+1,Jmax-1; Condition(3)=J.GT.JCmin.AND.J.LT.JCmax
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))
Sigm=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.+Sigm*Dt/4.+(Dt*(Perm1-1.)+4.*Tau*(Perm2-1.))/(Dt+4.*Tau)
Cee=1.-Sigm*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
TEzy=Ceh*((Hxy(i,j,k)+Hxz(i,j,k))-(Hxy(i,j-1,k)+Hxz(i,j-1,k)))
Ezy(i,j,k)=Cee*Ezy(i,j,k)+Cep*Pzy(i,j,k)-TEzy
ELSE
Ccc=1.+Sigm*Dt/4.
Cee=1.-Sigm*Dt/4.; Cee=Cee/Ccc
Ceh=SS; Ceh=Ceh/Ccc
TEzy=Ceh*((Hxy(i,j,k)+Hxz(i,j,k))-(Hxy(i,j-1,k)+Hxz(i,j-1,k)))
Ezy(i,j,k)=Cee*Ezy(i,j,k)-TEzy
ENDIF
ENDDO
ENDDO
ENDDO
Y1=FLOAT(JCmin)-0.5; Y2=FLOAT(JCmax)+0.5
DO K=KCmin,KCmax-1; Z=FLOAT(K)+0.5
DO I=ICmin,ICmax; X=FLOAT(I)
T1=X*Sxi*Cpi+Y1*Sxi*Spi+Z*Cxi; T1=Time2-(T1-Source)/S
T2=X*Sxi*Cpi+Y2*Sxi*Spi+Z*Cxi; T2=Time2-(T2-Source)/S
CALL IncWave(S1,T1); S1=SS*S1*Cin(4)
CALL IncWave(S2,T2); S2=SS*S2*Cin(4)
Ezy(i,jcmin,k)=Ezy(i,jcmin,k)+S1
Ezy(i,jcmax,k)=Ezy(i,jcmax,k)-S2
ENDDO
ENDDO
!*********** 直接完成 Pzy 场量迭代 ****************
DO K=KCmin,KCmax-1
DO I=ICmin+1,ICmax-1
DO J=JCmin+1,JCmax-1
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))
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
Pzy(i,j,k)=Cpe1*Ezy(i,j,k)+Cpe2*TmpE(i,j,k)-Cpp*Pzy(i,j,k)
ENDDO
ENDDO
ENDDO
!*********** 直接完成 Hxy 场量迭代 ****************
DO I=Imin+1,Imax-1
DO K=Kmin,Kmax-1
DO J=Jmin,Jmax-1
n=Object(i,j,k,7); Sigm=CMedia(n,3)
Ccc=1.+Sigm*Dt/4.
Chh=1.-Sigm*Dt/4.; Chh=Chh/Ccc
Che=SS; Che=Che/Ccc
THxy=Che*((Ezx(i,j+1,k)+TmpE(i,j+1,k))-(Ezx(i,j,k)+TmpE(i,j,k)))
Hxy(i,j,k)=Chh*Hxy(i,j,k)-THxy
ENDDO
ENDDO
ENDDO
Y1=FLOAT(JCmin); Y2=FLOAT(JCmax)
DO I=ICmin,ICmax; X=FLOAT(I)
DO K=KCmin,KCmax-1; Z=FLOAT(K)+0.5
T1=X*Sxi*Cpi+Y1*Sxi*Spi+Z*Cxi; T1=Time2-(T1-Source)/S
T2=X*Sxi*Cpi+Y2*Sxi*Spi+Z*Cxi; T2=Time2-(T2-Source)/S
CALL IncWave(S1,T1); S1=SS*S1*Cin(3)
CALL IncWave(S2,T2); S2=SS*S2*Cin(3)
Hxy(i,jcmin-1,k)=Hxy(i,jcmin-1,k)+S1
Hxy(i,jcmax,k)=Hxy(i,jcmax,k)-S2
ENDDO
ENDDO
!****** 通过求解线性方程组来完成Exy的场量迭代 *****
TmpE(ICmin:ICmax-1,JCmin+1:JCmax-1,KCmin+1:KCmax-1)=Exy(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 K=Kmin+1,Kmax-1
Z=FLOAT(K); Condition(2)=K.GT.KCmin.AND.K.LT.KCmax
DO J=Jmin+1,Jmax-1; Condition(3)=J.GT.JCmin.AND.J.LT.JCmax
JJ=J-Jmin
n1=Object(i,j-1,k-1,1); n2=Object(i,j,k-1,1)
n3=Object(i,j-1,k,1); n4=Object(i,j,k,1)
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
THz1=(Hzy(i,j,k)+Che*(Exz(i,j+1,k)-Exz(i,j,k)))+Hzx(i,j,k)
THz2=(Hzy(i,j-1,k)+Che*(Exz(i,j,k)-Exz(i,j-1,k)))+Hzx(i,j-1,k)
Ry(JJ)=Cee*Exy(i,j,k)+Cep*Pxy(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,12); 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,12); 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)
THz1=(Chh1*Hzy(i,j,k)+Che1*(Exz(i,j+1,k)-Exz(i,j,k)))+Hzx(i,j,k)
THz2=(Chh2*Hzy(i,j-1,k)+Che2*(Exz(i,j,k)-Exz(i,j-1,k)))+Hzx(i,j-1,k)
Ry(JJ)=Cee*Exy(i,j,k)+Ceh*(THz1-THz2)
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=Time3-(T1-Source)/S; CALL IncWave(S1,T1); S1=S1*Cin(1)
Ry(JJ)=Ry(JJ)-SS*SS*S1
CASE(JCmin)
Y=FLOAT(J); T1=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T1=Time3-(T1-Source)/S; CALL IncWave(S1,T1); S1=S1*Cin(1)
Y=FLOAT(J)-0.5; T2=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T2=Time3-(T2-Source)/S; CALL IncWave(S2,T2); S2=S2*Cin(6)
Ry(JJ)=Ry(JJ)+SS*(SS*S1-S2)
CASE(JCmax)
Y=FLOAT(J); T1=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T1=Time3-(T1-Source)/S; CALL IncWave(S1,T1); S1=S1*Cin(1)
Y=FLOAT(J)+0.5; T2=X*Sxi*Cpi+Y*Sxi*Spi+Z*Cxi
T2=Time3-(T2-Source)/S; CALL IncWave(S2,T2); S2=S2*Cin(6)
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=Time3-(T1-Source)/S; CALL IncWave(S1,T1); S1=S1*Cin(1)
Ry(JJ)=Ry(JJ)-SS*SS*S1
ENDSELECT
ENDIF
ENDIF
ENDDO
! CALL MyLSLTR(Jmax-Jmin-1,Ay,By,Cy,Ry)
Exy(i,jmin+1:jmax-1,k)=Ry(:)
ENDDO
ENDDO
!*********** 直接完成 Pxy 场量迭代 ****************
DO I=ICmin,ICmax-1
DO K=KCmin+1,KCmax-1
DO J=JCmin+1,JCmax-1
n1=Object(i,j-1,k-1,1); n2=Object(i,j,k-1,1)
n3=Object(i,j-1,k,1); n4=Object(i,j,k,1)
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
Pxy(i,j,k)=Cpe1*Exy(i,j,k)+Cpe2*TmpE(i,j,k)-Cpp*Pxy(i,j,k)
ENDDO
ENDDO
ENDDO
!*********** 直接完成 Hzy 场量迭代 ****************
DO K=Kmin+1,Kmax-1
DO I=Imin,Imax-1
DO J=Jmin,Jmax-1
n=Object(i,j,k,12); Sigm=CMedia(n,3)
Ccc=1.+Sigm*Dt/4.
Chh=1.-Sigm*Dt/4.; Chh=Chh/Ccc
Che=SS; Che=Che/Ccc
THzy=Che*((Exy(i,j+1,k)+Exz(i,j+1,k))-(Exy(i,j,k)+Exz(i,j,k)))
Hzy(i,j,k)=Chh*Hzy(i,j,k)+THzy
ENDDO
ENDDO
ENDDO
Y1=FLOAT(JCmin); Y2=FLOAT(JCmax)
DO K=KCmin,KCmax; Z=FLOAT(K)
DO I=ICmin,ICmax-1; X=FLOAT(I)+0.5
T1=X*Sxi*Cpi+Y1*Sxi*Spi+Z*Cxi; T1=Time3-(T1-Source)/S
T2=X*Sxi*Cpi+Y2*Sxi*Spi+Z*Cxi; T2=Time3-(T2-Source)/S
CALL IncWave(S1,T1); S1=SS*S1*Cin(1)
CALL IncWave(S2,T2); S2=SS*S2*Cin(1)
Hzy(i,jcmin-1,k)=Hzy(i,jcmin-1,k)-S1
Hzy(i,jcmax,k)=Hzy(i,jcmax,k)+S2
ENDDO
ENDDO
!****** 通过求解线性方程组来完成Eyz的场量迭代 *****
TmpE(ICmin+1:ICmax-1,JCmin:JCmax-1,KCmin+1:KCmax-1)=Eyz(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 I=Imin+1,Imax-1
X=FLOAT(I); Condition(2)=I.GT.ICmin.AND.I.LT.ICmax
DO K=Kmin+1,Kmax-1; Condition(3)=K.GT.KCmin.AND.K.LT.KCmax
KK=K-Kmin
n1=Object(i-1,j,k-1,3); n2=Object(i-1,j,k,3)
n3=Object(i,j,k-1,3); n4=Object(i,j,k,3)
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
THx1=(Hxz(i,j,k)+Che*(Eyx(i,j,k+1)-Eyx(i,j,k)))+Hxy(i,j,k)
THx2=(Hxz(i,j,k-1)+Che*(Eyx(i,j,k)-Eyx(i,j,k-1)))+Hxy(i,j,k-1)
Rz(KK)=Cee*Eyz(i,j,k)+Cep*Pyz(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,8); Sigmh=CMedia(n,3)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -