! TM波FDTD讨论平面波源的加入module data_moduleimplicit noneinteger,parameter::nx0=0,nx1=360,ny0=0,ny1=360,nz0=-100,nz1=1200integer,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.14159real,parameter::w=2*pi*f,s=-0.477369real,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.0real cez,chx,chyinteger,parameter::nt=2000,m0=200integer ncomplex 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_moduleimplicit noneinteger i,j,kreal t,dt=n*delttEi1=Ei0do k=nz0,nz1-1Hi0(k)=Hi0(k)-p1*(Ei1(k+1)-Ei1(k))end do!Ezido i=nxl1,nxl2do j=nyl1,nyl2d=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 doend dodo k=nz0+1,nz1-1Ei0(k)=Ei1(k)-p2*(Hi0(k)-Hi0(k-1)) ! 入射波的场量end doEi0(m0-2)=sin(w*t) !exp(-(4*pi*(t-t0)**2/Tal/tal))k=nz1Ei0(k)=Ei1(k-1)-1.0/3.0*(Ei0(k-1)-Ei1(k)) !入射波右吸收边界K=nz0Ei0(k)=Ei1(k+1)-1.0/3.0*(Ei0(k+1)-Ei1(k)) !入射波左吸收边界!Hxij=nyl1do i=nxl1,nxl2d=real(i-nxl1)*cos(fai)+(j-0.5-nyl1)*sin(fai)+0.5Hxi(i,j-1)=((d-int(d))*Hi0(m0-1+int(d)+1)+(1-(d-int(d)))*Hi0(m0-1+int(d)))*sin(fai) end doj=nyl2do i=nxl1,nxl2d=real(i-nxl1)*cos(fai)+(j+0.5-nyl1)*sin(fai)+0.5Hxi(i,j)=((d-int(d))*Hi0(m0-1+int(d)+1)+(1-(d-int(d)))*Hi0(m0-1+int(d)))*sin(fai) end do!Hyii=nxl1do j=nyl1,nyl2d=(i-0.5-nyl1)*cos(fai)+real(j-nyl1)*sin(fai)+0.5Hyi(i-1,j)=((d-int(d))*Hi0(m0-1+int(d)+1)+(1-(d-int(d)))*Hi0(m0-1+int(d)))*(-cos(fai)) enddoi=nxl2do j=nyl1,nyl2d=(i+0.5-nyl1)*cos(fai)+real(j-nyl1)*sin(fai)+0.5Hyi(i,j)=((d-int(d))*Hi0(m0-1+int(d)+1)+(1-(d-int(d)))*Hi0(m0-1+int(d)))*(-cos(fai)) enddowrite(13,*)n,Ei0(0),Ei0(20),Ei0(30),Ei0(100),sin(w*t)end subroutine inc!/////////////////////////////////////////////////////////////////////////////////////////////////program mainuse data_moduleimplicit noneinteger i,j,kreal topen(11,file='A.txt')open(12,file='fai.txt')open(13,file='Ez.txt',recl=300)Ez=0;Hx=0;Hy=0;Ez1=0Ei0=0;Hi0=0Ez4=0;Ez2=0;Ez3=cmplx(0.0,0.0)cez=1.0;chx=sin(fai);chy=-cos(fai)do n=0,ntcall inc()call FDTD()print *,'step',n,Ez(nx1/2,ny1/2)!write(13,*)n,Ei0(0),Ei0(20),Ei0(30),Ei0(100),sin(w*t)Ez1=Ezend doEz4=Ezdo n=nt+1,nt+42call inc()call FDTD()Ez1=Ezend doEz2=Ezcall AP()do j=ny0,ny1do i=nx0,nx1write(11,*)i,j,Ez0(i,j)write(12,*)i,j,sita(i,j)end doend doclose(11)close(12)close(13)end program main!///////////////////////////////////////////////////////////////////////////////////////////////// subroutine FDTD()use data_moduleimplicit noneinteger i,jcall cal_H()call cal_E()end subroutine FDTD!///////////////////////////////////////////////////////////////////////////////////////////////// subroutine cal_H()use data_moduleimplicit noneinteger i,jreal,external::Eido j=ny0,ny1-1Hx(i,j)=Hx(i,j)-p1*(Ez(i,j+1)-Ez(i,j))end doend dodo j=ny0,ny1do i=nx0,nx1-1Hy(i,j)=Hy(i,j)+p1*(Ez(i+1,j)-Ez(i,j))end doend do!连接边界i=nxl1-1do j=nyl1,nyl2Hy(i,j)=Hy(i,j)-p1*Ezi(i+1,j)end doi=nxl2do j=nyl1,nyl2Hy(i,j)=Hy(i,j)+p1*Ezi(i,j)end doj=nyl1-1do i=nxl1,nxl2Hx(i,j)=Hx(i,j)+p1*Ezi(i,j+1)end doj=nyl2do i=nxl1,nxl2Hx(i,j)=Hx(i,j)-p1*Ezi(i,j)end doend subroutine cal_H!///////////////////////////////////////////////////////////////////////////////////////////////// subroutine cal_E()use data_moduleimplicit noneinteger i,jreal,external::Hido j=ny0+1,ny1-1do i=nx0+1,nx1-1Ez(i,j)=Ez(i,j)+p2*(Hy(i,j)-Hy(i-1,j)-Hx(i,j)+Hx(i,j-1))end doend do!Ez(nx1/2,ny1/2)=Ez(nx1/2,ny1/2)-p2/delt*sin(w*(n+0.5)*deltt) call mur()! 连接边界!四条边do j=nyl1+1,nyl2-1Ez(i,j)=Ez(i,j)-p2*Hyi(i-1,j)end doi=nxl2do j=nyl1+1,nyl2-1Ez(i,j)=Ez(i,j)+p2*Hyi(i,j)end doj=nyl1do i=nxl1+1,nxl2-1Ez(i,j)=Ez(i,j)+p2*Hxi(i,j-1)end doj=nyl2do i=nxl1+1,nxl2-1Ez(i,j)=Ez(i,j)-p2*Hxi(i,j)end do!四个角i=nxl1j=nyl1Ez(i,j)=Ez(i,j)+p2*(Hxi(i,j-1)-Hyi(i-1,j))i=nxl2Ez(i,j)=Ez(i,j)+p2*(Hxi(i,j-1)+Hyi(i,j))j=nyl2Ez(i,j)=Ez(i,j)+p2*(-Hxi(i,j)+Hyi(i,j))i=nxl1Ez(i,j)=Ez(i,j)+p2*(-Hxi(i,j)-Hyi(i-1,j))end subroutine cal_E!/////////////////////////////////////////////////////////////////////////////////////////////////subroutine mur()use data_moduleimplicit noneinteger i,j!四条边i=nx0do j=ny0+1,ny1-1Ez(i,j)=Ez1(i+1,j)+p*(Ez(i+1,j)-Ez1(i,j))+q*(Hx(i+1,j)+Hx(i,j)-Hx(i+1,j-1)-Hx(i,j-1)) end doi=nx1-1do j=ny0+1,ny1-1Ez(i+1,j)=Ez1(i,j)+p*(Ez(i,j)-Ez1(i+1,j))+q*(Hx(i+1,j)+Hx(i,j)-Hx(i+1,j-1)-Hx(i,j-1)) end doj=ny0do i=nx0+1,nx1-1Ez(i,j)=Ez1(i,j+1)+p*(Ez(i,j+1)-Ez1(i,j))-q*(Hy(i,j+1)+Hy(i,j)-Hy(i-1,j+1)-Hy(i-1,j)) end doj=ny1-1do i=nx0+1,nx1-1Ez(i,j+1)=Ez1(i,j)+p*(Ez(i,j)-Ez1(i,j+1))-q*(Hy(i,j+1)+Hy(i,j)-Hy(i-1,j+1)-Hy(i-1,j)) end do!角点i=nx0j=ny0Ez(i,j)=Ez1(i+1,j+1)+s*(Ez(i+1,j+1)-Ez1(i,j))i=nx1Ez(i,j)=Ez1(i-1,j+1)+s*(Ez(i-1,j+1)-Ez1(i,j))j=ny1Ez(i,j)=Ez1(i-1,j-1)+s*(Ez(i-1,j-1)-Ez1(i,j))i=nx0Ez(i,j)=Ez1(i+1,j-1)+s*(Ez(i+1,j-1)-Ez1(i,j))end subroutine mur!/////////////////////////////////////////////////////////////////////////////////////////////////////// subroutine AP()use data_moduleimplicit noneinteger i,jreal mreal tmp1,tmp2,tmp3complex tmptmp=cmplx(Ez2(nx1/2,ny1/2),Ez4(nx1/2,ny1/2))do j=ny0,ny1do i=nx0,nx1Ez3(i,j)=cmplx(Ez2(i,j),Ez4(i,j))Ez0(i,j)=sqrt(Ez2(i,j)**2+Ez4(i,j)**2)m=Imag(Ez3(i,j)/tmp)/real(Ez3(i,j)/tmp)if(real(Ez3(i,j)/tmp)==0 .and. Imag(Ez3(i,j)/tmp)>0) thensita(i,j)=pi/2elseif(real(Ez3(i,j)/tmp)==0 .and. Imag(Ez3(i,j)/tmp)<0) thensita(i,j)=-pi/2elsesita(I,j)=atan(m)end ifif (abs(sita(i,j))<=1e-6) thensita(i,j)=0end ifend doend doend subroutine AP!///////////////////////////////////////////////////////////////////////////////////////////////////////。