mur_3dfdtd.f90

fortran程序计算光子晶体

!===========================================================
! 3D FDTD with 1st Mur boundary conditions
! Written by Bing Wang, Department of Physics, 
! Wuhan University, Wuhan 430072
! Don't hesitate to contact me at [e-mail] if 
! there are any questions in calculation with this program.
!===========================================================  
program main
	use constants
	use subroutines
	implicit none
	
	integer i,j,k,n,si,sj,sk,x1,x2,y1,y2,z1,z2,width,mwidth,thick,mthick
	real(kind) U,time_begin,time_begin_loop,time_end
	complex(kind) S
	character(256) filename
	integer xc,yc,zc ! the origin of the coordinates 

!==>	

	print *,"==============3D-FDTD=================="
	
	lamda=539.1d0*nm
	w=2d0*pi*c/lamda

	ds=10d0*nm
	dt=ds/(2d0*c)
	dx=ds
	dy=ds
	dz=ds*4
	
	nx=81
	ny=101
	nz=101
	nt=2000
	
	xc=(1+nx)/2
	yc=(1+ny)/2
	zc=11


	width=3
	!mwidth=6
	thick=3
	mthick=4

	si=8
	sj=8
	sk=10
	
	filename="coupling"

	mat_num=3
	
	call initialization() 
	
	!material 0
	e_p(0)=e0
	u_p(0)=u0
	g_p(0)=g0
	!material 1
	e_p(1)=e0
	u_p(1)=u0
	g_p(1)=cmplx(-42.13d0,11.96d0)*cmplx(0d0,-w*e0) ! metal 1 <-> Al
	!material 2
	e_p(2)=e0
	u_p(2)=u0
	g_p(2)=cmplx(-10.55d0, 0.84d0)*cmplx(0d0,-w*e0) ! metal 2 <-> Ag

	e_p(3)=1.5d0**2*e0
	u_p(3)=u0
	g_p(3)=g0 !                                       metal 3 <-> Si

    do i=0,mat_num
		c_p(i)=const_c(e_p(i),u_p(i),g_p(i))
		d_p(i)=const_d(e_p(i),u_p(i),g_p(i))
	end do

	CEx=c_p(0);DEx=d_p(0);
	CEy=c_p(0);DEy=d_p(0);
	CEz=c_p(0);DEz=d_p(0);

	!=========================================================== 
	! set the structure
	!===========================================================
	z1=11
	z2=nz-10
	x1=xc-10-thick;  x2=xc-10+thick; 
	call rectangle(x1,x2,5,ny-4,z1,z2,3)
	
	y1=yc-width;  y2=yc+width;
	call rectangle(x1,x2,y1,y2,z1,z2,2)

	x1=xc+10-thick;  x2=xc+10+thick; 
	call rectangle(x1,x2,5,ny-4,z1,z2,1)
	
	y1=yc-width;  y2=yc+width;
	call rectangle(x1,x2,y1,y2,z1,z2,2)
	!-----------------------------------------------------------


	call cpu_time(time_begin)
	time_begin_loop=time_begin

	print *,"initialization completed"
	
	n=0
	fdtd_loop: do

		n=n+1

		!======================================================
		!	set open-close function
		!======================================================
		if(n<pi/w/dt) then
			U=0.5d0*(1d0-cos(pi*n/(pi/w/dt)))
		else 
			U=1d0;
		end if
		!======================================================
		!	mov on
		!======================================================
		call movon()

		!======================================================
		!	set wave source
		!======================================================
		S=cmplx(cos(n*w*dt),-sin(n*w*dt))

		Hx(x1:x2,y1:y2,sk)=Hx(x1:x2,y1:y2,sk)+S*U


		!======================================================
		
		if(mod(n,100)==0)	then 
			call cpu_time(time_end)
			write(*,"('step',I6,' has completed: ', F8.3,' s')") n, time_end-time_begin_loop
			time_begin_loop=time_end
		endif
		if (n>=nt) exit fdtd_loop

	end do fdtd_loop
	!=====================================================
	!    save data
	!=====================================================
	call cpu_time(time_end)

	write(*,"('Total operation time is ',F8.3,' minutes!')") (time_end-time_begin)/60d0
	
	!call makedirqq(trim(filename))
	!call savedata(trim(filename//"\result.mat"))

	call savedata(filename)
	
	write(*,"(/,'result has beens saved',/)")
	
	call freespace()
	
end program main