618阅读
0回复

[求助]求助啊。。。谁帮我改下程序啊

谁帮我改下程序啊。。。万分感谢

UID ：81280

注册：2011-08-09
登录：2013-01-14
发帖：13
等级：仿真新人

0楼发表于: 2011-08-09 17:40:41

!******************************************************************************************************
!        二维 FDTD区域中设置一随时间为正弦周期变化的点波源
!       对于Ez的TM波(Hx.Hy.Ez一组) 设介质是无耗介质(自由空间)
!        在区域截断边界设为二阶Mur吸收边界，角点亦处理
!        画出TM波情形是正玹变化点源辐射时的相位分布
!       language:Fortran90
!           彭欢 2011.6.1
!*******************************************************************************************************
program FDTD2D
implicit none
!********************************************************************************************************
!         定义变量
!********************************************************************************************************
integer i,j,m,n                           !循环变量与时间步数
real*8 Ez(0:200,0:200),Hx(0:200,0:200),Hy(0:200,0:200)
real*8 Eac(4,2,0:1),Eab(1:199,1:4)        !Eab吸收边界条件临时变量 Eac吸收边界条件角点临时变量
real*8 CA,CB,CP,CQ        !FDTD迭代公式系数
real*8 CM,CN,C1                           !TM波Mur吸收边界条件系数
real*8 Omega,dt,dx,dy,delta       !角频率
integer Imax,Imin,Jmax,Jmin               !FDTD 区域总边界（吸收边界）
integer WL          !数值波长
real*8 Z          !波阻抗
real*8,parameter:: C=3.0e+8               !光速
real*8,parameter:: Pi=3.14159263          !圆周率
real*8,parameter:: wavelength=0.01        !波长
real*8,parameter:: Epso=8.85e-12          !介电常数
real*8,parameter:: Mu0=4*pi*1.e-7         !自由空间电磁参数：磁导率

write(*,*)"时间步数"
read(*,*)n

!open(1,file='Ez.dat')
!********************************************************************************************************
!           参数赋初始值
!********************************************************************************************************
Z=sqrt(Mu0/Epso)        !自由空间中波阻抗为377
WL=40
dx=wavelength/WL
dy=wavelength/WL
delta=dx
dt=0.5*dx/C                            !需要满足Courant二维稳定性条件
Omega=2*Pi*(C/wavelength)
CA=1.
CB=dt/Epso          !由自有空间介质中电损耗与磁损耗分别为0，求出的迭代系数CA,CB,CP,CQ
CP=1.
CQ=dt/Mu0
Hx=0.
Hy=0.
Ez=0.
Imax=200
Imin=0
Jmax=200
Jmin=0
Ez(100,100)=sin(Omega*dt)      !设置一随时间为正弦周期变化的点波源
CM=(c*dt-dx)/(c*dt+dx)
CN=c**2*Mu0*dt/(2*(c*dt+dx))/Z        !吸收边界条件系数。注意：CN计算时候徐除以自由空间的波阻抗常数377
C1=(c*dt-sqrt(2.)*delta)/(c*dt+sqrt(2.)*delta) !角点处理的系数
!*********************************************************************************************************
!
!            Ez的迭代
!
!*********************************************************************************************************
do m=1,n
!*********************************************************************************************************
!            迭代求解电场Ez
!*********************************************************************************************************
do i=Imin+1,Imax-1
  do j=Jmin+1,Jmax-1
   Ez(i,j)=CA*Ez(i,j)+CB*((Hy(i,j)-Hy(i-1,j))/dx-(Hx(i,j)-Hx(i,j-1))/dy) !（见课本2-3-12式）求Ez
  end do
enddo
!*********************************************************************************************************
!            迭代求解磁场Hx
!*********************************************************************************************************
!Ez(Imax,Jmax)
do i=Imin,Imax
  do j=Jmin,Jmax-1
   Hx(i,j)=CP*Hx(i,j)-CQ*(Ez(i,j+1)-Ez(i,j))/dy                         !（见课本2-3-10式）求Hx
  end do
end do
!*********************************************************************************************************
!            迭代求解磁场Hy
!*********************************************************************************************************
do i=Imin,Imax-1
  do j=Jmin,Jmax
   Hy(i,j)=CP*Hy(i,j)+CQ*(Ez(i+1,j)-Ez(i,j))/dx                         !（见课本2-3-11式）求Hy
  end do
end do

!以上三式为二维情况FDTD电磁场的时间推进计算公式
Ez(100,100)=sin(Omega*m*dt)            !在（100，100）处加入正弦点波源
!************************************************************************************************************
!            二阶Mur吸收边界
!************************************************************************************************************

!************************************************************************************************************
!            y=0与y=200边界
!************************************************************************************************************
do i=Imin+1,Imax-1
  Ez(i,Jmin)=Eab(i,1)+CM*(Ez(i,Jmin+1)-Ez(i,Jmin))+CN*(Hy(i,Jmin)-Hy(i-1,Jmin)+Hy(i,Jmin+1)-Hy(i-1,Jmin+1))
  Eab(i,1)=Ez(i,Jmin+1)
  Ez(i,Jmax)=Eab(i,3)+CM*(Ez(i,Jmax-1)-Ez(i,Jmax))+CN*(Hy(i,Jmax)-Hy(i-1,Jmax)+Hy(i,Jmax-1)-Hy(i-1,Jmax-1))
  Eab(i,3)=Ez(i,Jmax-1)
end do
!************************************************************************************************************
!            x=0与x=200边界
!************************************************************************************************************
do j=Jmin+1,Jmax-1
  Ez(Imin,j)=Eab(j,4)+CM*(Ez(Imin+1,j)-Ez(Imin,j))-CN*(Hx(Imin,j)-Hx(Imin,j-1)+Hx(Imin+1,j)-Hx(Imin+1,j-1))
  Eab(j,4)=Ez(1,j)
  Ez(Imax,j)=Eab(j,2)+CM*(Ez(Imax-1,j)-Ez(Imax,j))-CN*(Hx(Imax,j)-Hx(Imax,j-1)+Hx(Imax-1,j)-Hx(Imax-1,j-1))
  Eab(j,2)=Ez(Imax-1,j)

    end do
!*************************************************************************************************************
!            吸收边界角点处
!*************************************************************************************************************
Ez(Imin,Jmin)=Ez(Imin+1,Jmin+1)+C1*(Ez(Imin+1,Jmin+1)-Ez(Imin,Jmin))
Ez(Imax,Jmin)=Ez(Imax-1,Jmin+1)+C1*(Ez(Imax-1,Jmin+1)-Ez(Imax,Jmin))
Ez(Imax,Jmax)=Ez(Imax-1,Jmax-1)+C1*(Ez(Imax-1,Jmax-1)-Ez(Imax,Jmax))
Ez(Imax,Jmax)=Ez(Imin+1,Jmax-1)+C1*(Ez(Imin+1,Jmax-1))-Ez(Imax,Jmax)

!*********************************************************************************************************
!            多文件输出
!*********************************************************************************************************
100   format(1350f16.5)
if(m==300)THEN

      OPEN(1,File='Ez300.dat')
        do j=Jmin,Jmax
          write(1,100) (Ez(i,j),i=Imin,Imax)
        end do
      close(1)

elseif(m==500)then
  OPEN(1,File='Ez500.dat')
        do j=Jmin,Jmax
          write(1,100) (Ez(i,j),i=Imin,Imax)
        end do
  close(1)

elseif(m==1000)then
  OPEN(1,File='Ez1000.dat')
        do j=Jmin,Jmax
          write(1,100) (Ez(i,j),i=Imin,Imax)
        end do
      close(1)

elseif(m==3000)then
  OPEN(1,File='Ez3000.dat')
        do j=Jmin,Jmax
          write(1,100) (Ez(i,j),i=Imin,Imax)
        end do
      close(1)

elseif(m==5000)then
  OPEN(1,File='Ez5000.dat')
        do j=Jmin,Jmax
          write(1,100) (Ez(i,j),i=Imin,Imax)
        end do
      close(1)

elseif(m==6000)then
  OPEN(1,File='Ez6000.dat')
        do j=Jmin,Jmax
          write(1,100) (Ez(i,j),i=Imin,Imax)
        end ..

未注册仅能浏览部分内容，查看全部内容及附件请先登录或注册

发帖回复

[求助]求助啊 。。。谁帮我改下程序啊

[求助]求助啊。。。谁帮我改下程序啊