fdtd_3d.f90

一个基于葛德彪书后程序拓展的3D-FDTD的Fortran程序代码

	Program main

	real Ez(-30:30,-30:30,-30:30)
	real Ex(-30:30,-30:30,-30:30)
	real Ey(-30:30,-30:30,-30:30)
	real Hx(-30:30,-30:30,-30:30)
	real Hy(-30:30,-30:30,-30:30)
	real Hz(-30:30,-30:30,-30:30)

!	real Ein(-800:800)
!	real Hin(-800:800)
!	integer IncidentStart,IncidentEnd,Isource
!	real EBin(4)

!	real ECB(-30:30,-30:30,6)
!	real HCB(-30:30,-30:30,6)

	real Ab1x(-30:30,1:4,-30:30)
	real Ab1z(-30:30,1:4,-30:30)     !front
	real Ab2x(-30:30,1:4,-30:30)
	real Ab2z(-30:30,1:4,-30:30)     !back
	real Ab3x(-30:30,-30:30,1:4)
	real Ab3y(-30:30,-30:30,1:4)     !up
	real Ab4x(-30:30,-30:30,1:4)
	real Ab4y(-30:30,-30:30,1:4)     !down
	real Ab5z(1:4,-30:30,-30:30)
	real Ab5y(1:4,-30:30,-30:30)     !left
	real Ab6z(1:4,-30:30,-30:30)
	real Ab6y(1:4,-30:30,-30:30)     !right

	integer ob(-30:30,-30:30,-30:30)

	real FE1(0:10,0:10,0:10,0:10)
	real FE2(0:10,0:10,0:10,0:10)
	real FH1(0:10,0:10)
	real FH2(0:10,0:10)
!	integer TEMflag

!	real FE
!	real FH

	real MU0
	real EPS0
!	real SMU0
!	real SEPS0

	real PI
	real OMIGA
	integer WL
!	real Phi,Sita
	integer TimeStop

	real MaxX
	real MinX
	real MaxY
	real MinY
	real MaxZ
	real MinZ
	real Step

	integer XNum
	integer YNum
	integer ZNum

	integer Imax
	integer Imin
	integer Jmax
	integer Jmin
	integer Kmax
	integer Kmin

	integer Itmax
	integer Itmin
	integer Jtmax
	integer Jtmin
	integer Ktmax
	integer Ktmin

	integer Iomax
	integer Iomin
	integer Jomax
	integer Jomin
	integer Komax
	integer Komin

	integer Ismax
	integer Ismin
	integer Jsmax
	integer Jsmin
	integer Ksmax
	integer Ksmin

	character * 30 FileName
!	character * 1 Blank
	integer FileNameLength
!	logical FileExist

	character * 1 OUTflag
!	character * 1 Respond
!	character * 2 WaveMode
!	character * 30 rem

	real k1,Z,WaveLength

	real Temp,T1,T2,T3,T4,T5,T6,EI0,ER0
	complex CT1,CT2,CT3,CT4
	integer I,J,K,L,N,II,III,IFlag,NN

	real Media(1:4,0:10)
	integer MediaNumber

!*************************************Initializing**************************************************************

	WaveLength=1
	WL=5

!	TEMFlag=1			!1TM,-1TE

	Ismin=-10
	Ismax=10
	Jsmin=-10
	Jsmax=10
	Ksmin=-10
	Ksmax=10

	Itmin=-15
	Itmax=15
	Jtmin=-15
	Jtmax=15
	Ktmin=-15
	Ktmax=15

!	T1=Itmax*Itmax+Jtmax*Jtmax+Ktmax*Ktmax
!	Isource=-sqrt(T1)-3		

	Imin=-25
	Imax=25
	Jmin=-25
	Jmax=25
	Kmin=-25
	Kmax=25

	Iomin=-20
	Iomax=20
	Jomin=-20
	Jomax=20
	Komin=-20
	Komax=20

!	Phi=0.0						!degree
!	Phi=Phi/180.0*Pi			!rad
!	Sita=0.0
!	Sita=Sita/180.0*Pi	

	TimeStop=200

!	write(*,*) char(7)
!	write(*,'(11X,23h Now initializing)')

	MediaNumber=1

	Media(1,1)=1.0
	Media(2,1)=1.0
	Media(3,1)=0.0
	Media(4,1)=0.0

	Media(1,0)=1.0
	Media(2,0)=1.0
	Media(3,0)=0.0
	Media(4,0)=0.0

	do I=Imin,Imax
		do J=Jmin,Jmax
			do K=Kmin,Kmax
				ob(I,J,K)=0
			end do
		end do
	end do

	do I=Ismin,Ismax
		do J=Jsmin,Jsmax
			do K=Ksmin,Ksmax
				if((I*I+J*J+K*K).le.30)then
					ob(I,J,K)=1
				end if
			end do
		end do
	end do

!	IncidentStart=-710
!	IncidentEnd=710

	Pi=3.14159265
	MU0=Pi*0.0000004
	EPS0=8.85e-12
	Z=sqrt(MU0/EPS0)
	k1=2*Pi/WaveLength

	do I=1,MediaNumber
		Media(3,I)=-Media(3,I)*Z/k1
		Media(4,I)=-Media(4,I)/Z/k1
	end do

!	if(TEMFlag.eq.-1)then				
!		Temp=EPS0
!		EPS0=MU0
!		MU0=Temp
!		Z=sqrt(MU0/EPS0)
!	end if

!	SMU0=sqrt(MU0)
!	SEPS0=sqrt(EPS0)
	OMIGA=Pi/WL
!	cp=cos(Phi)
!	sp=sin(Phi)

!	FE=TEMFlag*0.5
!	FH=TEMFlag*0.5

!	if(TEMFlag.eq.1)then
	do I=0,MediaNumber
		do J=0,MediaNumber
			do K=0,MediaNumber
				do L=0,MediaNumber
					FE1(I,J,K,L)=((Media(1,I)+Media(1,J)+Media(1,K)+Media(1,L))/2+(Media(3,I)+Media(3,J)+Media(3,K)+Media(3,L))/4*Pi/WL)/((Media(1,I)+Media(1,J)+Media(1,K)+Media(1,L))/2-(Media(3,I)+Media(3,J)+Media(3,K)+Media(3,L))/4*Pi/WL)
					FE2(I,J,K,L)=1.0/((Media(1,I)+Media(1,J)+Media(1,K)+Media(1,L))/2-(Media(3,I)+Media(3,J)+Media(3,K)+Media(3,L))/4*Pi/WL)
				end do
			end do
			FH1(I,J)=((Media(2,I)+Media(2,J))+(Media(4,I)+Media(4,J))/2*Pi/WL)/((Media(2,I)+Media(2,J))-(Media(4,I)+Media(4,J))/2*Pi/WL)
			FH2(I,J)=1.0/((Media(2,I)+Media(2,J))-(Media(4,I)+Media(4,J))/2*Pi/WL)
		end do
	end do
!	else							!???????????????????????????
!	do I=0,MediaNumber
!		do J=0,MediaNumber
!			FH1(I,J)=((Media(1,I)+Media(1,J))+(Media(3,I)+Media(3,J))/2*Pi/WL)/((Media(1,I)+Media(1,J))-(Media(3,I)+Media(3,J))/2*Pi/WL)
!			FH2(I,J)=-1.0/((Media(1,I)+Media(1,J))-(Media(3,I)+Media(3,J))/2*Pi/WL)
!			FE1(I,J)=((Media(2,I)+Media(2,J))+(Media(4,I)+Media(4,J))/2*Pi/WL)/((Media(2,I)+Media(2,J))-(Media(4,I)+Media(4,J))/2*Pi/WL)
!			FE2(I,J)=-1.0/((Media(2,I)+Media(2,J))-(Media(4,I)+Media(4,J))/2*Pi/WL)
!			end do
!		end do
!	end if

!************************************Main Loop***************************************************************

	OUTflag='I'
	N=0
999	Continue
	do while (N.LT.TimeStop)
		N=N+1
		write(*,1)N
1		format(1H+,10X,'Time step',I5,' is in processing......')


!********************E FDTD*****************************

		do I=Imin+1,Imax-1
			do J=Jmin+1,Jmax-1
				do K=Kmin,Kmax-1
					T1=Hy(I,J,K)-Hy(I-1,J,K)-Hx(I,J,K)+Hx(I,J-1,K)
					Ez(I,J,K)=FE1(ob(I,J,K),ob(I-1,J,K),ob(I,J-1,K),ob(I-1,J-1,K))*Ez(I,J,K)+FE2(ob(I,J,K),ob(I-1,J,K),ob(I,J-1,K),ob(I-1,J-1,K))*T1
				end do
			end do
		end do

		do I=Imin+1,Imax-1
			do J=Jmin,Jmax-1
				do K=Kmin+1,Kmax-1
					T1=Hx(I,J,K)-Hx(I,J,K-1)-Hz(I,J,K)+Hz(I-1,J,K)
					Ey(I,J,K)=FE1(ob(I,J,K),ob(I-1,J,K),ob(I,J,K-1),ob(I-1,J,K-1))*Ey(I,J,K)+FE2(ob(I,J,K),ob(I-1,J,K),ob(I,J,K-1),ob(I-1,J,K-1))*T1
				end do
			end do
		end do

		do I=Imin,Imax-1
			do J=Jmin+1,Jmax-1
				do K=Kmin+1,Kmax-1
					T1=Hz(I,J,K)-Hz(I,J-1,K)-Hy(I,J,K)+Hy(I,J,K-1)
					Ex(I,J,K)=FE1(ob(I,J,K),ob(I,J,K-1),ob(I,J-1,K),ob(I,J-1,K-1))*Ex(I,J,K)+FE2(ob(I,J,K),ob(I,J,K-1),ob(I,J-1,K),ob(I,J-1,K-1))*T1
				end do
			end do
		end do

!************************Incidentwave**************************