fdtd_3d.f90

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

		Ez(0,3,0)=sin(OMIGA*N)
!		Ez(0,-3,0)=-sin(OMIGA*N)
		if(N.le.WL)then 
			Ez(0,3,0)=Ez(0,3,0)/2*(1.0-cos(OMIGA*N)) 
!			Ez(0,-3,0)=Ez(0,-3,0)/2*(1.0-cos(OMIGA*N)) 
		end if

!**********************E Absorbing Boundary*****************
!-------------------------------------6 Planes-----------------------------------------
!---------------Front(y=Jmin,absorb Ez,Ex)----------------------

		do i=Imin+2,Imax-2
			do j=Kmin+2,Kmax-2
				T1=Ez(i,Jmin+1,j)+Ab1z(i,2,j)
				T2=Ab1z(i,3,j)+Ab1z(i,1,j)
				T3=Ab1z(i+1,1,j)+Ab1z(i-1,1,j)+Ab1z(i,1,j+1)+Ab1z(i,1,j-1)+Ab1z(i+1,3,j)+Ab1z(i-1,3,j)+Ab1z(i,3,j+1)+Ab1z(i,3,j-1)
				Ez(i,Jmin,j)=-Ab1z(i,4,j)-T1/3+T2+T3/12
				T1=Ex(i,Jmin+1,j)+Ab1x(i,2,j)
				T2=Ab1x(i,3,j)+Ab1x(i,1,j)
				T3=Ab1x(i+1,1,j)+Ab1x(i-1,1,j)+Ab1x(i,1,j+1)+Ab1x(i,1,j-1)+Ab1x(i+1,3,j)+Ab1x(i-1,3,j)+Ab1x(i,3,j+1)+Ab1x(i,3,j-1)
				Ex(i,Jmin,j)=-Ab1x(i,4,j)-T1/3+T2+T3/12
			end do
		end do
		do i=Imin+1,Imax-1
			Ez(i,Jmin,Kmin+1)=Ab1z(I,3,Kmin+1)-(Ez(I,Jmin+1,Kmin+1)-Ab1z(I,1,Kmin+1))/3
			Ex(i,Jmin,Kmin+1)=Ab1x(I,3,Kmin+1)-(Ex(I,Jmin+1,Kmin+1)-Ab1x(I,1,Kmin+1))/3
			Ez(i,Jmin,Kmax-1)=Ab1z(I,3,Kmax-1)-(Ez(I,Jmin+1,Kmax-1)-Ab1z(I,1,Kmax-1))/3
			Ex(i,Jmin,Kmax-1)=Ab1x(I,3,Kmax-1)-(Ex(I,Jmin+1,Kmax-1)-Ab1x(I,1,Kmax-1))/3
		end do
		do i=Kmin+1,Kmax-1
			Ez(Imin+1,Jmin,i)=Ab1z(Imin+1,3,i)-(Ez(Imin+1,Jmin+1,i)-Ab1z(Imin+1,1,i))/3
			Ex(Imin+1,Jmin,i)=Ab1x(Imin+1,3,i)-(Ex(Imin+1,Jmin+1,i)-Ab1x(Imin+1,1,i))/3
			Ez(Imax-1,Jmin,i)=Ab1z(Imax-1,3,i)-(Ez(Imax-1,Jmin+1,i)-Ab1z(Imax-1,1,i))/3
			Ex(Imax-1,Jmin,i)=Ab1x(Imax-1,3,i)-(Ex(Imax-1,Jmin+1,i)-Ab1x(Imax-1,1,i))/3
		end do
		do i=Imin,Imax
			do j=Kmin,Kmax
				Ab1z(i,2,j)=Ab1z(i,1,j)
				Ab1x(i,2,j)=Ab1x(i,1,j)
				Ab1z(i,4,j)=Ab1z(i,3,j)
				Ab1x(i,4,j)=Ab1x(i,3,j)
				Ab1z(i,1,j)=Ez(i,Jmin,j)     
				Ab1x(i,1,j)=Ex(i,Jmin,j)
				Ab1z(i,3,j)=Ez(i,Jmin+1,j)  
				Ab1x(i,3,j)=Ex(i,Jmin+1,j)
			end do
		end do
!---------------Back(y=Jmax,absorb Ez,Ex)----------------------

		do i=Imin+2,Imax-2
			do j=Kmin+2,Kmax-2
				T1=Ez(i,Jmax-1,j)+Ab2z(i,2,j)
				T2=Ab2z(i,3,j)+Ab2z(i,1,j)
				T3=Ab2z(i+1,1,j)+Ab2z(i-1,1,j)+Ab2z(i,1,j+1)+Ab2z(i,1,j-1)+Ab2z(i+1,3,j)+Ab2z(i-1,3,j)+Ab2z(i,3,j+1)+Ab2z(i,3,j-1)
				Ez(i,Jmax,j)=-Ab2z(i,4,j)-T1/3+T2+T3/12
				T1=Ex(i,Jmax-1,j)+Ab2x(i,2,j)
				T2=Ab2x(i,3,j)+Ab2x(i,1,j)
				T3=Ab2x(i+1,1,j)+Ab2x(i-1,1,j)+Ab2x(i,1,j+1)+Ab2x(i,1,j-1)+Ab2x(i+1,3,j)+Ab2x(i-1,3,j)+Ab2x(i,3,j+1)+Ab2x(i,3,j-1)
				Ex(i,Jmax,j)=-Ab2x(i,4,j)-T1/3+T2+T3/12
			end do
		end do
		do i=Imin+1,Imax-1
			Ez(i,Jmax,Kmin+1)=Ab2z(I,3,Kmin+1)-(Ez(I,Jmax-1,Kmin+1)-Ab2z(I,1,Kmin+1))/3
			Ex(i,Jmax,Kmin+1)=Ab2x(I,3,Kmin+1)-(Ex(I,Jmax-1,Kmin+1)-Ab2x(I,1,Kmin+1))/3
			Ez(i,Jmax,Kmax-1)=Ab2z(I,3,Kmax-1)-(Ez(I,Jmax-1,Kmax-1)-Ab2z(I,1,Kmax-1))/3
			Ex(i,Jmax,Kmax-1)=Ab2x(I,3,Kmax-1)-(Ex(I,Jmax-1,Kmax-1)-Ab2x(I,1,Kmax-1))/3
		end do
		do i=Kmin+1,Kmax-1
			Ez(Imin+1,Jmax,i)=Ab2z(Imin+1,3,i)-(Ez(Imin+1,Jmax-1,i)-Ab2z(Imin+1,1,i))/3
			Ex(Imin+1,Jmax,i)=Ab2x(Imin+1,3,i)-(Ex(Imin+1,Jmax-1,i)-Ab2x(Imin+1,1,i))/3
			Ez(Imax-1,Jmax,i)=Ab2z(Imax-1,3,i)-(Ez(Imax-1,Jmax-1,i)-Ab2z(Imax-1,1,i))/3
			Ex(Imax-1,Jmax,i)=Ab2x(Imax-1,3,i)-(Ex(Imax-1,Jmax-1,i)-Ab2x(Imax-1,1,i))/3
		end do
		do i=Imin,Imax
			do j=Kmin,Kmax
				Ab2z(i,2,j)=Ab2z(i,1,j)
				Ab2x(i,2,j)=Ab2x(i,1,j)
				Ab2z(i,4,j)=Ab2z(i,3,j)
				Ab2x(i,4,j)=Ab2x(i,3,j)
				Ab2z(i,1,j)=Ez(i,Jmax,j)     
				Ab2x(i,1,j)=Ex(i,Jmax,j)
				Ab2z(i,3,j)=Ez(i,Jmax-1,j)  
				Ab2x(i,3,j)=Ex(i,Jmax-1,j)
			end do
		end do

!---------------Down(z=Kmin,absorb Ey,Ex)----------------------

		do i=Imin+2,Imax-2
			do j=Jmin+2,Jmax-2
				T1=Ey(i,j,Kmin+1)+Ab3y(i,j,2)
				T2=Ab3y(i,j,3)+Ab3y(i,j,1)
				T3=Ab3y(i+1,j,1)+Ab3y(i-1,j,1)+Ab3y(i,j+1,1)+Ab3y(i,j-1,1)+Ab3y(i+1,j,3)+Ab3y(i-1,j,3)+Ab3y(i,j+1,3)+Ab3y(i,j-1,3)
				Ey(i,j,Kmin)=-Ab3y(i,j,4)-T1/3+T2+T3/12
				T1=Ex(i,j,Kmin+1)+Ab3x(i,j,2)
				T2=Ab3x(i,j,3)+Ab3x(i,j,1)
				T3=Ab3x(i+1,j,1)+Ab3x(i-1,j,1)+Ab3x(i,j+1,1)+Ab3x(i,j-1,1)+Ab3x(i+1,j,3)+Ab3x(i-1,j,3)+Ab3x(i,j+1,3)+Ab3x(i,j-1,3)
				Ex(i,j,Kmin)=-Ab3x(i,j,4)-T1/3+T2+T3/12
			end do
		end do
		do i=Imin+1,Imax-1
			Ey(i,Jmin+1,Kmin)=Ab3y(I,Jmin+1,3)-(Ey(I,Jmin+1,Kmin+1)-Ab3y(I,Jmin+1,1))/3
			Ex(i,Jmin+1,Kmin)=Ab3x(I,Jmin+1,3)-(Ex(I,Jmin+1,Kmin+1)-Ab3x(I,Jmin+1,1))/3
			Ey(i,Jmax-1,Kmin)=Ab3y(I,Jmax-1,3)-(Ey(I,Jmax-1,Kmin+1)-Ab3y(I,Jmax-1,1))/3
			Ex(i,Jmax-1,Kmin)=Ab3x(I,Jmax-1,3)-(Ex(I,Jmax-1,Kmin+1)-Ab3x(I,Jmax-1,1))/3
		end do
		do i=Jmin+1,Jmax-1
			Ey(Imin+1,i,Kmin)=Ab3y(Imin+1,i,3)-(Ey(Imin+1,i,Kmin+1)-Ab3y(Imin+1,i,1))/3
			Ex(Imin+1,i,Kmin)=Ab3x(Imin+1,i,3)-(Ex(Imin+1,i,Kmin+1)-Ab3x(Imin+1,i,1))/3
			Ey(Imax-1,i,Kmin)=Ab3y(Imax-1,i,3)-(Ey(Imax-1,i,Kmin+1)-Ab3y(Imax-1,i,1))/3
			Ex(Imax-1,i,Kmin)=Ab3x(Imax-1,i,3)-(Ex(Imax-1,i,Kmin+1)-Ab3x(Imax-1,i,1))/3
		end do
		do i=Imin,Imax
			do j=Jmin,Jmax
				Ab3y(i,j,2)=Ab3y(i,j,1)
				Ab3x(i,j,2)=Ab3x(i,j,1)
				Ab3y(i,j,4)=Ab3y(i,j,3)
				Ab3x(i,j,4)=Ab3x(i,j,3)
				Ab3y(i,j,1)=Ey(i,j,Kmin)     
				Ab3x(i,j,1)=Ex(i,j,Kmin)
				Ab3y(i,j,3)=Ey(i,j,Kmin+1)  
				Ab3x(i,j,3)=Ex(i,j,Kmin+1)
			end do
		end do

!---------------Up(z=Kmax,absorb Ey,Ex)----------------------

		do i=Imin+2,Imax-2
			do j=Jmin+2,Jmax-2
				T1=Ey(i,j,Kmax-1)+Ab4y(i,j,2)
				T2=Ab4y(i,j,3)+Ab4y(i,j,1)
				T3=Ab4y(i+1,j,1)+Ab4y(i-1,j,1)+Ab4y(i,j+1,1)+Ab4y(i,j-1,1)+Ab4y(i+1,j,3)+Ab4y(i-1,j,3)+Ab4y(i,j+1,3)+Ab4y(i,j-1,3)
				Ey(i,j,Kmax)=-Ab4y(i,j,4)-T1/3+T2+T3/12
				T1=Ex(i,j,Kmax-1)+Ab4x(i,j,2)
				T2=Ab4x(i,j,3)+Ab4x(i,j,1)
				T3=Ab4x(i+1,j,1)+Ab4x(i-1,j,1)+Ab4x(i,j+1,1)+Ab4x(i,j-1,1)+Ab4x(i+1,j,3)+Ab4x(i-1,j,3)+Ab4x(i,j+1,3)+Ab4x(i,j-1,3)
				Ex(i,j,Kmax)=-Ab4x(i,j,4)-T1/3+T2+T3/12
			end do
		end do
		do i=Imin+1,Imax-1
			Ey(i,Jmin+1,Kmax)=Ab4y(I,Jmin+1,3)-(Ey(I,Jmin+1,Kmax-1)-Ab4y(I,Jmin+1,1))/3
			Ex(i,Jmin+1,Kmax)=Ab4x(I,Jmin+1,3)-(Ex(I,Jmin+1,Kmax-1)-Ab4x(I,Jmin+1,1))/3
			Ey(i,Jmax-1,Kmax)=Ab4y(I,Jmax-1,3)-(Ey(I,Jmax-1,Kmax-1)-Ab4y(I,Jmax-1,1))/3
			Ex(i,Jmax-1,Kmax)=Ab4x(I,Jmax-1,3)-(Ex(I,Jmax-1,Kmax-1)-Ab4x(I,Jmax-1,1))/3
		end do
		do i=Jmin+1,Jmax-1
			Ey(Imin+1,i,Kmax)=Ab4y(Imin+1,i,3)-(Ey(Imin+1,i,Kmax-1)-Ab4y(Imin+1,i,1))/3
			Ex(Imin+1,i,Kmax)=Ab4x(Imin+1,i,3)-(Ex(Imin+1,i,Kmax-1)-Ab4x(Imin+1,i,1))/3
			Ey(Imax-1,i,Kmax)=Ab4y(Imax-1,i,3)-(Ey(Imax-1,i,Kmax-1)-Ab4y(Imax-1,i,1))/3
			Ex(Imax-1,i,Kmax)=Ab4x(Imax-1,i,3)-(Ex(Imax-1,i,Kmax-1)-Ab4x(Imax-1,i,1))/3
		end do
		do i=Imin,Imax
			do j=Jmin,Jmax
				Ab4y(i,j,2)=Ab4y(i,j,1)
				Ab4x(i,j,2)=Ab4x(i,j,1)
				Ab4y(i,j,4)=Ab4y(i,j,3)
				Ab4x(i,j,4)=Ab4x(i,j,3)
				Ab4y(i,j,1)=Ey(i,j,Kmax)     
				Ab4x(i,j,1)=Ex(i,j,Kmax)
				Ab4y(i,j,3)=Ey(i,j,Kmax-1)  
				Ab4x(i,j,3)=Ex(i,j,Kmax-1)
			end do
		end do

!---------------Left(x=Imin,absorb Ez,Ey)----------------------

		do i=Jmin+2,Jmax-2
			do j=Kmin+2,Kmax-2
				T1=Ey(Imin+1,i,j)+Ab5y(2,i,j)
				T2=Ab5y(3,i,j)+Ab5y(1,i,j)
				T3=Ab5y(1,i+1,j)+Ab5y(1,i-1,j)+Ab5y(1,i,j+1)+Ab5y(1,i,j-1)+Ab5y(3,i+1,j)+Ab5y(3,i-1,j)+Ab5y(3,i,j+1)+Ab5y(3,i,j-1)
				Ey(Imin,i,j)=-Ab5y(4,i,j)-T1/3+T2+T3/12
				T1=Ez(Imin+1,i,j)+Ab5z(2,i,j)
				T2=Ab5z(3,i,j)+Ab5z(1,i,j)
				T3=Ab5z(1,i+1,j)+Ab5z(1,i-1,j)+Ab5z(1,i,j+1)+Ab5z(1,i,j-1)+Ab5z(3,i+1,j)+Ab5z(3,i-1,j)+Ab5z(3,i,j+1)+Ab5z(3,i,j-1)
				Ez(Imin,i,j)=-Ab5z(4,i,j)-T1/3+T2+T3/12
			end do
		end do
		do i=Jmin+1,Jmax-1
			Ey(Imin,i,Kmin+1)=Ab5y(3,I,Kmin+1)-(Ey(Imin+1,i,Kmin+1)-Ab5y(1,I,Kmin+1))/3
			Ez(Imin,i,Kmin+1)=Ab5z(3,I,Kmin+1)-(Ez(Imin+1,i,Kmin+1)-Ab5z(1,I,Kmin+1))/3
			Ey(Imin,i,Kmax-1)=Ab5y(3,I,Kmax-1)-(Ey(Imin+1,i,Kmax-1)-Ab5y(1,I,Kmax-1))/3
			Ez(Imin,i,Kmax-1)=Ab5z(3,I,Kmax-1)-(Ez(Imin+1,i,Kmax-1)-Ab5z(1,I,Kmax-1))/3
		end do
		do i=Kmin+1,Kmax-1
			Ey(Imin,Jmin+1,i)=Ab5y(3,Jmin+1,i)-(Ey(Imin+1,Jmin+1,i)-Ab5y(1,Jmin+1,i))/3
			Ez(Imin,Jmin+1,i)=Ab5z(3,Jmin+1,i)-(Ez(Imin+1,Jmin+1,i)-Ab5z(1,Jmin+1,i))/3
			Ey(Imin,Jmax-1,i)=Ab5y(3,Jmax-1,i)-(Ey(Imin+1,Jmax-1,i)-Ab5y(1,Jmax-1,i))/3
			Ez(Imin,Jmax-1,i)=Ab5z(3,Jmax-1,i)-(Ez(Imin+1,Jmax-1,i)-Ab5z(1,Jmax-1,i))/3
		end do
		do i=Jmin,Jmax
			do j=Kmin,Kmax
				Ab5y(2,i,j)=Ab5y(1,i,j)
				Ab5z(2,i,j)=Ab5z(1,i,j)
				Ab5y(4,i,j)=Ab5y(3,i,j)
				Ab5z(4,i,j)=Ab5z(3,i,j)
				Ab5y(1,i,j)=Ey(Imin,i,j)     
				Ab5z(1,i,j)=Ez(Imin,i,j)
				Ab5y(3,i,j)=Ey(Imin+1,i,j)  
				Ab5z(3,i,j)=Ez(Imin+1,i,j)
			end do
		end do

!---------------Right(x=Imax,absorb Ez,Ey)----------------------

		do i=Jmin+2,Jmax-2
			do j=Kmin+2,Kmax-2
				T1=Ey(Imax-1,i,j)+Ab6y(2,i,j)
				T2=Ab6y(3,i,j)+Ab6y(1,i,j)
				T3=Ab6y(1,i+1,j)+Ab6y(1,i-1,j)+Ab6y(1,i,j+1)+Ab6y(1,i,j-1)+Ab6y(3,i+1,j)+Ab6y(3,i-1,j)+Ab6y(3,i,j+1)+Ab6y(3,i,j-1)
				Ey(Imax,i,j)=-Ab6y(4,i,j)-T1/3+T2+T3/12
				T1=Ez(Imax-1,i,j)+Ab6z(2,i,j)
				T2=Ab6z(3,i,j)+Ab6z(1,i,j)
				T3=Ab6z(1,i+1,j)+Ab6z(1,i-1,j)+Ab6z(1,i,j+1)+Ab6z(1,i,j-1)+Ab6z(3,i+1,j)+Ab6z(3,i-1,j)+Ab6z(3,i,j+1)+Ab6z(3,i,j-1)
				Ez(Imax,i,j)=-Ab6z(4,i,j)-T1/3+T2+T3/12
			end do
		end do
		do i=Jmin+1,Jmax-1
			Ey(Imax,i,Kmin+1)=Ab6y(3,I,Kmin+1)-(Ey(Imax-1,i,Kmin+1)-Ab6y(1,I,Kmin+1))/3
			Ez(Imax,i,Kmin+1)=Ab6z(3,I,Kmin+1)-(Ez(Imax-1,i,Kmin+1)-Ab6z(1,I,Kmin+1))/3
			Ey(Imax,i,Kmax-1)=Ab6y(3,I,Kmax-1)-(Ey(Imax-1,i,Kmax-1)-Ab6y(1,I,Kmax-1))/3
			Ez(Imax,i,Kmax-1)=Ab6z(3,I,Kmax-1)-(Ez(Imax-1,i,Kmax-1)-Ab6z(1,I,Kmax-1))/3
		end do
		do i=Kmin+1,Kmax-1
			Ey(Imax,Jmin+1,i)=Ab6y(3,Jmin+1,i)-(Ey(Imax-1,Jmin+1,i)-Ab6y(1,Jmin+1,i))/3
			Ez(Imax,Jmin+1,i)=Ab6z(3,Jmin+1,i)-(Ez(Imax-1,Jmin+1,i)-Ab6z(1,Jmin+1,i))/3
			Ey(Imax,Jmax-1,i)=Ab6y(3,Jmax-1,i)-(Ey(Imax-1,Jmax-1,i)-Ab6y(1,Jmax-1,i))/3
			Ez(Imax,Jmax-1,i)=Ab6z(3,Jmax-1,i)-(Ez(Imax-1,Jmax-1,i)-Ab6z(1,Jmax-1,i))/3
		end do
		do i=Jmin,Jmax
			do j=Kmin,Kmax
				Ab6y(2,i,j)=Ab6y(1,i,j)
				Ab6z(2,i,j)=Ab6z(1,i,j)
				Ab6y(4,i,j)=Ab6y(3,i,j)
				Ab6z(4,i,j)=Ab6z(3,i,j)
				Ab6y(1,i,j)=Ey(Imax,i,j)     
				Ab6z(1,i,j)=Ez(Imax,i,j)
				Ab6y(3,i,j)=Ey(Imax-1,i,j)  
				Ab6z(3,i,j)=Ez(Imax-1,i,j)
			end do
		end do

!*******************H FDTD*********************

		do I=Imin,Imax
			do J=Jmin,Jmax-1
				do K=Kmin,Kmax-1
					T1=Ey(I,J,K+1)-Ey(I,J,K)-Ez(I,J+1,K)+Ez(I,J,K)
					Hx(I,J,K)=FH1(ob(I,J,K),ob(I-1,J,K))*Hx(I,J,K)+FH2(ob(I,J,K),ob(I-1,J,K))*T1
				end do
			end do
		end do

		do I=Imin,Imax-1
			do J=Jmin,Jmax
				do K=Kmin,Kmax-1
					T1=Ez(I+1,J,K)-Ez(I,J,K)-Ex(I,J,K+1)+Ex(I,J,K)
					Hy(I,J,K)=FH1(ob(I,J,K),ob(I,J-1,K))*Hy(I,J,K)+FH2(ob(I,J,K),ob(I,J-1,K))*T1
				end do
			end do
		end do

		do I=Imin,Imax-1
			do J=Jmin,Jmax-1
				do K=Kmin,Kmax
					T1=Ex(I,J+1,K)-Ex(I,J,K)-Ey(I+1,J,K)+Ey(I,J,K)
					Hz(I,J,K)=FH1(ob(I,J,K),ob(I,J,K-1))*Hz(I,J,K)+FH2(ob(I,J,K),ob(I,J,K-1))*T1
				end do
			end do
		end do

!************************Output*************************************

		write(FileName,'(i3)' )N
		Open(1,File='FDTD_3D_'//trim(adjustl(FileName))//'.dat')
		do i=Imin,Imax
			do j=Jmin,Jmax
					write(1,*)Ez(i,j,0)
			end do
		end do
		Close(1)
	end do
end