subroutines2.f90

三维FDTD,matlal编程,mur边界条件,平面波光源

module subroutines2
	use constants
	use matdata
	implicit none
	contains
	!==========================================================================
	!                        initialization
	!==========================================================================
	subroutine initialization()
		allocate(	e_p(0:mat_num),u_p(0:mat_num),g_p(0:mat_num),&
					c_p(0:mat_num),d_p(0:mat_num)	)
		allocate(	CEx(1:nx,1:ny,1:nz),DEx(1:nx,1:ny,1:nz),&
					CEy(1:nx,1:ny,1:nz),DEy(1:nx,1:ny,1:nz),&
					CEz(1:nx,1:ny,1:nz),DEz(1:nx,1:ny,1:nz),& 
					Ex(1:nx-1,1:ny,1:nz),Ey(1:nx,1:ny-1,1:nz),Ez(1:nx,1:ny,1:nz-1),&
					Exy_old(1:nx-1,1:4,1:nz),Exz_old(1:nx-1,1:ny,1:4),&
					Eyx_old(1:4,1:ny-1,1:nz),Eyz_old(1:nx,1:ny-1,1:4),&
					Ezx_old(1:4,1:ny,1:nz-1),Ezy_old(1:nx,1:4,1:nz-1),&
					Hx(1:nx,1:ny-1,1:nz-1),Hy(1:nx-1,1:ny,1:nz-1),Hz(1:nx-1,1:ny-1,1:nz),&	
					outdata(1:nx-1,1:ny-1,1:nz-1)	)
		e_p=0
		u_p=0
		g_p=0
		
		c_p=0
		d_p=0
		
		CEx=0;DEx=0;CEy=0;DEy=0;CEz=0;DEz=0;

		Ex=0
		Ey=0
		Ez=0
		Hx=0
		Hy=0
		Hz=0
		
		Exy_old=0; Exz_old=0
		Eyx_old=0; Eyz_old=0
		Ezx_old=0; Ezy_old=0

		outdata=0

	end subroutine initialization
	!==========================================================================
	!                       freespace
	!==========================================================================	
	subroutine freespace()
		deallocate(	e_p,u_p,g_p,&
					c_p,d_p	)
		deallocate(	CEx,DEx,CEy,DEy,CEz,DEz,& 
					Ex,Ey,Ez,&
					Hx,Hy,Hz,&	
					outdata	)
	end subroutine freespace
	!==========================================================================
	!          const_c                    const_d
	!==========================================================================
	complex(kind) function const_c(ee,uu,gg)
		complex(kind) :: ee,uu,gg
		const_c=1d0/(sqrt(u0/e0)*(gg/2d0+ee/dt))
		return
	end function const_c

	complex(kind) function const_d(ee,uu,gg)
		complex(kind) :: ee,uu,gg
		const_d=sqrt(u0/e0)*(gg/2d0-ee/dt)
		return
	end function const_d
	!==========================================================================
	!          movon
	!==========================================================================
	subroutine movon()

		!FDTD计算程序
		Ex(:,2:ny-1,2:nz-1)=CEx(1:nx-1,2:ny-1,2:nz-1)*(-DEx(1:nx-1,2:ny-1,2:nz-1)*Ex(:,2:ny-1,2:nz-1) &
			+(Hz(:,2:ny-1,2:nz-1)-Hz(:,1:ny-2,2:nz-1))/dy-(Hy(:,2:ny-1,2:nz-1)-Hy(:,2:ny-1,1:nz-2))/dz)
			!   1D Mur         
		Ex(:,1,:)=Exy_old(:,3,:)+(c*dt-dy)/(c*dt+dy)*(Ex(:,2,:)-Exy_old(:,1,:))
		Ex(:,ny,:)=Exy_old(:,4,:)+(c*dt-dy)/(c*dt+dy)*(Ex(:,ny-1,:)-Exy_old(:,2,:))
		Ex(:,:,1)=Exz_old(:,:,3)+(c*dt-dz)/(c*dt+dz)*(Ex(:,:,2)-Exz_old(:,:,1))
		Ex(:,:,nz)=Exz_old(:,:,4)+(c*dt-dz)/(c*dt+dz)*(Ex(:,:,nz-1)-Exz_old(:,:,2))

		!-----------------------------------------------------------------------------------------------
		Ey(2:nx-1,:,2:nz-1)=CEy(2:nx-1,1:ny-1,2:nz-1)*(-DEy(2:nx-1,1:ny-1,2:nz-1)*Ey(2:nx-1,:,2:nz-1) &
			+(Hx(2:nx-1,:,2:nz-1)-Hx(2:nx-1,:,1:nz-2))/dz-(Hz(2:nx-1,:,2:nz-1)-Hz(1:nx-2,:,2:nz-1))/dx)
		!   1D Mur             
		Ey(1,:,:)=Eyx_old(3,:,:)+(c*dt-dx)/(c*dt+dx)*(Ey(2,:,:)-Eyx_old(1,:,:))
		Ey(nx,:,:)=Eyx_old(4,:,:)+(c*dt-dx)/(c*dt+dx)*(Ey(nx-1,:,:)-Eyx_old(2,:,:))
		Ey(:,:,1)=Eyz_old(:,:,3)+(c*dt-dz)/(c*dt+dz)*(Ey(:,:,2)-Eyz_old(:,:,1))
		Ey(:,:,nz)=Eyz_old(:,:,4)+(c*dt-dz)/(c*dt+dz)*(Ey(:,:,nz-1)-Eyz_old(:,:,2))
	
		!-----------------------------------------------------------------------------------------------
		Ez(2:nx-1,2:ny-1,:)=CEz(2:nx-1,2:ny-1,1:nz-1)*(-DEz(2:nx-1,2:ny-1,1:nz-1)*Ez(2:nx-1,2:ny-1,:) &  
			+(Hy(2:nx-1,2:ny-1,:)-Hy(1:nx-2,2:ny-1,:))/dx-(Hx(2:nx-1,2:ny-1,:)-Hx(2:nx-1,1:ny-2,:))/dy)
		!  1D Mur             
		Ez(:,1,:)=Ezy_old(:,3,:)+(c*dt-dy)/(c*dt+dy)*(Ez(:,2,:)-Ezy_old(:,1,:))
		Ez(:,ny,:)=Ezy_old(:,4,:)+(c*dt-dy)/(c*dt+dy)*(Ez(:,ny-1,:)-Ezy_old(:,2,:))
		Ez(1,:,:)=Ezx_old(3,:,:)+(c*dt-dx)/(c*dt+dx)*(Ez(2,:,:)-Ezx_old(1,:,:))
		Ez(nx,:,:)=Ezx_old(4,:,:)+(c*dt-dx)/(c*dt+dx)*(Ez(nx-1,:,:)-Ezx_old(2,:,:))
		
		!---save electric fileds of the fomer step--------------------------------------------------------------------
		Exy_old(1:nx-1,1,1:nz)=Ex(1:nx-1,1,1:nz); Exy_old(1:nx-1,2,1:nz)=Ex(1:nx-1,ny,1:nz)
		Exy_old(1:nx-1,3,1:nz)=Ex(1:nx-1,2,1:nz); Exy_old(1:nx-1,4,1:nz)=Ex(1:nx-1,ny-1,1:nz)

		Exz_old(1:nx-1,1:ny,1)=Ex(1:nx-1,1:ny,1); Exz_old(1:nx-1,1:ny,2)=Ex(1:nx-1,1:ny,nz)
		Exz_old(1:nx-1,1:ny,3)=Ex(1:nx-1,1:ny,2); Exz_old(1:nx-1,1:ny,4)=Ex(1:nx-1,1:ny,nz-1)

		
		Eyx_old(1,1:ny-1,1:nz)=Ey(1,1:ny-1,1:nz); Eyx_old(2,1:ny-1,1:nz)=Ey(nx,1:ny-1,1:nz)
		Eyx_old(3,1:ny-1,1:nz)=Ey(2,1:ny-1,1:nz); Eyx_old(4,1:ny-1,1:nz)=Ey(nx-1,1:ny-1,1:nz)

		Eyz_old(1:nx,1:ny-1,1)=Ey(1:nx,1:ny-1,1); Eyz_old(1:nx,1:ny-1,2)=Ey(1:nx,1:ny-1,nz)
		Eyz_old(1:nx,1:ny-1,3)=Ey(1:nx,1:ny-1,2); Eyz_old(1:nx,1:ny-1,4)=Ey(1:nx,1:ny-1,nz-1)
		

		Ezx_old(1,1:ny,1:nz-1)=Ez(1,1:ny,1:nz-1); Ezx_old(2,1:ny,1:nz-1)=Ez(nx,1:ny,1:nz-1)
		Ezx_old(3,1:ny,1:nz-1)=Ez(2,1:ny,1:nz-1); Ezx_old(4,1:ny,1:nz-1)=Ez(nx-1,1:ny,1:nz-1)

		Ezy_old(1:nx,1,1:nz-1)=Ez(1:nx,1,1:nz-1); Ezy_old(1:nx,2,1:nz-1)=Ez(1:nx,ny,1:nz-1)
		Ezy_old(1:nx,3,1:nz-1)=Ez(1:nx,2,1:nz-1); Ezy_old(1:nx,4,1:nz-1)=Ez(1:nx,ny-1,1:nz-1)

		!------------------------------------------------------------------------------------------------
		Hx(:,1:ny-1,1:nz-1)=Hx(:,1:ny-1,1:nz-1) &
			+(Ey(:,1:ny-1,2:nz)-Ey(:,1:ny-1,1:nz-1))*ds/2/dz-(Ez(:,2:ny,1:nz-1)-Ez(:,1:ny-1,1:nz-1))*ds/2/dy;              
		!   STWBC               
    
		!------------------------------------------------------------------------------------------------
		Hy(1:nx-1,:,1:nz-1)=Hy(1:nx-1,:,1:nz-1) &
			+(Ez(2:nx,:,1:nz-1)-Ez(1:nx-1,:,1:nz-1))*ds/2/dx-(Ex(1:nx-1,:,2:nz)-Ex(1:nx-1,:,1:nz-1))*ds/2/dz;


		!------------------------------------------------------------------------------------------------
		Hz(1:nx-1,1:ny-1,:)=Hz(1:nx-1,1:ny-1,:) &
			+(Ex(1:nx-1,2:ny,:)-Ex(1:nx-1,1:ny-1,:))*ds/2/dy-(Ey(2:nx,1:ny-1,:)-Ey(1:nx-1,1:ny-1,:))*ds/2/dx;


		!Hz(sx,sy,sz)=Hz(sx,sy,sz)+S;  !设置z方向的波源
		!stop   
		!------------------------------------------------------------------------------------------------

	end subroutine movon

	subroutine savedata(fname)
		implicit none
		character(256) fname
		outdata(1:nx-1,1:ny-1,1:nz-1)=real(Hx(1:nx-1,1:ny-1,1:nz-1)) 
		call matwrite(trim(fname),'Hxr',outdata,'w')

		outdata(1:nx-1,1:ny-1,1:nz-1)=imag(Hx(1:nx-1,1:ny-1,1:nz-1)) 
		call matwrite(trim(fname),'Hxi',outdata,'a')
		outdata(1:nx-1,1:ny-1,1:nz-1)=real(Hy(1:nx-1,1:ny-1,1:nz-1)) 
		call matwrite(trim(fname),'Hyr',outdata,'a')
		outdata(1:nx-1,1:ny-1,1:nz-1)=imag(Hy(1:nx-1,1:ny-1,1:nz-1)) 
		call matwrite(trim(fname),'Hyi',outdata,'a')

		outdata(1:nx-1,1:ny-1,1:nz-1)=real(Hz(1:nx-1,1:ny-1,1:nz-1)) 
		call matwrite(trim(fname),'Hzr',outdata,'a')
		outdata(1:nx-1,1:ny-1,1:nz-1)=imag(Hz(1:nx-1,1:ny-1,1:nz-1)) 
		call matwrite(trim(fname),'Hzi',outdata,'a')

		outdata(1:nx-1,1:ny-1,1:nz-1)=real(Ex(1:nx-1,1:ny-1,1:nz-1)) 
		call matwrite(trim(fname),'Exr',outdata,'a')
		outdata(1:nx-1,1:ny-1,1:nz-1)=imag(Ex(1:nx-1,1:ny-1,1:nz-1)) 
		call matwrite(trim(fname),'Exi',outdata,'a')

		outdata(1:nx-1,1:ny-1,1:nz-1)=real(Ey(1:nx-1,1:ny-1,1:nz-1)) 
		call matwrite(trim(fname),'Eyr',outdata,'a')
		outdata(1:nx-1,1:ny-1,1:nz-1)=imag(Ey(1:nx-1,1:ny-1,1:nz-1)) 
		call matwrite(trim(fname),'Eyi',outdata,'a')

		outdata(1:nx-1,1:ny-1,1:nz-1)=real(Ez(1:nx-1,1:ny-1,1:nz-1)) 
		call matwrite(trim(fname),'Ezr',outdata,'a')
		outdata(1:nx-1,1:ny-1,1:nz-1)=imag(Ez(1:nx-1,1:ny-1,1:nz-1)) 
		call matwrite(trim(fname),'Ezi',outdata,'a')
	end subroutine savedata


end module subroutines2