medium.f90

利用时域有限差分法计算有介质板的高斯波源

! ========================================
! ==============高斯波源有介质板的FDTD====
! =================厚度为9.00厘米=========
! ========================================
program main    
implicit none 
integer,parameter::Imax=250,Tmax=160
integer,parameter::obStart=100
complex(8)::cei(tmax),cext(tmax)
real(8)::ex(Imax),hy(Imax-1),mu0,empsi0,mur1,empsir1
real(8)::pi=3.14159265d0,mu1,empsi1,c1,c0,ex1,ex2 !上时刻的Ex(Imax-1)值
real(8)::fa,FH(0:1),FE(0:1),dt,dz,tao	!中间系数
real(8)::obz,fmax !介质板的厚度,关心的最大频率
integer::i,n,m,ob(Imax),obEnd
open(1,file='高斯波ns源ext.dat')
open(2,file='高斯波ns源ex.dat')
open(3,file='高斯波ns源hy.dat')
open(5,file='反射频谱.dat')

m=40
fmax=10d9
obz=9.d-2
mu0=4*pi*1.d-7
empsi0=8.85*1.d-12
mur1=1.
empsir1=4.
mu1=mu0*mur1
empsi1=empsi0*empsir1
c0=3d8   !真空中速度
c1=c0/dsqrt(mur1*empsir1)	  !介质板中的速度
dz=min(c0,c1)/fmax
dz=dz/20.
obend=obz/dz
print*,obz/dz
if(obz/dz-obend>=0.5) obend=obend+1
obend=obstart+obend
dt=dz/max(c0,c1)
tao=m*dt

if(obend>Imax-5) then
  print*,'空间网格太短了,请扩大!'
  write(*,*) char(7)
  write(*,*) char(7)
  write(*,*) char(7)
  stop
endif
print*,obstart,obend
ob=0
do i=1,Imax
 if(i>=obstart.and.i<=obend) ob(i)=1
enddo

FH(0)=dt/(mu0*dz)
FH(1)=FH(0)/mur1
FE(0)=dt/(empsi0*dz)
FE(1)=FE(0)/empsir1
fa=(c0*dt-dz)/(c0*dt+dz)
ex2=ex(Imax-1)

do n=1,Tmax
 
 do i=1,Imax-1
   hy(i)=hy(i)-FH(ob(i))*(ex(i+1)-ex(i))
 enddo
 
 do i=2,Imax-1
   if(i==obstart.or.i==obend) then
	ex(i)=ex(i)-2*dt/((empsi0+empsi1)*dz)*(hy(i)-hy(i-1))
   else
    ex(i)=ex(i)-FE(ob(i))*(hy(i)-hy(i-1))
   endif	 
 enddo
 ex(Imax)=ex2+fa*(ex(Imax-1)-ex(Imax)) !右边吸收边界
 ex2=ex(Imax-1)
 
 if(n>3*m) then      !将波源处变为吸收边界
  ex(1)=ex1+fa*(ex(2)-ex(1))
  ex1=ex(2)
 else
   ex(1)=1.0*dexp(-4.*pi*(float(n)/float(m)-1.5)**2)
   cei(n)=ex(1)
   ex1=ex(2)
 endif
 
 cext(n)=ex(obstart-m)
 write(1,*) n,ex(obstart-m)

enddo !时间迭代结束
call dft(cext,tmax)
!call dft(cei,tmax)
do i=1,tmax/2
   write(5,*) i,abs(cext(i))
enddo
do i=1,Imax
 write(2,*) i,ex(i)
enddo

do i=1,Imax-1
 write(3,*) i,hy(i)
enddo

close(1)
close(2)
close(3)
end program main

subroutine dft(p,nn)
implicit none
integer::n,m,nn
complex(8)::p(nn),q(nn),ii,csum
real(8)::pi=3.1415927d0
ii=(0,1.)
do m=1,nn
   csum=0
   do n=0,nn-1
      csum=csum+p(n+1)*exp(-ii*2*pi*m*n/nn)
   enddo
   q(m)=csum
enddo
p=q
end subroutine dft