[资料共享]发一个我自己在matlab下编写的高斯制fdtd+pml程序代码

        前段时间一直在论坛里面向大家求助有关于高斯制下完全匹配层的差分公式的问题，在大家的热心帮助下，现在我基本上在高斯制下实现了同国际制下同等的吸收效果，这里把我所用的控制方程和程序代码传上来回馈论坛中各位热心坛友的帮助！由于我的编程思路不是很规范，所以在吸收效果上与理论上的完全匹配层吸收效果还有差距，但是已经基本达到了我的预期目的，我在这里也附上几张效果图。目前存在的问题是角点的反射比较严重，也恳请有熟悉这方面的编程高手帮我看看，该怎么改动才能取掉角点的反射，实现更好的吸收。请多多指教！！ mwutv8?  
    程序代码如下： TVF:z_M9  
% 二维FDTD  TE波仿真 'cqY-64CJZ  
clear all; 26}fB  
% 定义常数 :_|Xr'n`A  
pi=3.1415; ! u@JH`  
c=3.0e10;                     %高斯制下光速 y@(EGfI  
f=1.0e15;                     %频率 wf_ $#.;m  
lambda=c/f;                   %波长 pJJO
y  
nmax=400;                     %时间步数 .%+anVXS  
del_s=lambda/20;              %每最小波长20个采样点           h_Ssm{C\  
del_t=0.5*del_s/c;            %迭代时间步长 'p[*2J"K4  
n=182;                         %真空区域网格数 8 ,W*)Q  
np=9;                         %pml层数 %\L{Ud%7  
N1=n+2*np;                    %总网格数 J"-/ok(<@  
N=N1+1;                       %采样点数 ITy/h]0  
M=4;                          %导电率渐变指数 XAN.Plk  
sigma_max=(M+1)/1.50/pi/del_s;%最大导电率 @RZbo@{~  
rZ&li/Z
 
% TE波的分量初始化 vdh[%T,&  
tic; @&x'.2[nv  
figure(1); DzIV5FG  
axis([0 N 0 N -0.5 0.5]); 2sj: &][R  
Ex=zeros(N1,N);               %x方向为横向，采样点为网格的横向边，故行数+1 JS/~6'uB  
Ey=zeros(N,N1);               %y方向为纵向，采样点为网格的纵向边，故列数+1 mQ9y{}t=4  
Bz=zeros(N1,N1);              %矩阵行为纵向网格数，矩阵列为横向网格数,循环中用j表示行数，i表示列数 Pk;1q?tGw  
Bzx=zeros(N1,N1); w"O{@2B3:H  
Bzy=zeros(N1,N1); FX!Qd&kl1  
Bzxx=zeros(nmax,2); {v"Y!/ [z  
A2!pbeG  
%进入电磁场迭代计算 
26p_fKY  
for tt=1:nmax <'48mip  
  for i=1:N1 a ?\:,5=  
        if i>=np+1&&i<=N1-np
D`fIw` _  
           di=0; i91k0q*di  
        elseif i<=np {e1sq^>|  
           di=np-i+0.5; l.(|&U~  
        elseif i>=N1-np+1 ,/qY 9eh  
           di=np+i-N1-0.5; BMtk/r/  
        end                                 %di是采样点横向距PML内边界的距离 Nza@6nI"  
         sigma_mx=sigma_max*(di/np)^M; :`>tCYy;  
     for j=1:N1   u2IU/z8 ^  
         if j>=np+1&&j<=N1-np \]g51U!'  
                dj=0;  @{Dfro  
         elseif j<=np _S,UpR~2W  
                dj=np-j+0.5; G#~U\QlG-  
         elseif j>=N1-np+1 \E'z+0  
                dj=np+j-N1-0.5; _U9.u#>sV  
         end                                 %dj是采样点纵向距PML内边界的距离 R8 LHwRQ  
            sigma_my=sigma_max*(dj/np)^M;                                     C~WWuju'  
         5jv*C]z  
              if sigma_mx+sigma_my==0        %真空区 !DUOi4I  
                  if j==100&&i==100 7uf5w0]
  
                      t=30;                %可选择高斯脉冲 `(/xj{"Fr}  
                      term=(tt-t); &fsk ESV0  
                      pulse=exp(-4*pi*term^2/20^2) ; 1p$*N  
%                       pulse=sin(2*pi*tt/40); %可选正弦时谐源 =?_:h`}  
                      c_miu=c*del_t/del_s; SW5n?Qj3-  
                      Eterm1=c_miu*(Ex(i,j+1)-Ex(i,j)); nvA7eTO6C  
                      Eterm2=c_miu*(Ey(i+1,j)-Ey(i,j)); ,4=mlte"  
                      Bz(i,j)=Bz(i,j)+Eterm1-Eterm2+pulse;%加入脉冲源 <rvM)EJv|  
                  else At'M? Q@v  
                      c_miu=c*del_t/del_s; P4LiU2C  
                      Eterm1=c_miu*(Ex(i,j+1)-Ex(i,j)); b`JS&E  
                      Eterm2=c_miu*(Ey(i+1,j)-Ey(i,j)); v4K! BW  
                      Bz(i,j)=Bz(i,j)+Eterm1-Eterm2; \}\# fg  
                  end v(6[z)A0  
              else                            %PML区 m.<or?l'y>  
                    cpm=(1-2*c*sigma_mx*del_t)/(1+2*c*sigma_mx*del_t); j{johV+`8  
                    cqm=c*del_t/(1+2*c*sigma_mx*del_t)/del_s; F8#MI G   
                    Bzx(i,j)=cpm*Bzx(i,j)-cqm*(Ey(i+1,j)-Ey(i,j));               / U!xh3  
                    cpm=(1-2*c*sigma_my*del_t)/(1+2*c*sigma_my*del_t); RN}joKV  
                    cqm=c*del_t/(1+2*c*sigma_my*del_t)/del_s; B)/c]"@89  
                    Bzy(i,j)=cpm*Bzy(i,j)+cqm*(Ex(i,j+1)-Ex(i,j)); C$$Zwgy  
                    Bz(i,j)=Bzx(i,j)+Bzy(i,j); S<jiy<|`  
              end \qTp#sF  
     end  6@"E*-z$  
  end Hy -)yR  
  for i=2:N1 EPA 2_  
      if i>=np+1&&i<=N1-np B%`|W@v  
         di=0; -0TI7 @  
      elseif i<=np H?FiZy*[Y  
         di=np-i+1; n/`!G?kvI  
      elseif i>=N1-np+1 %zo 6A1Q;  
         di=np+i-N1-1; 76]Z~^Y  
      end                                   %di是采样点横向距PML内边界的距离 O.Y|},F  
      sigma_ex=sigma_max*(di/np)^M; v V;]?  
     for j=1:N1 XB@i
{/6K  
         cam=(1-2*c*sigma_ex*del_t)/(1+2*c*sigma_ex*del_t); A'Z!l20_  
         cbm=c*del_t/(1+2*c*sigma_ex*del_t)/del_s; ko|M2\  
         Ey(i,j)=cam*Ey(i,j)-cbm*(Bz(i,j)-Bz(i-1,j)); WYP\J1sy  
     end "a(e2H2&T4  
  end (N/-blto  
  for i=1:N1 2b
G92  
      for j=2:N1 :S}ZF$ $j%  
          if j>=np+1&&j<=N1-np q(csZ\e=  
             dj=0; &g>MZ"Z|  
          elseif j<=np 4"\x#  
             dj=np-j+1; @BPQ >  
          elseif j>=N1-np+1 C(UWir3mW?  
             dj=np+j-N1-1; Vk2%yw>  
          end                               %dj是采样点纵向距PML内边界的距离 L,!\PV|  
          sigma_ey=sigma_max*(dj/np)^M; Dm$SW<!l|  
          cam=(1-2*c*sigma_ey*del_t)/(1+2*c*sigma_ey*del_t); #DARZhU)  
          cbm=c*del_t/(1+2*c*sigma_ey*del_t)/del_s;< .. 7HQL^Q  
I8 [ *  

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

角点反射比较严重，这个情况在国际单位制和高斯制中都存在吗？

       嗯，我编的程序国际制下角点的反射也比较严重，可能是我的分区方法不对的缘故，我看看换个思路再尝试一下吧。 ;:w0%>X^  
   前段时间我采用的差分格式和现在差不多，但是一直实现不了较好的吸收效果，原因就是出在sigma_max的设置上面，在高斯制下长度单位变为cm，所以同样的空间步长在高斯制下数量级就扩大了100倍，而sigma_max的计算式里面分母含有空间步长，所以如果不加以修正，sigma_max在高斯制下就变成了国际制下的1%，太大的导电率和太小的导电率都不能实现较好的吸收效果，现在我把它修正后就差不多了

呵呵，和我最开始预想的一样，长度单位还是有影响的。 ExZ|_7^<  
角点的反射率大概有多少？1%?

那个我还没有进行量级的分析，看图就可以看得出来，怕是不只1%吧

很好，多谢分享！研究一下！

留个脚印，编程的时候参考

图形很好看啊，我的MRTD什么时候能达到这个程度就好了！

坚持+努力 ，加油！ _Ud!tK*H  

引用第7楼ghost32于2008-12-11 21:42发表的  : ?1r<`o3l\  
图形很好看啊，我的MRTD什么时候能达到这个程度就好了！

角点反射比较严重，这个情况在国际单位制和高斯制中都存在吗？ lXEnm-_  
您正在看的文章来自RFEDA微波社区--微波仿真论坛 原文地址:http://bbs.rfeda.cn/read.php?tid=17082