te1.f90

这是一个计算电磁场的程序 该程序已经调试通过了 可以直接计算 欢迎交流工作经验

program ex0101
implicit none

integer,parameter::imax=100                              !FDTD吸收边界
integer,parameter::imin=-100
integer,parameter::jmax=100
integer,parameter::jmin=-100

!integer::iflag=0                               !交替存储标志

real::ex(imin:imax,jmin:jmax-1)=0.0            !TE波分量的采样点的设置
real::ey(imin:imax-1,jmin:jmax)=0.0
real::hz(imin:imax,jmin:jmax)=0.0
 
real::hac(1:4,1)=0.0                           !吸收边界条件角点临时变量


real,parameter::pi=3.14159265
real,parameter::mu0=pi*4.e-7
real,parameter::eps0=8.85e-12
real z                                         !自由空间阻抗
real,parameter::wavelength=632.8e-9            !波长
!real,parameter::k=2*pi/wavelength              !波矢
real,parameter::v=3.0e+8                       !自由空间电磁波传播速度
real,parameter::omiga=2*pi*v/wavelength        !角频率
real,parameter::sqrt2=1.41421356

real cp
real sp

real,parameter::dx=10.0e-9                     !空间步长
real,parameter::dy=10.0e-9
real,parameter::dd=10.0e-9                     
real,parameter::dt=dd/(2*3.e+8)                !时间步长
real::t1,t2                                    !临时存储变量
integer i,j,n,ii,iii,nn
integer::k
integer,parameter::timestop=2000               !总时间步数
real::media(1:4)                               !介质参数数组
character*1 outflag                            !输出标志

real ca(-100:100,-100:100)                                                                                                                                                                                                                                                                                                                                                                                                    
real cb(-100:100,-100:100)
real cr(-100:100,-100:100)
real cq(-100:100,-100:100)

real::hz1(imin:imax,jmin:jmax)=0.0             !前一时间步电场值
real::hz2(imin:imax,jmin:jmax)=0.0
real::hz3(imin:imax,jmin:jmax)=0.0
real::hz4(imin:imax,jmin:jmax)=0.0


z=sqrt(mu0/eps0)


media(1)=1.                                   !自由空间参数设置
media(2)=1.
media(3)=0.
media(4)=0.

!***************************上半空间****************************

do i=-100,100,1                                
do j=-100,100,1
ca(i,j)=(1.-media(3)*dt/(2*media(1)*eps0))/&
                      (1.+media(3)*dt/(2*media(1)*eps0))
cb(i,j)=(dt/(media(1)*eps0))/(1.+media(3)*dt/(2*media(1)*eps0))
cr(i,j)=(1.-media(4)*dt/(2*media(2)*mu0))/&
                      (1.+media(4)*dt/(2*media(2)*mu0))
cq(i,j)=(dt/(media(2)*mu0))/(1.+media(4)*dt/(2*media(2)*mu0))
end do
end do

outflag='i'

!**************************************************************************

n=0
do while(n<=timestop)             !主循环
n=n+1
if(mod(n,100)==0)then
write(*,*) n
end if

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

!* * * * * * * * * * * * * HZ的FDTD迭代 * * * * * * * * * * * * * *

do i=imin+1,imax-1
do j=jmin+1,jmax-1
hz(i,j)=cr(i,j)*hz(i,j)-cq(i,j)*((ey(i,j)-ey(i-1,j))/dd-(ex(i,j)-ex(i,j-1))/dd)
end do
end do

!* * * * * * * * * * 吸收边界条件Hz * * * * * * * * * * * *

do i=imin+1,imax-1                         
hz(i,jmin)=hz1(i,jmin+1)+(v*dt-dd)/(v*dt+dd)*(hz(i,jmin+1)-hz(i,jmin))-v*v*eps0*dt/(2*(v*dt+dd))&
           *(ey(i,jmin+1)-ey(i-1,jmin+1)+ey(i,jmin)-ey(i-1,jmin))           !下边界
hz1(i,jmin+1)=hz(i,jmin+1)
hz(i,jmax)=hz3(i,jmax-1)+(v*dt-dd)/(v*dt+dd)*(hz(i,jmax-1)-hz(i,jmax))-v*v*eps0*dt/(2*(v*dt+dd))&
           *(ey(i,jmax-1)-ey(i-1,jmax-1)+ey(i,jmax)-ey(i-1,jmax))           !上边界
hz3(i,jmax-1)=hz(i,jmax-1)
end do

do j=jmin+1,jmax-1                          
hz(imin,j)=hz4(imin+1,j)+(v*dt-dd)/(v*dt+dd)*(hz(imin+1,j)-hz(imin,j))+v*v*eps0*dt/(2*(v*dt+dd))&
           *(ex(imin,j)-ex(imin,j-1)+ex(imin+1,j)-ex(imin+1,j-1))           !左边界
hz4(imin+1,j)=hz(imin+1,j)
hz(imax,j)=hz2(imax-1,j)+(v*dt-dd)/(v*dt+dd)*(hz(imax-1,j)-hz(imax,j))+v*v*eps0*dt/(2*(v*dt+dd))&
           *(ex(imax,j)-ex(imax,j-1)+ex(imax-1,j)-ex(imax-1,j-1))           !右边界
hz2(imax-1,j)=hz(imax-1,j)
end do

!* * * * * * * * * * 吸收边界角点Hz * * * * * * * * * * 

hz(imin,jmin)=hac(4,1)+(v*dt-sqrt2*dd)/(v*dt+sqrt2*dd)*(hz(imin+1,jmin+1)-hz(imin,jmin))
hac(4,1)=hz(imin+1,jmin+1)

hz(imax,jmax)=hac(2,1)+(v*dt-sqrt2*dd)/(v*dt+sqrt2*dd)*(hz(imax-1,jmax-1)-hz(imax,jmax))
hac(2,1)=hz(imax-1,jmax-1)

hz(imin,jmax)=hac(3,1)+(v*dt-sqrt2*dd)/(v*dt+sqrt2*dd)*(hz(imin+1,jmax-1)-hz(imin,jmax))
hac(3,1)=hz(imin+1,jmax-1)

hz(imax,jmin)=hac(1,1)+(v*dt-sqrt2*dd)/(v*dt+sqrt2*dd)*(hz(imax-1,jmin+1)-hz(imax,jmin))
hac(1,1)=hz(imax-1,jmin+1)

!***************************Ex,Ey的FDTD迭代**************************

do i=imin,imax,1
do j=jmin,jmax-1,1
ex(i,j)=ca(i,j)*ex(i,j)+cb(i,j)*(hz(i,j+1)-hz(i,j))/dy
end do
end do

do i=imin,imax-1,1
do j=jmin,jmax,1
ey(i,j)=ca(i,j)*ey(i,j)-cb(i,j)*(hz(i+1,j)-hz(i,j))/dx
end do
end do


end do                                      !时间步主循环结束




open(1,file='ex.txt',access='sequential')
open(2,file='ey.txt',access='sequential')
open(3,file='hz.txt',access='sequential')

do i=imin,imax,1
do j=jmin,jmax-1,1
write(1,*)ex(i,j)
end do
end do

do i=imin,imax-1,1
do j=jmin,jmax,1
write(2,*)ey(i,j)
end do
end do

do i=imin,imax,1
do j=jmin,jmax,1
write(3,*)hz(i,j)
end do
end do

close(1)
close(2)
close(3)

stop
end