[求助]谁有FDTD模拟TM传播的程序啊

谁有FDTD模拟TM传播的程序啊 Io8h 8N-  
谢谢 .. m6'VMW  
vUlGE  

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

  今天刚做好的二维TM   UPML边界的程序。 \*r]v;NcP  
（注释掉的部分是Mur吸收边界及散射边界处理程序） dp=#|!jc  
yrCY-'%  
rd0[(-  
Two_dimension_TM.rar (4 K) 下载次数:52 <ZwmXD.VD  
4/+P7.}ea-  

图片:Ex.jpg

引用第1楼shuangzh于2010-11-22 19:59发表的  : NTHy!y<!h  
  今天刚做好的二维TM   UPML边界的程序。 W|PKcZ ]Uc  
（注释掉的部分是Mur吸收边界及散射边界处理程序）
[图片] nj7wc9z4  
....... u,d@oF(=  

-}Jf4k#G  
;!CYp;_  
t^=6czk  
bIP'(B#1K  
谢谢，我也刚刚做了一个Mur二阶的吸收边界，稍后把图发上来，UPML还没有做，参考一下你的啊 ;
plzJ6>  
3q [S}o[v\  
;btH[a iV  
Q&Rj)1!  
!~{AF|2f  

不错，不错！

为什么我下载不了程序啊

我编的二维FDTD   TM波点源辐射为什么将结果不好，哪位高手能帮我指点下 啊 9H:5XR
 
!****************************************************************************************************** 4mSL*1j  
!                             二维 FDTD区域中设置一随时间为正弦周期变化的点波源 b{qeu$GR  
!                            对于Ez的TM波(Hx.Hy.Ez一组)     设介质是无耗介质(自由空间) C-4NiXa  
!                             在区域截断边界设为二阶Mur吸收边界，角点亦处理 {^gbS  
!                             画出TM波情形是正玹变化点源辐射时的相位分布 jXq~ x"(  
!                            language:Fortran90 |j53'>N[  
!                                        !******************************************************************************************************* e_s9E{(  
program  FDTD2D K}ACZT)Wp  
implicit none 6T{Zee  
+N1oOcPC>C  
!******************************************************************************************************** 5\1Z"?  
!                                    定义变量 g{wIdV  
!******************************************************************************************************** <r]7xsr  
 ^! /7  
integer i,j,m,n                           !循环变量与时间步数 MVHj?  
real*8 Ez(0:200,0:200),Hx(0:200,0:200),Hy(0:200,0:200) c\~H_ ~F  
real*8 Eac(4,2,0:1),Eab(1:199,1:4)        !Eab吸收边界条件临时变量  Eac吸收边界条件角点临时变量 }LQ*vD-Jj  
real*8 CA,CB,CP,CQ                          !FDTD迭代公式系数 +cv7]  
real*8 CM,CN,C1                           !TM波Mur吸收边界条件系数 e\i K  
real*8 Omega,dt,dx,dy,delta                      !角频率 T5_z^7d  
integer Imax,Imin,Jmax,Jmin               !FDTD    区域总边界（吸收边界） 2xRb$QF  
integer WL                                  !数值波长 %0T/>:1[E  
real*8 Z                                  !波阻抗 iOz<n z  
real*8,parameter:: C=3.0e+8               !光速 <Ur(< WTV  
real*8,parameter:: Pi=3.14159263          !圆周率 2Cn^<(F^4I  
real*8,parameter:: wavelength=0.01        !波长 JG&`l{c9  
real*8,parameter:: Epso=8.85e-12          !介电常数 %||}WT-wv  
real*8,parameter:: Mu0=4*pi*1.e-7         !自由空间电磁参数：磁导率 B%!z7AT  
Z0T{1YEJ  
|5)~WoV/G  
3;=nQ{0b  
f'aQ T  
write(*,*)"时间步数"                                 ;;'b;,/  
read(*,*)n 9i@AOU  
5zG6V2  
i$[wgvJIV  
R_J=x  
!open(1,file='Ez.dat') 5(bG  
m)9N9Ii#)  
~d6_  
!******************************************************************************************************** 7
?j$Lwt  
!                                            参数赋初始值                                                             l!e8=QlJ  
!******************************************************************************************************** "Q9S<O8)  
Z=sqrt(Mu0/Epso)                       !自由空间中波阻抗为377 4S|! iOY  
WL=40 'JY*K:-  
dx=wavelength/WL                       jQ.]m  
dy=wavelength/WL ,}F2l|x_  
delta=dx         9PhdoREb  
dt=0.5*dx/C                            !需要满足Courant二维稳定性条件 tVQfR*=  
Omega=2*Pi*(C/wavelength) p3O%|)yV  
CA=1. }/BwFB+(/  
CB=dt/Epso                               !由自有空间介质中电损耗与磁损耗分别为0，求出的迭代系数CA,CB,CP,CQ s`Fv!  
CP=1. m! &bK5+*  
CQ=dt/Mu0 K6=-Zf  
Hx=0. R-]QU`c  
Hy=0. vfzGRr  
Ez=0.                                     u.iFlU  
Imax=200 ^qqP):0y1V  
Imin=0 ;>[).fX>/  
Jmax=200 M`\c'|i/  
Jmin=0 XPXC7_fV  
Ez(100,100)=sin(Omega*dt)                        !设置一随时间为正弦周期变化的点波源 YXXUYi~!f  
CM=(c*dt-dx)/(c*dt+dx) p%Ae"#_X%  
CN=c**2*Mu0*dt/(2*(c*dt+dx))/Z           !吸收边界条件系数。注意：CN计算时候徐除以自由空间的波阻抗常数377 e3YZ-w^W~h  
C1=(c*dt-sqrt(2.)*delta)/(c*dt+sqrt(2.)*delta) !角点处理的系数 K !8+~[  
!********************************************************************************************************* bT-(lIU  
!         0K7]<\)  
!                                                Ez的迭代 At%g^  
!                                                             =N01!?{  
!********************************************************************************************************* (FZL>  
^ cN-   
do m=1,n *JGm  
!********************************************************************************************************* b_ Sh#d&  
!                                                迭代求解电场Ez _
dppUUm  
!********************************************************************************************************* Pgf$GXE  
ba|x?kz  
    do i=Imin+1,Imax-1 M]!\X6<_  
        do j=Jmin+1,Jmax-1 ;&e5.K+.Z  
            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 w(`X P  
        end do Mo &I
a6^  
    enddo ,HS\(Z  
!*********************************************************************************************************
!.iu_xJ  
!                                                迭代求解磁场Hx
5b9_6L6  
!********************************************************************************************************* &q1(v3cOO  
1iaNb[:QX  
!Ez(Imax,Jmax) 1JgnuBX"  
    do i=Imin,Imax UV)[a%/SB&  
        do j=Jmin,Jmax-1 Q@%VJPLv.  
            Hx(i,j)=CP*Hx(i,j)-CQ*(Ez(i,j+1)-Ez(i,j))/dy                         !（见课本2-3-10式）求Hx lT$Vv=M  
        end do C]*9:lK  
    end do srx`" :  
!********************************************************************************************************* "m})~va  
!                                                迭代求解磁场Hy TJ7on
.;  
!********************************************************************************************************* Svn|vH  
?yddr`?W  
    do i=Imin,Imax-1 gj\r>~S  
        do j=Jmin,Jmax KJ,{w?p~ )  
            Hy(i,j)=CP*Hy(i,j)+CQ*(Ez(i+1,j)-Ez(i,j))/dx                         !（见课本2-3-11式）求Hy '1ff|c!x9  
        end do h0Acpd2  
    end do L':;Vv~-  
EiI3$y3;  
s['F?GWg  
    !以上三式为二维情况FDTD电磁场的时间推进计算公式 En+4@BC  
*kgbcUf8  
    Ez(100,100)=sin(Omega*m*dt)                                                !在（100，100）处加入正弦点波源 ?v#t{e0eQ  
!************************************************************************************************************ _]1dm)%  
!                                                二阶Mur吸收边界                                                             n5 @H  
!************************************************************************************************************ hnL"f[p@gC  
xZtA
) Bp  
-`]B4Nt6  
!************************************************************************************************************ j9%u&  
!                                                y=0与y=200边界 Ts0.Ck  
!************************************************************************************************************ $J[h(>-X  
    do i=Imin+1,Imax-1     :g'"*VXYB  
        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)) ?)Z~H,Q(z  
        Eab(i,1)=Ez(i,Jmin+1) )8ctNpQt  
        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)) |k)h' ?  
        Eab(i,3)=Ez(i,Jmax-1) + }(  
    end do ,nMLua\  
!************************************************************************************************************ )<DL'  
!                                                x=0与x=200边界 =w".B[r  
!************************************************************************************************************ Xo(K*eIN  
&AGV0{NMh]  
    do j=Jmin+1,Jmax-1 RK/SeS  
        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)) : i3-7k  
        Eab(j,4)=Ez(1,j) gTW(2?xYf  
        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)) T9{94Ra  
        Eab(j,2)=Ez(Imax-1,j) vb9G_Pfz  
         $b\Gl=YX^  
    end do {#=q[jVi%1  
!************************************************************************************************************* -#3B>VY  
!                                                吸收边界角点处                 Mz40([{  
!************************************************************************************************************* ^ft_1d[  
?OYu BZF  
    Ez(Imin,Jmin)=Ez(Imin+1,Jmin+1)+C1*(Ez(Imin+1,Jmin+1)-Ez(Imin,Jmin)) vrGRZa  
    Ez(Imax,Jmin)=Ez(Imax-1,Jmin+1)+C1*(Ez(Imax-1,Jmin+1)-Ez(Imax,Jmin)) fcr\XCG7U  
    Ez(Imax,Jmax)=Ez(Imax-1,Jmax-1)+C1*(Ez(Imax-1,Jmax-1)-Ez(Imax,Jmax)) SWO$#X /  
    Ez(Imax,Jmax)=Ez(Imin+1,Jmax-1)+C1*(Ez(Imin+1,Jmax-1))-Ez(Imax,Jmax) M^mS#<!y  
     Cf<i"   
!********************************************************************************************************* S`K8e^]
 
!                                                多文件输出 ~8`r.1aUO  
!********************************************************************************************************* RNcHU  
100   format(1350f16.5) 1P3^il7  
    if(m==300)THEN & @^|=>L  
         pb$U~TvzhM  
      OPEN(1,File='Ez300.dat') %l,p />r  
0mH>fs 4  
        do j=Jmin,Jmax q3T'rw%Eh  
          write(1,100) (Ez(i,j),i=Imin,Imax) H1 n`A#6?  
        end do nW{). P  
vNd4Fn)H  
      close(1) E$4\Yc)(AL  
,S:g5n>M  
~(-B%Az  
    elseif(m==500)then w80g)4V+  
|6"zIHvtc  
        OPEN(1,File='Ez500.dat') 0#G&8*FMN  
q,^^c
1f  
        do j=Jmin,Jmax ;,JCA# N  
          write(1,100) (Ez(i,j),i=Imin,Imax) 8>T '  
        end do B:&/*HU  
        close(1) @wYuc
{%S  
d.% Vm&3  
\.9-:\'(  
    elseif(m==1000)then ;l &mA1+  
        OPEN(1,File='Ez1000.dat')
PZf^r  
lk%rE  
        do j=Jmin,Jmax F
,P,dc  
          write(1,100) (Ez(i,j),i=Imin,Imax) $vTU|o>|  
        end do `+\+  
      close(1) B<H5WI  
%S]5wR6;_  
|mc!v*O  
    elseif(m==3000)then + 3aAL&  
        OPEN(1,File='Ez3000.dat') 1 BAnf9  
 Sl  
        do j=Jmin,Jmax S3P;@Rm  
          write(1,100) (Ez(i,j),i=Imin,Imax) "So+  
        end do Cc]t*;nU_  
      close(1) {E@Vh  
iU~oPp[e  
+smPR  
SE;Jl[PgcL  
    elseif(m==5000)then pI( OI>~3  
        OPEN(1,File='Ez5000.dat') mmu{K$9}I  
|bO}|X  
        do j=Jmin,Jmax RoJ{ ou@cs  
          write(1,100) (Ez(i,j),i=Imin,Imax) }Rt?p8p  
        end do =eDVgOZ)  
      close(1) :jT1=PfL  
Hb#8?{  
wg<DV!GZ  
    elseif(m==6000)then ]Yp;8#:1  
        OPEN(1,File='Ez6000.dat') H8mmmt6g  
mKvk6OC  
        do j=Jmin,Jmax 3*/y<Z'H  
          write(1,100) (Ez(i,j),i=Imin,Imax) tOn/r@Fd^E  
        end do u1~H1 ]Ii  
      close(1) <omSK- T
-  
}(hx$G^M  
0AZ Vc  
      elseif(m==7000)then GghZ".O  
        OPEN(1,File='Ez7000.dat') nkG1&wiX  
jRmv~]  
        do j=Jmin,Jmax ~Z=Q+'Hu0  
          write(1,100) (Ez(i,j),i=Imin,Imax) 2h@/Q)z  
        end do >j4;{r+eQw  
      close(1) q$yTG!q*  
_<+!  
N41)?-7F  
&j<B22t!  
    elseif(m==8000)then Xa
t>d>nJ]  
        OPEN(1,File='Ez8000.dat') kOfbO'O9  
~dkS-6q~Q  
        do j=Jmin,Jmax ;/XWX$G@  
          write(1,100) (Ez(i,j),i=Imin,Imax) MXh0a@*]  
        end do QgqR93Ic  
      close(1) 2TXrVaM  
P6!c-\  
H!y1&  
    endif *Ew`Fm H  
end do /vYuwaWG=  
mp*&{[XoVC  
end p`fUpARA!  
!*********************************************************************************************************** Sp]"Xr)  
!                                           end   program TTZb.  
!*********************************************************************************************************** }R=n!Y$F  
                 M2W4 RovfR  
     =2F;'T\6  
_D1Uc|  
)\s{\u \  

能不能把你的这个程序给贴出来呢？我写的这个程序在角点总反射，也不知道什么原因···谢谢了···

出差时间较长，好久没有上来，现在才看到，实在抱歉，请问还需要吗，需要的话我回去找找贴上来

还需要啊，楼主，我也做了一个三维的PML程序，但是四个角点的反射很大，也不知道是什么原因。

同楼上一样是三维PML出了问题 纠结了好久也找不到错误