二维TM波讨论平面波源(使用直接算方法)的加入
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
! TM波FDTD讨论平面波源的加入
module data_module
implicit none
integer,parameter::nx0=0,nx1=360,ny0=0,ny1=360,nz0=-100,nz1=1200
integer,parameter::nxl1=nx0+80,nxl2=nx1-80,nyl1=ny0+80,nyl2=ny1-80 !连接边界
real,parameter::f=2.0e8,c=3.0e8,delt=0.0177,deltt=delt/6.0e8,eps0=8.85e-12,miu0=1.2566e-6,pi= 3.14159
real,parameter::w=2*pi*f,s=-0.477369
real,parameter::p=-1.0/3.0,q=-miu0*c/6,r=-miu0*c/2,p1=1/(2*miu0*c),p2=1/(2*eps0*c)
real,parameter::tal=2e-9,t0=0.8*tal,fai=pi/3.0
real cez,chx,chy
integer,parameter::nt=2000,m0=200
integer n
complex Ez3(nx0:nx1,ny0:ny1)
real Ez4(nx0:nx1,ny0:ny1),Ez2(nx0:nx1,ny0:ny1) !记录幅值提取时的实部和虚部
real sita(nx0:nx1,ny0:ny1),Ez0(nx0:nx1,ny0:ny1) !记录幅值和相位
real Ez(nx0:nx1,ny0:ny1),Hx(nx0:nx1,ny0:ny1),Hy(nx0:nx1,ny0:ny1),Ez1(nx0:nx1,ny0:ny1) real Ei0(nz0:nz1),Hi0(nz0:nz1),Ei1(nz0:nz1)
real Ezi(nx0:nx1,ny0:ny1),Hxi(nx0:nx1,ny0:ny1),Hyi(nx0:nx1,ny0:ny1)
end module data_module
!/////////////////////////////////////////////////////////////////////////////////////////////////
subroutine inc()
use data_module
implicit none
integer i,j,k
real t,d
t=n*deltt
Ei1=Ei0
do k=nz0,nz1-1
Hi0(k)=Hi0(k)-p1*(Ei1(k+1)-Ei1(k))
end do
!Ezi
do i=nxl1,nxl2
do j=nyl1,nyl2
d=real(i-nxl1)*cos(fai)+real(j-nyl1)*sin(fai)
Ezi(i,j)=(d-int(d))*Ei0(m0+int(d)+1)+(1-(d-int(d)))*Ei0(m0+int(d))
end do
end do
do k=nz0+1,nz1-1
Ei0(k)=Ei1(k)-p2*(Hi0(k)-Hi0(k-1)) ! 入射波的场量
end do
Ei0(m0-2)=sin(w*t) !exp(-(4*pi*(t-t0)**2/Tal/tal))
k=nz1
Ei0(k)=Ei1(k-1)-1.0/3.0*(Ei0(k-1)-Ei1(k)) !入射波右吸收边界
K=nz0
Ei0(k)=Ei1(k+1)-1.0/3.0*(Ei0(k+1)-Ei1(k)) !入射波左吸收边界
!Hxi
j=nyl1
do i=nxl1,nxl2
d=real(i-nxl1)*cos(fai)+(j-0.5-nyl1)*sin(fai)+0.5
Hxi(i,j-1)=((d-int(d))*Hi0(m0-1+int(d)+1)+(1-(d-int(d)))*Hi0(m0-1+int(d)))*sin(fai) end do
j=nyl2
do i=nxl1,nxl2
d=real(i-nxl1)*cos(fai)+(j+0.5-nyl1)*sin(fai)+0.5
Hxi(i,j)=((d-int(d))*Hi0(m0-1+int(d)+1)+(1-(d-int(d)))*Hi0(m0-1+int(d)))*sin(fai) end do
!Hyi
i=nxl1
do j=nyl1,nyl2
d=(i-0.5-nyl1)*cos(fai)+real(j-nyl1)*sin(fai)+0.5
Hyi(i-1,j)=((d-int(d))*Hi0(m0-1+int(d)+1)+(1-(d-int(d)))*Hi0(m0-1+int(d)))*(-cos(fai)) enddo
i=nxl2
do j=nyl1,nyl2
d=(i+0.5-nyl1)*cos(fai)+real(j-nyl1)*sin(fai)+0.5
Hyi(i,j)=((d-int(d))*Hi0(m0-1+int(d)+1)+(1-(d-int(d)))*Hi0(m0-1+int(d)))*(-cos(fai)) enddo
write(13,*)n,Ei0(0),Ei0(20),Ei0(30),Ei0(100),sin(w*t)
end subroutine inc
!/////////////////////////////////////////////////////////////////////////////////////////////////
program main
use data_module
implicit none
integer i,j,k
real t
open(11,file='A.txt')
open(12,file='fai.txt')
open(13,file='Ez.txt',recl=300)
Ez=0;Hx=0;Hy=0;Ez1=0
Ei0=0;Hi0=0
Ez4=0;Ez2=0;Ez3=cmplx(0.0,0.0)