📄 dadi-fdtd.f90
字号:
! 计算等离子体涂覆导体球的散射截面
PARAMETER(NCon=2,NOut=2,NPML=8,NFFT=16,NT=2**NFFT,NTime=1500)
PARAMETER(ICmax=33,IOmax=ICmax+NCon,IPmax=IOmax+NOut,Imax=IPmax+NPml)
PARAMETER(JCmax=33,JOmax=JCmax+NCon,JPmax=JOmax+NOut,Jmax=JPmax+NPml)
PARAMETER(KCmax=33,KOmax=KCmax+NCon,KPmax=KOmax+NOut,Kmax=KPmax+NPml)
PARAMETER(ICmin=-ICmax,IOmin=-IOmax,IPmin=-IPmax,Imin=-Imax)
PARAMETER(JCmin=-JCmax,JOmin=-JOmax,JPmin=-JPmax,Jmin=-Jmax)
PARAMETER(KCmin=-KCmax,KOmin=-KOmax,KPmin=-KPmax,Kmin=-Kmax)
parameter(NXta=180)
REAL Exy(Imin:Imax-1,Jmin:Jmax,Kmin:Kmax),Exz(Imin:Imax-1,Jmin:Jmax,Kmin:Kmax)
REAL Eyz(Imin:Imax,Jmin:Jmax-1,Kmin:Kmax),Eyx(Imin:Imax,Jmin:Jmax-1,Kmin:Kmax)
REAL Ezx(Imin:Imax,Jmin:Jmax,Kmin:Kmax-1),Ezy(Imin:Imax,Jmin:Jmax,Kmin:Kmax-1)
REAL Hxy(Imin+1:Imax-1,Jmin:Jmax-1,Kmin:Kmax-1),Hxz(Imin+1:Imax-1,Jmin:Jmax-1,Kmin:Kmax-1)
REAL Hyz(Imin:Imax-1,Jmin+1:Jmax-1,Kmin:Kmax-1),Hyx(Imin:Imax-1,Jmin+1:Jmax-1,Kmin:Kmax-1)
REAL Hzx(Imin:Imax-1,Jmin:Jmax-1,Kmin+1:Kmax-1),Hzy(Imin:Imax-1,Jmin:Jmax-1,Kmin+1:Kmax-1)
REAL Pxy(ICmin:ICmax-1,JCmin+1:JCmax-1,KCmin+1:KCmax-1),Pxz(ICmin:ICmax-1,JCmin+1:JCmax-1,KCmin+1:KCmax-1)
REAL Pyz(ICmin+1:ICmax-1,JCmin:JCmax-1,KCmin+1:KCmax-1),Pyx(ICmin+1:ICmax-1,JCmin:JCmax-1,KCmin+1:KCmax-1)
REAL Pzx(ICmin+1:ICmax-1,JCmin+1:JCmax-1,KCmin:KCmax-1),Pzy(ICmin+1:ICmax-1,JCmin+1:JCmax-1,KCmin:KCmax-1)
REAL TmpE(Imin:Imax,Jmin:Jmax,Kmin:Kmax)
REAL Ax(Imax-Imin-1),Bx(Imax-Imin-1),Cx(Imax-Imin-1),Rx(Imax-Imin-1)
REAL Ay(Jmax-Jmin-1),By(Jmax-Jmin-1),Cy(Jmax-Jmin-1),Ry(Jmax-Jmin-1)
REAL Az(Kmax-Kmin-1),Bz(Kmax-Kmin-1),Cz(Kmax-Kmin-1),Rz(Kmax-Kmin-1)
REAL CMedia(0:150,3),Cin(6)
REAL Ei(NT),Es1(NT),Es2(NT),Ercs(0:NXta,0:2*NTime,2)
REAL CJx(2),CJy(2),CJz(2),CMx(2),CMy(2),CMz(2)
REAL Jx1,Jx2,Jy1,Jy2,Jz1,Jz2,Mx1,Mx2,My1,My2,Mz1,Mz2
INTEGER Object(Imin:Imax-1,Jmin:Jmax-1,Kmin:Kmax-1,12)
LOGICAL Condition(3)
EXTERNAL Timef
COMMON /COMMONDSDT/DS,DT
COMMON /Scale/S
COMMON /PlasmaParameter/TEleDen,b,Weff,CPlasma,Tau
!******************** 设定初值 ********************
PI=3.14159265358979323846; CLight=3.e8; Eps0=1.e-9/(36.*Pi)
DS=5.e-3; S=0.5; DT=S*DS/CLight; SS=S/2.
Exy=0.; Exz=0.; Eyz=0.; Eyx=0.; Ezx=0.; Ezy=0.
Hxy=0.; Hxz=0.; Hyz=0.; Hyx=0.; Hzx=0.; Hzy=0.
Pxy=0.; Pxz=0.; Pyz=0.; Pyx=0.; Pzx=0.; Pzy=0.
!******************** 等离子体参数 ****************
TEleDen=1.e16; b=-50.; Weff=1.e10; CPlasma=56.4
!******************** 入射波相关信息 *************
Exta=1.; Ephi=0.
Xtai=0.; Cxi=COSD(Xtai); Sxi=SIND(Xtai)
Phii=0.; Cpi=COSD(Phii); Spi=SIND(Phii)
Source=(FLOAT(ICmin)-0.5)*ABS(Sxi*Cpi)+(FLOAT(JCmin)-0.5)*ABS(Sxi*Spi)+(FLOAT(KCmin)-0.5)*ABS(Cxi)
Cin(1)=Exta*Cxi*Cpi-Ephi*Spi
Cin(2)=Exta*Cxi*Spi+Ephi*Cpi
Cin(3)=-Exta*Sxi
Cin(4)=-Exta*Spi-Ephi*Cxi*Cpi
Cin(5)=Exta*Cpi-Ephi*Cxi*Spi
Cin(6)=Ephi*Sxi
!***** 判定不同Phi角是的输出边界的基准点 **********
Phio=0.; Cpo=COSD(Phio); Spo=SIND(Phio)
IF(Phio.LT.90.) THEN
SX=FLOAT(IOmax); SY=FLOAT(JOmax)
ELSEIF(Phio.LT.180.) THEN
SX=FLOAT(IOmin); SY=FLOAT(JOmax)
ELSEIF(Phio.LT.270.) THEN
SX=FLOAT(IOmin); SY=FLOAT(JOmin)
ELSE
SX=FLOAT(IOmax); SY=FLOAT(JOmin)
ENDIF
!***** 设定整个区域的电磁参数 *********************
CALL MediaCODE(Object,IPmin,IPmax,JPmin,JPmax,KPmin,KPmax,Imin,Imax,Jmin,Jmax,Kmin,Kmax)
CALL DeMediaCode(CMedia,NPml)
NTime2=2*NTime
DO ITime=0,NTime
WRITE(*,100) NTime,ITime,timef(),Exz(0,0,0)+Exy(0,0,0)
100 FORMAT(i4,2x,i5,2x,f13.5,4x,f15.5)
Time1=FLOAT(ITime); Time2=FLOAT(ITime)+0.5; Time3=FLOAT(ITime)+1.
!*********** 提取入射波 ***************************
T1=Time1+Source/S; CALL IncWave(S1,T1); Ei(2*ITime+1)=S1
T2=Time2+Source/S; CALL IncWave(S2,T2); Ei(2*ITime+2)=S2
!*********** 直接完成 Exy 场量迭代 ****************
TmpE(Imin:Imax-1,Jmin:Jmax,Kmin:Kmax)=Exy(Imin:Imax-1,Jmin:Jmax,Kmin:Kmax)
DO I=Imin,Imax-1; Condition(1)=I.GE.ICmin.AND.I.LT.ICmax
DO K=Kmin+1,Kmax-1; Condition(2)=K.GT.KCmin.AND.K.LT.KCmax
DO J=Jmin+1,Jmax-1; Condition(3)=J.GT.JCmin.AND.J.LT.JCmax
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))
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
TExy=Ceh*((Hzx(i,j,k)+Hzy(i,j,k))-(Hzx(i,j-1,k)+Hzy(i,j-1,k)))
Exy(i,j,k)=Cee*Exy(i,j,k)+Cep*Pxy(i,j,k)+TExy
ELSE
Ccc=1.+Sigm*Dt/4.
Cee=1.-Sigm*Dt/4.; Cee=Cee/Ccc
Ceh=SS; Ceh=Ceh/Ccc
TExy=Ceh*((Hzx(i,j,k)+Hzy(i,j,k))-(Hzx(i,j-1,k)+Hzy(i,j-1,k)))
Exy(i,j,k)=Cee*Exy(i,j,k)+TExy
ENDIF
ENDDO
ENDDO
ENDDO
Y1=FLOAT(JCmin)-0.5; Y2=FLOAT(JCmax)+0.5
DO I=ICmin,ICmax-1; X=FLOAT(I)+0.5
DO K=KCmin,KCmax; Z=FLOAT(K)
T1=X*Sxi*Cpi+Y1*Sxi*Spi+Z*Cxi; T1=Time1-(T1-Source)/S
T2=X*Sxi*Cpi+Y2*Sxi*Spi+Z*Cxi; T2=Time1-(T2-Source)/S
CALL IncWave(S1,T1); S1=SS*S1*Cin(6)
CALL IncWave(S2,T2); S2=SS*S2*Cin(6)
Exy(i,jcmin,k)=Exy(i,jcmin,k)-S1
Exy(i,jcmax,k)=Exy(i,jcmax,k)+S2
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*((TmpE(i,j+1,k)+Exz(i,j+1,k))-(TmpE(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=Time1-(T1-Source)/S
T2=X*Sxi*Cpi+Y2*Sxi*Spi+Z*Cxi; T2=Time1-(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(Imin:Imax,Jmin:Jmax-1,Kmin:Kmax)=Eyz(Imin:Imax,Jmin:Jmax-1,Kmin:Kmax)
DO J=Jmin,Jmax-1; Condition(1)=J.GE.JCmin.AND.J.LT.JCmax
DO I=Imin+1,Imax-1; Condition(2)=I.GT.ICmin.AND.I.LT.ICmax
DO K=Kmin+1,Kmax-1; Condition(3)=K.GT.KCmin.AND.K.LT.KCmax
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))
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
TEyz=Ceh*((Hxy(i,j,k)+Hxz(i,j,k))-(Hxy(i,j,k-1)+Hxz(i,j,k-1)))
Eyz(i,j,k)=Cee*Eyz(i,j,k)+Cep*Pyz(i,j,k)+TEyz
ELSE
Ccc=1.+Sigm*Dt/4.
Cee=1.-Sigm*Dt/4.; Cee=Cee/Ccc
Ceh=SS; Ceh=Ceh/Ccc
TEyz=Ceh*((Hxy(i,j,k)+Hxz(i,j,k))-(Hxy(i,j,k-1)+Hxz(i,j,k-1)))
Eyz(i,j,k)=Cee*Eyz(i,j,k)+TEyz
ENDIF
ENDDO
ENDDO
ENDDO
Z1=FLOAT(KCmin)-0.5; Z2=FLOAT(KCmax)+0.5
DO J=JCmin,JCmax-1; Y=FLOAT(J)+0.5
DO I=ICmin,ICmax; X=FLOAT(I)
T1=X*Sxi*Cpi+Y*Sxi*Spi+Z1*Cxi; T1=Time1-(T1-Source)/S
T2=X*Sxi*Cpi+Y*Sxi*Spi+Z2*Cxi; T2=Time1-(T2-Source)/S
CALL IncWave(S1,T1); S1=SS*S1*Cin(4)
CALL IncWave(S2,T2); S2=SS*S2*Cin(4)
Eyz(i,j,kcmin)=Eyz(i,j,kcmin)-S1
Eyz(i,j,kcmax)=Eyz(i,j,kcmax)+S2
ENDDO
ENDDO
!*********** 直接完成 Pyz 场量迭代 ****************
DO J=JCmin,JCmax-1
DO I=ICmin+1,ICmax-1
DO K=KCmin+1,KCmax-1
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))
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
Pyz(i,j,k)=Cpe1*Eyz(i,j,k)+Cpe2*TmpE(i,j,k)-Cpp*Pyz(i,j,k)
ENDDO
ENDDO
ENDDO
!*********** 直接完成 Hxz 场量迭代 ****************
DO I=Imin+1,Imax-1
DO J=Jmin,Jmax-1
DO K=Kmin,Kmax-1
n=Object(i,j,k,8); Sigm=CMedia(n,3)
Ccc=1.+Sigm*Dt/4.
Chh=1.-Sigm*Dt/4.; Chh=Chh/Ccc
Che=SS; Che=Che/Ccc
THxz=Che*((TmpE(i,j,k+1)+Eyx(i,j,k+1))-(TmpE(i,j,k)+Eyx(i,j,k)))
Hxz(i,j,k)=Chh*Hxz(i,j,k)+THxz
ENDDO
ENDDO
ENDDO
Z1=FLOAT(KCmin); Z2=FLOAT(KCmax)
DO I=ICmin,ICmax; X=FLOAT(I)
DO J=JCmin,JCmax-1; Y=FLOAT(J)+0.5
T1=X*Sxi*Cpi+Y*Sxi*Spi+Z1*Cxi; T1=Time1-(T1-Source)/S
T2=X*Sxi*Cpi+Y*Sxi*Spi+Z2*Cxi; T2=Time1-(T2-Source)/S
CALL IncWave(S1,T1); S1=SS*S1*Cin(2)
CALL IncWave(S2,T2); S2=SS*S2*Cin(2)
Hxz(i,j,kcmin-1)=Hxz(i,j,kcmin-1)-S1
Hxz(i,j,kcmax)=Hxz(i,j,kcmax)+S2
ENDDO
ENDDO
!********** 直接完成 Ezx 场量迭代 ******************
TmpE(Imin:Imax,Jmin:Jmax,Kmin:Kmax-1)=Ezx(Imin:Imax,Jmin:Jmax,Kmin:Kmax-1)
DO K=Kmin,Kmax-1; Condition(1)=K.GE.KCmin.AND.K.LT.KCmax
DO J=Jmin+1,Jmax-1; Condition(2)=J.GT.JCmin.AND.J.LT.JCmax
DO I=Imin+1,Imax-1; Condition(3)=I.GT.ICmin.AND.I.LT.ICmax
n1=Object(i-1,j-1,k,5); n2=Object(i,j-1,k,5)
n3=Object(i-1,j,k,5); n4=Object(i,j,k,5)
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
TEzx=Ceh*((Hyz(i,j,k)+Hyx(i,j,k))-(Hyz(i-1,j,k)+Hyx(i-1,j,k)))
Ezx(i,j,k)=Cee*Ezx(i,j,k)+Cep*Pzx(i,j,k)+TEzx
ELSE
Ccc=1.+Sigm*Dt/4.
Cee=1.-Sigm*Dt/4.; Cee=Cee/Ccc
Ceh=SS; Ceh=Ceh/Ccc
TEzx=Ceh*((Hyz(i,j,k)+Hyx(i,j,k))-(Hyz(i-1,j,k)+Hyx(i-1,j,k)))
Ezx(i,j,k)=Cee*Ezx(i,j,k)+TEzx
ENDIF
ENDDO
ENDDO
ENDDO
X1=FLOAT(ICmin)-0.5; X2=FLOAT(ICmax)+0.5
DO K=KCmin,KCmax-1; Z=FLOAT(K)+0.5
DO J=JCmin,JCmax; Y=FLOAT(J)
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(5)
CALL IncWave(S2,T2); S2=SS*S2*Cin(5)
Ezx(icmin,j,k)=Ezx(icmin,j,k)-S1
Ezx(icmax,j,k)=Ezx(icmax,j,k)+S2
ENDDO
ENDDO
!*********** 直接完成 Pzx 场量迭代 ****************
DO K=KCmin,KCmax-1
DO J=JCmin+1,JCmax-1
DO I=ICmin+1,ICmax-1
n1=Object(i-1,j-1,k,5); n2=Object(i,j-1,k,5)
n3=Object(i-1,j,k,5); n4=Object(i,j,k,5)
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
Pzx(i,j,k)=Cpe1*Ezx(i,j,k)+Cpe2*TmpE(i,j,k)-Cpp*Pzx(i,j,k)
ENDDO
ENDDO
ENDDO
!*********** 直接完成 Hyx 场量迭代 ****************
DO J=Jmin+1,Jmax-1
DO K=Kmin,Kmax-1
DO I=Imin,Imax-1
n=Object(i,j,k,10); Sigm=CMedia(n,3)
Ccc=1.+Sigm*Dt/4.
Chh=1.-Sigm*Dt/4.; Chh=Chh/Ccc
Che=SS; Che=Che/Ccc
THyx=Che*((TmpE(i+1,j,k)+Ezy(i+1,j,k))-(TmpE(i,j,k)+Ezy(i,j,k)))
Hyx(i,j,k)=Chh*Hyx(i,j,k)+THyx
ENDDO
ENDDO
ENDDO
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -