📄 dadi-fdtd.f90
字号:
Ry(JJ)=Ry(JJ)-SS*SS*S1
ENDSELECT
ENDIF
ENDIF
ENDDO
!CALL MyLSLTR(Jmax-Jmin-1,Ay,By,Cy,Ry)
Ezy(i,jmin+1:jmax-1,k)=Ry(:)
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)+Ezy(i,j+1,k))-(Ezx(i,j,k)+Ezy(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
write(99,*) time1,exy(0,0,0)+exz(0,0,0),exy(0,0,komax)+exz(0,0,komax)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
ITime2=2*ITime
DO IXta=0,NXta
Xtao=FLOAT(IXta); Cxo=COSD(Xtao); Sxo=SIND(Xtao)
SZ=FLOAT(KOmin); IF(IXta.LT.90) SZ=FLOAT(KOmax)
CJx(1)=Cxo*Cpo; CJx(2)=-Spo; CMx(1)=CJx(2); CMx(2)=-CJx(1)
CJy(1)=Cxo*Spo; CJy(2)=Cpo; CMy(1)=CJy(2); CMy(2)=-CJy(1)
CJz(1)=-Sxo; CJz(2)=0.; CMz(1)=CJz(2); CMz(2)=-CJz(1)
DJ1=FLOAT(JOmin); DJ2=FLOAT(JOmax)
DO I=IOmin,IOmax-1; DI=FLOAT(I)+0.5
DO K=KOmin,KOmax; DK=FLOAT(K)
Jx1=-0.5*(Hzx(I,JOmin,K)+Hzy(I,JOmin,K)+Hzx(I,JOmin-1,K)+Hzy(I,JOmin-1,K))
Jx2=0.5*(Hzx(I,JOmax,K)+Hzy(I,JOmax,K)+Hzx(I,JOmax-1,K)+Hzy(I,JOmax-1,K))
Mz1=-Exy(I,JOmin,K)-Exz(I,JOmin,K)
Mz2=Exy(I,JOmax,K)+Exz(I,JOmax,K)
IF(K.EQ.KOmin.OR.K.EQ.KOmax) THEN
Jx1=0.5*Jx1; Jx2=0.5*Jx2; Mz1=0.5*Mz1; Mz2=0.5*Mz2
ENDIF
T1=(SX-DI)*Sxo*Cpo+(SY-DJ1)*Sxo*Spo+(SZ-DK)*Cxo; T1=T1/SS; N1=FLOOR(T1); T1=T1-N1
T2=(SX-DI)*Sxo*Cpo+(SY-DJ2)*Sxo*Spo+(SZ-DK)*Cxo; T2=T2/SS; N2=FLOOR(T2); T2=T2-N2
IF((ITime2+N1+1).LE.NTime2) THEN
Ercs(IXta,ITime2+N1,1)=(1.-T1)*(Jx1*CJx(1)+Mz1*CMz(1))+Ercs(IXta,ITime2+N1,1)
Ercs(IXta,ITime2+N1,2)=(1.-T1)*(Jx1*CJx(2)+Mz1*CMz(2))+Ercs(IXta,ITime2+N1,2)
Ercs(IXta,ITime2+N1+1,1)=T1*(Jx1*CJx(1)+Mz1*CMz(1))+Ercs(IXta,ITime2+N1+1,1)
Ercs(IXta,ITime2+N1+1,2)=T1*(Jx1*CJx(2)+Mz1*CMz(2))+Ercs(IXta,ITime2+N1+1,2)
ENDIF
IF((ITime2+N2+1).LE.NTime2) THEN
Ercs(IXta,ITime2+N2,1)=(1.-T1)*(Jx2*CJx(1)+Mz2*CMz(1))+Ercs(IXta,ITime2+N2,1)
Ercs(IXta,ITime2+N2,2)=(1.-T1)*(Jx2*CJx(2)+Mz2*CMz(2))+Ercs(IXta,ITime2+N2,2)
Ercs(IXta,ITime2+N2+1,1)=T1*(Jx2*CJx(1)+Mz2*CMz(1))+Ercs(IXta,ITime2+N2+1,1)
Ercs(IXta,ITime2+N2+1,2)=T1*(Jx2*CJx(2)+Mz2*CMz(2))+Ercs(IXta,ITime2+N2+1,2)
ENDIF
ENDDO
ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DK1=FLOAT(KOmin); DK2=FLOAT(KOmax)
DO I=IOmin,IOmax-1; DI=FLOAT(I)+0.5
DO J=JOmin,JOmax; DJ=FLOAT(J)
Jx1=0.5*(Hyz(I,J,KOmin)+Hyx(I,J,KOmin)+Hyz(I,J,KOmin-1)+Hyx(I,J,KOmin-1))
Jx2=-0.5*(Hyz(I,J,KOmax)+Hyx(I,J,KOmax)+Hyz(I,J,KOmax-1)+Hyx(I,J,KOmax-1))
My1=Exy(I,J,KOmin)+Exz(I,J,KOmin)
My2=-Exy(I,J,KOmax)-Exz(I,J,KOmax)
IF(J.EQ.JOmin.OR.J.EQ.JOmax) THEN
Jx1=0.5*Jx1; Jx2=0.5*Jx2; My1=0.5*My1; My2=0.5*My2
ENDIF
T1=(SX-DI)*Sxo*Cpo+(SY-DJ)*Sxo*Spo+(SZ-DK1)*Cxo; T1=T1/SS; N1=FLOOR(T1); T1=T1-N1
T2=(SX-DI)*Sxo*Cpo+(SY-DJ)*Sxo*Spo+(SZ-DK2)*Cxo; T2=T2/SS; N2=FLOOR(T2); T2=T2-N2
IF((ITime2+N1+1).LE.NTime2) THEN
Ercs(IXta,ITime2+N1,1)=(1.-T1)*(Jx1*CJx(1)+My1*CMy(1))+Ercs(IXta,ITime2+N1,1)
Ercs(IXta,ITime2+N1,2)=(1.-T1)*(Jx1*CJx(2)+My1*CMy(2))+Ercs(IXta,ITime2+N1,2)
Ercs(IXta,ITime2+N1+1,1)=T1*(Jx1*CJx(1)+My1*CMy(1))+Ercs(IXta,ITime2+N1+1,1)
Ercs(IXta,ITime2+N1+1,2)=T1*(Jx1*CJx(2)+My1*CMy(2))+Ercs(IXta,ITime2+N1+1,2)
ENDIF
IF((ITime2+N2+1).LE.NTime2) THEN
Ercs(IXta,ITime2+N2,1)=(1.-T1)*(Jx2*CJx(1)+My2*CMy(1))+Ercs(IXta,ITime2+N2,1)
Ercs(IXta,ITime2+N2,2)=(1.-T1)*(Jx2*CJx(2)+My2*CMy(2))+Ercs(IXta,ITime2+N2,2)
Ercs(IXta,ITime2+N2+1,1)=T1*(Jx2*CJx(1)+My2*CMy(1))+Ercs(IXta,ITime2+N2+1,1)
Ercs(IXta,ITime2+N2+1,2)=T1*(Jx2*CJx(2)+My2*CMy(2))+Ercs(IXta,ITime2+N2+1,2)
ENDIF
ENDDO
ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DK1=FLOAT(KOmin); DK2=FLOAT(KOmax)
DO J=JOmin,JOmax-1; DJ=FLOAT(J)+0.5
DO I=IOmin,IOmax; DI=FLOAT(I)
Jy1=-0.5*(Hxy(I,J,KOmin)+Hxz(I,J,KOmin)+Hxy(I,J,KOmin-1)+Hxz(I,J,KOmin-1))
Jy2=0.5*(Hxy(I,J,KOmax)+Hxz(I,J,KOmax)+Hxy(I,J,KOmax-1)+Hxz(I,J,KOmax-1))
Mx1=-Eyz(I,J,KOmin)-Eyx(I,J,KOmin)
Mx2=Eyz(I,J,KOmax)+Eyx(I,J,KOmax)
IF(I.EQ.IOmin.OR.I.EQ.IOmax) THEN
Jy1=0.5*Jy1; Jy2=0.5*Jy2; Mx1=0.5*Mx1; Mx2=0.5*Mx2
ENDIF
T1=(SX-DI)*Sxo*Cpo+(SY-DJ)*Sxo*Spo+(SZ-DK1)*Cxo; T1=T1/SS; N1=FLOOR(T1); T1=T1-N1
T2=(SX-DI)*Sxo*Cpo+(SY-DJ)*Sxo*Spo+(SZ-DK2)*Cxo; T2=T2/SS; N2=FLOOR(T2); T2=T2-N2
IF((ITime2+N1+1).LE.NTime2) THEN
Ercs(IXta,ITime2+N1,1)=(1.-T1)*(Jy1*CJy(1)+Mx1*CMx(1))+Ercs(IXta,ITime2+N1,1)
Ercs(IXta,ITime2+N1,2)=(1.-T1)*(Jy1*CJy(2)+Mx1*CMx(2))+Ercs(IXta,ITime2+N1,2)
Ercs(IXta,ITime2+N1+1,1)=T1*(Jy1*CJy(1)+Mx1*CMx(1))+Ercs(IXta,ITime2+N1+1,1)
Ercs(IXta,ITime2+N1+1,2)=T1*(Jy1*CJy(2)+Mx1*CMx(2))+Ercs(IXta,ITime2+N1+1,2)
ENDIF
IF((ITime2+N2+1).LE.NTime2) THEN
Ercs(IXta,ITime2+N2,1)=(1.-T1)*(Jy2*CJy(1)+Mx2*CMx(1))+Ercs(IXta,ITime2+N2,1)
Ercs(IXta,ITime2+N2,2)=(1.-T1)*(Jy2*CJy(2)+Mx2*CMx(2))+Ercs(IXta,ITime2+N2,2)
Ercs(IXta,ITime2+N2+1,1)=T1*(Jy2*CJy(1)+Mx2*CMx(1))+Ercs(IXta,ITime2+N2+1,1)
Ercs(IXta,ITime2+N2+1,2)=T1*(Jy2*CJy(2)+Mx2*CMx(2))+Ercs(IXta,ITime2+N2+1,2)
ENDIF
ENDDO
ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DI1=FLOAT(IOmin); DI2=FLOAT(IOmax)
DO J=JOmin,JOmax-1; DJ=FLOAT(J)+0.5
DO K=KOmin,KOmax; DK=FLOAT(K)
Jy1=0.5*(Hzx(IOmin,J,K)+Hzy(IOmin,J,K)+Hzx(IOmin-1,J,K)+Hzy(IOmin-1,J,K))
Jy2=-0.5*(Hzx(IOmax,J,K)+Hzy(IOmax,J,K)+Hzx(IOmax-1,J,K)+Hzy(IOmax-1,J,K))
Mz1=Eyz(IOmin,J,K)+Eyx(IOmin,J,K)
Mz2=-Eyz(IOmax,J,K)-Eyx(IOmax,J,K)
IF(K.EQ.KOmin.OR.K.EQ.KOmax) THEN
Jy1=0.5*Jy1; Jy2=0.5*Jy2; Mz1=0.5*Mz1; Mz2=0.5*Mz2
ENDIF
T1=(SX-DI1)*Sxo*Cpo+(SY-DJ)*Sxo*Spo+(SZ-DK)*Cxo; T1=T1/SS; N1=FLOOR(T1); T1=T1-N1
T2=(SX-DI2)*Sxo*Cpo+(SY-DJ)*Sxo*Spo+(SZ-DK)*Cxo; T2=T2/SS; N2=FLOOR(T2); T2=T2-N2
IF((ITime2+N1+1).LE.NTime2) THEN
Ercs(IXta,ITime2+N1,1)=(1.-T1)*(Jy1*CJy(1)+Mz1*CMz(1))+Ercs(IXta,ITime2+N1,1)
Ercs(IXta,ITime2+N1,2)=(1.-T1)*(Jy1*CJy(2)+Mz1*CMz(2))+Ercs(IXta,ITime2+N1,2)
Ercs(IXta,ITime2+N1+1,1)=T1*(Jy1*CJy(1)+Mz1*CMz(1))+Ercs(IXta,ITime2+N1+1,1)
Ercs(IXta,ITime2+N1+1,2)=T1*(Jy1*CJy(2)+Mz1*CMz(2))+Ercs(IXta,ITime2+N1+1,2)
ENDIF
IF((ITime2+N2+1).LE.NTime2) THEN
Ercs(IXta,ITime2+N2,1)=(1.-T1)*(Jy2*CJy(1)+Mz2*CMz(1))+Ercs(IXta,ITime2+N2,1)
Ercs(IXta,ITime2+N2,2)=(1.-T1)*(Jy2*CJy(2)+Mz2*CMz(2))+Ercs(IXta,ITime2+N2,2)
Ercs(IXta,ITime2+N2+1,1)=T1*(Jy2*CJy(1)+Mz2*CMz(1))+Ercs(IXta,ITime2+N2+1,1)
Ercs(IXta,ITime2+N2+1,2)=T1*(Jy2*CJy(2)+Mz2*CMz(2))+Ercs(IXta,ITime2+N2+1,2)
ENDIF
ENDDO
ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DI1=FLOAT(IOmin); DI2=FLOAT(IOmax)
DO K=KOmin,KOmax-1; DK=FLOAT(K)+0.5
DO J=JOmin,JOmax; DJ=FLOAT(J)
Jz1=-0.5*(Hyz(IOmin,J,K)+Hyx(IOmin,J,K)+Hyz(IOmin-1,J,K)+Hyx(IOmin-1,J,K))
Jz2=0.5*(Hyz(IOmax,J,K)+Hyx(IOmax,J,K)+Hyz(IOmax-1,J,K)+Hyx(IOmax-1,J,K))
My1=-Ezx(IOmin,J,K)-Ezy(IOmin,J,K)
My2=Ezx(IOmax,J,K)+Ezy(IOmax,J,K)
IF(J.EQ.JOmin.OR.J.EQ.JOmax) THEN
Jz1=0.5*Jz1; Jz2=0.5*Jz2; My1=0.5*My1; My2=0.5*My2
ENDIF
T1=(SX-DI1)*Sxo*Cpo+(SY-DJ)*Sxo*Spo+(SZ-DK)*Cxo; T1=T1/SS; N1=FLOOR(T1); T1=T1-N1
T2=(SX-DI2)*Sxo*Cpo+(SY-DJ)*Sxo*Spo+(SZ-DK)*Cxo; T2=T2/SS; N2=FLOOR(T2); T2=T2-N2
IF((ITime2+N1+1).LE.NTime2) THEN
Ercs(IXta,ITime2+N1,1)=(1.-T1)*(Jz1*CJz(1)+My1*CMy(1))+Ercs(IXta,ITime2+N1,1)
Ercs(IXta,ITime2+N1,2)=(1.-T1)*(Jz1*CJz(2)+My1*CMy(2))+Ercs(IXta,ITime2+N1,2)
Ercs(IXta,ITime2+N1+1,1)=T1*(Jz1*CJz(1)+My1*CMy(1))+Ercs(IXta,ITime2+N1+1,1)
Ercs(IXta,ITime2+N1+1,2)=T1*(Jz1*CJz(2)+My1*CMy(2))+Ercs(IXta,ITime2+N1+1,2)
ENDIF
IF((ITime2+N2+1).LE.NTime2) THEN
Ercs(IXta,ITime2+N2,1)=(1.-T1)*(Jz2*CJz(1)+My2*CMy(1))+Ercs(IXta,ITime2+N2,1)
Ercs(IXta,ITime2+N2,2)=(1.-T1)*(Jz2*CJz(2)+My2*CMy(2))+Ercs(IXta,ITime2+N2,2)
Ercs(IXta,ITime2+N2+1,1)=T1*(Jz2*CJz(1)+My2*CMy(1))+Ercs(IXta,ITime2+N2+1,1)
Ercs(IXta,ITime2+N2+1,2)=T1*(Jz2*CJz(2)+My2*CMy(2))+Ercs(IXta,ITime2+N2+1,2)
ENDIF
ENDDO
ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
DJ1=FLOAT(JOmin); DJ2=FLOAT(JOmax)
DO K=KOmin,KOmax-1; DK=FLOAT(K)+0.5
DO I=IOmin,IOmax; DI=FLOAT(I)
Jz1=0.5*(Hxy(I,JOmin,K)+Hxz(I,JOmin,K)+Hxy(I,JOmin-1,K)+Hxz(I,JOmin-1,K))
Jz2=-0.5*(Hxy(I,JOmax,K)+Hxz(I,JOmax,K)+Hxy(I,JOmax-1,K)+Hxz(I,JOmax-1,K))
Mx1=Ezx(I,JOmin,K)+Ezy(I,JOmin,K)
Mx2=-Ezx(I,JOmax,K)-Ezy(I,JOmax,K)
IF(I.EQ.IOmin.OR.I.EQ.IOmax) THEN
Jz1=0.5*Jz1; Jz2=0.5*Jz2; Mx1=0.5*Mx1; Mx2=0.5*Mx2
ENDIF
T1=(SX-DI)*Sxo*Cpo+(SY-DJ1)*Sxo*Spo+(SZ-DK)*Cxo; T1=T1/SS; N1=FLOOR(T1); T1=T1-N1
T2=(SX-DI)*Sxo*Cpo+(SY-DJ2)*Sxo*Spo+(SZ-DK)*Cxo; T2=T2/SS; N2=FLOOR(T2); T2=T2-N2
IF((ITime2+N1+1).LE.NTime2) THEN
Ercs(IXta,ITime2+N1,1)=(1.-T1)*(Jz1*CJz(1)+Mx1*CMx(1))+Ercs(IXta,ITime2+N1,1)
Ercs(IXta,ITime2+N1,2)=(1.-T1)*(Jz1*CJz(2)+Mx1*CMx(2))+Ercs(IXta,ITime2+N1,2)
Ercs(IXta,ITime2+N1+1,1)=T1*(Jz1*CJz(1)+Mx1*CMx(1))+Ercs(IXta,ITime2+N1+1,1)
Ercs(IXta,ITime2+N1+1,2)=T1*(Jz1*CJz(2)+Mx1*CMx(2))+Ercs(IXta,ITime2+N1+1,2)
ENDIF
IF((ITime2+N2+1).LE.NTime2) THEN
Ercs(IXta,ITime2+N2,1)=(1.-T1)*(Jz2*CJz(1)+Mx2*CMx(1))+Ercs(IXta,ITime2+N2,1)
Ercs(IXta,ITime2+N2,2)=(1.-T1)*(Jz2*CJz(2)+Mx2*CMx(2))+Ercs(IXta,ITime2+N2,2)
Ercs(IXta,ITime2+N2+1,1)=T1*(Jz2*CJz(1)+Mx2*CMx(1))+Ercs(IXta,ITime2+N2+1,1)
Ercs(IXta,ITime2+N2+1,2)=T1*(Jz2*CJz(2)+Mx2*CMx(2))+Ercs(IXta,ITime2+N2+1,2)
ENDIF
ENDDO
ENDDO
ENDDO
!*********** 直接完成 Exz 场量迭代 ****************
TmpE(Imin:Imax-1,Jmin:Jmax,Kmin:Kmax)=Exz(Imin:Imax-1,Jmin:Jmax,Kmin:Kmax)
DO I=Imin,Imax-1; Condition(1)=I.GE.ICmin.AND.I.LT.ICmax
DO J=Jmin+1,Jmax-1; Condition(2)=J.GT.JCmin.AND.J.LT.JCmax
DO K=Kmin+1,Kmax-1; Condition(3)=K.GT.KCmin.AND.K.LT.KCmax
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))
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
TExz=Ceh*((Hyz(i,j,k)+Hyx(i,j,k))-(Hyz(i,j,k-1)+Hyx(i,j,k-1)))
Exz(i,j,k)=Cee*Exz(i,j,k)+Cep*Pxz(i,j,k)-TExz
ELSE
Ccc=1.+Sigm*Dt/4.
Cee=1.-Sigm*Dt/4.; Cee=Cee/Ccc
Ceh=SS; Ceh=Ceh/Ccc
TExz=Ceh*((Hyz(i,j,k)+Hyx(i,j,k))-(Hyz(i,j,k-1)+Hyx(i,j,k-1)))
Exz(i,j,k)=Cee*Exz(i,j,k)-TExz
ENDIF
ENDDO
ENDDO
ENDDO
Z1=FLOAT(KCmin)-0.5; Z2=FLOAT(KCmax)+0.5
DO I=ICmin,ICmax-1; X=FLOAT(I)+0.5
DO J=JCmin,JCmax; Y=FLOAT(J)
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(5)
CALL IncWave(S2,T2); S2=SS*S2*Cin(5)
Exz(i,j,kcmin)=Exz(i,j,kcmin)+S1
Exz(i,j,kcmax)=Exz(i,j,kcmax)-S2
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)+TmpE(i,j,k+1))-(Exy(i,j,k)+TmpE(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
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -