二维传播.f90

利用时域有限差分法模拟二维坐标中波传播

program main  ! 	TM mode	仅水平入射的情况
implicit none 
real,parameter::pi=3.1415927
real::Ez(-160:160,-160:160)	 ! Ez采样点
real::Hx(-160:160,-160:159)	 ! Hx采样点
real::Hy(-160:159,-160:160)	 ! Hy采样点
real::FE,FH                  ! 电场磁场迭代系数
real::Eab(-160:160,4)        ! 吸收边界临时变量
real::Ezcornor(4)			 ! 吸收边界角点临时变量
integer::timestop,wl		 ! 总时间步,波长划分数目
integer::Imax,Imin,Jmax,Jmin ! 区域总边界(吸收边界) 
real::wavelength	         ! 自由空间波数，波阻抗，波长
real::t1 		    	     ! 临时变量
real::fa1,fa2,fcornor		 ! 吸收系数
real::mu0,empsi0,omiga
real::dx,dt,c0
integer::i,j,n
open(1,file='二维传播Ez.dat')
open(2,file='二维传播Hx.dat')
open(3,file='二维传播Hy.dat')
open(16,file='二维Ezt.dat')
open(18,file='二维Ez2.dat')
Imax=150
Imin=-150
Jmax=150
Jmin=-150
timestop=302
mu0=4*pi*1e-7
empsi0=8.85*1e-12
c0=1./sqrt(mu0*empsi0)
wl=10
omiga=1.e10
wavelength=2*pi*c0/omiga
dx=wavelength/wl	   !也是y的离散
dt=0.5*dx/c0
fa1=(c0*dt-dx)/(c0*dt+dx)		    ! 吸收边界的电场迭代系数
fa2=0.5*(c0*c0*mu0*dt)/(c0*dt+dx)	! 吸收边界的磁场迭代系数
fcornor=(c0*dt-1.4142*dx)/(c0*dt+1.4142*dx)
FE=dt/(empsi0*dx)
FH=dt/(mu0*dx)

! *********Main loop **************
n=0
do while(n<timestop)
t1=n
if(int((t1+1)/timestop*100)>int(t1/timestop*100)) &
 write(*,*) 'It is computing',int(t1/timestop*100),'%'

!=========Ez的FDTD迭代===========
do i=imin+1,imax-1
 do j=jmin+1,jmax-1
 t1=hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1)
 ez(i,j)=ez(i,j)+FE*t1
 enddo
enddo

  ! -------- 吸收边界条件Ez  ---------------
do i=imin+1,imax-1
 t1=hy(i,jmin)-hy(i-1,jmin)+hy(i,jmin+1)-hy(i-1,jmin+1)
 ez(i,jmin)=eab(i,1)+fa1*(ez(i,jmin+1)-ez(i,jmin))+fa2*t1
 eab(i,1)=ez(i,jmin+1)
 t1=hy(i,jmax)-hy(i-1,jmax)+hy(i,jmax-1)-hy(i-1,jmax-1)
 ez(i,jmax)=eab(i,3)+fa1*(ez(i,jmax-1)-ez(i,jmax))+fa2*t1
 eab(i,3)=ez(i,jmax-1)
enddo

do j=jmin+1,jmax-1
 t1=hx(imax,j)-hx(imax,j-1)+hx(imax-1,j)-hx(imax-1,j-1)
 ez(imax,j)=eab(j,2)+fa1*(ez(imax-1,j)-ez(imax,j))-fa2*t1
 eab(j,2)=ez(imax-1,j)
 t1=hx(imin,j)-hx(imin,j-1)+hx(imin+1,j)-hx(imin+1,j-1)
 ez(imin,j)=eab(j,4)+fa1*(ez(imin+1,j)-ez(imin,j))-fa2*t1
 eab(j,4)=ez(imin+1,j)
enddo
  ! ------------吸收边界角点Ez-----------
Ez(imin,jmin)=ezcornor(1)+fcornor*(ez(imin+1,jmin+1)-ez(imin,jmin))
Ez(imax,jmin)=ezcornor(2)+fcornor*(ez(imax-1,jmin+1)-ez(imax,jmin))
Ez(imax,jmax)=ezcornor(3)+fcornor*(ez(imax-1,jmax-1)-ez(imax,jmax))
Ez(imin,jmax)=ezcornor(4)+fcornor*(ez(imin+1,jmax-1)-ez(imin,jmax))
ezcornor(1)=ez(imin+1,jmin+1)
ezcornor(2)=ez(imax-1,jmin+1)
ezcornor(3)=ez(imax-1,jmax-1)
ezcornor(4)=ez(imin+1,jmax-1)

Ez(0,0)=sin(omiga*n*dt)

! =========Hx,Hy 的FDTD迭代==========
do i=imin,imax
 do j=jmin,jmax-1
 t1=ez(i,j+1)-ez(i,j)
 hx(i,j)=hx(i,j)-FH*t1
 enddo
enddo

do i=imin,imax-1
 do j=jmin,jmax
 t1=ez(i+1,j)-ez(i,j)
 hy(i,j)=hy(i,j)+FH*t1
 enddo
enddo  

write(16,*) n,Ez(0,0),Ez(Imax,0),Ez(20,0),Ez(0,20)  
n=n+1
enddo !endwhile     时间步主循环结束

write(*,*) char(7)

do i=Imin,Imax
 write(18,*) i,Ez(i,20),Ez(0,i)
enddo

! ========== 记录空间全部的数据	=======
do i=-40,40
 do j=-40,40
  write(1,*) i,j,ez(i,j)
 enddo
enddo

do i=imin,imax
 do j=jmin,jmax-1
  write(2,*) i,j,hx(i,j)
 enddo
enddo
 
do i=imin,imax-1
 do j=jmin,jmax
  write(3,*) i,j,hy(i,j)
 enddo
enddo
 
close(1)
close(2)
close(3)
close(18)
close(16)
end program main