fdtd_2d_te.m

二维高斯制式的TE波基于FDFT方法的Matlab仿真程序！

% 二维FDTD  TE波仿真 
clear all;
% 定义常数
pi=3.1415;
c=3.0e10;                     %高斯制下光速
f=1.0e15;                     %频率
lambda=c/f;                   %波长
nmax=400;                     %时间步数
del_s=lambda/20;              %每最小波长20个采样点          
del_t=0.5*del_s/c;            %迭代时间步长
n=182;                         %真空区域网格数
np=9;                         %pml层数
N1=n+2*np;                    %总网格数
N=N1+1;                       %采样点数
M=4;                          %导电率渐变指数
sigma_max=(M+1)/1.50/pi/del_s;%最大导电率

% TE波的分量初始化
tic;
figure(1);
axis([0 N 0 N -0.5 0.5]);
Ex=zeros(N1,N);               %x方向为横向，采样点为网格的横向边，故行数+1
Ey=zeros(N,N1);               %y方向为纵向，采样点为网格的纵向边，故列数+1
Bz=zeros(N1,N1);              %矩阵行为纵向网格数，矩阵列为横向网格数,循环中用j表示行数，i表示列数
Bzx=zeros(N1,N1);
Bzy=zeros(N1,N1);
Bzxx=zeros(nmax,2);

%进入电磁场迭代计算
for tt=1:nmax
  for i=1:N1
        if i>=np+1&&i<=N1-np
           di=0;
        elseif i<=np
           di=np-i+0.5;
        elseif i>=N1-np+1
           di=np+i-N1-0.5;
        end                                 %di是采样点横向距PML内边界的距离
         sigma_mx=sigma_max*(di/np)^M;
     for j=1:N1  
         if j>=np+1&&j<=N1-np
                dj=0;
         elseif j<=np
                dj=np-j+0.5;
         elseif j>=N1-np+1
                dj=np+j-N1-0.5;
         end                                 %dj是采样点纵向距PML内边界的距离
            sigma_my=sigma_max*(dj/np)^M;                                     
         
              if sigma_mx+sigma_my==0        %真空区
                  if j==100&&i==100
                      t=30;                %可选择高斯脉冲
                      term=(tt-t);
                      pulse=exp(-4*pi*term^2/20^2) ;
%                       pulse=sin(2*pi*tt/40); %可选正弦时谐源
                      c_miu=c*del_t/del_s;
                      Eterm1=c_miu*(Ex(i,j+1)-Ex(i,j));
                      Eterm2=c_miu*(Ey(i+1,j)-Ey(i,j));
                      Bz(i,j)=Bz(i,j)+Eterm1-Eterm2+pulse;%加入脉冲源
                  else
                      c_miu=c*del_t/del_s;
                      Eterm1=c_miu*(Ex(i,j+1)-Ex(i,j));
                      Eterm2=c_miu*(Ey(i+1,j)-Ey(i,j));
                      Bz(i,j)=Bz(i,j)+Eterm1-Eterm2;
                  end
              else                            %PML区
                    cpm=(1-2*c*sigma_mx*del_t)/(1+2*c*sigma_mx*del_t);
                    cqm=c*del_t/(1+2*c*sigma_mx*del_t)/del_s;
                    Bzx(i,j)=cpm*Bzx(i,j)-cqm*(Ey(i+1,j)-Ey(i,j));               
                    cpm=(1-2*c*sigma_my*del_t)/(1+2*c*sigma_my*del_t);
                    cqm=c*del_t/(1+2*c*sigma_my*del_t)/del_s;
                    Bzy(i,j)=cpm*Bzy(i,j)+cqm*(Ex(i,j+1)-Ex(i,j));
                    Bz(i,j)=Bzx(i,j)+Bzy(i,j);
              end
     end
  end
  for i=2:N1
      if i>=np+1&&i<=N1-np
         di=0;
      elseif i<=np
         di=np-i+1;
      elseif i>=N1-np+1
         di=np+i-N1-1;
      end                                   %di是采样点横向距PML内边界的距离
      sigma_ex=sigma_max*(di/np)^M;
     for j=1:N1
         cam=(1-2*c*sigma_ex*del_t)/(1+2*c*sigma_ex*del_t);
         cbm=c*del_t/(1+2*c*sigma_ex*del_t)/del_s;
         Ey(i,j)=cam*Ey(i,j)-cbm*(Bz(i,j)-Bz(i-1,j));
     end
  end
  for i=1:N1
      for j=2:N1 
          if j>=np+1&&j<=N1-np
             dj=0;
          elseif j<=np
             dj=np-j+1;
          elseif j>=N1-np+1
             dj=np+j-N1-1;
          end                               %dj是采样点纵向距PML内边界的距离
          sigma_ey=sigma_max*(dj/np)^M;
          cam=(1-2*c*sigma_ey*del_t)/(1+2*c*sigma_ey*del_t);
          cbm=c*del_t/(1+2*c*sigma_ey*del_t)/del_s;
          Ex(i,j)=cam*Ex(i,j)+cbm*(Bz(i,j)-Bz(i,j-1));   
      end 
  end    
  Bzxx(tt,1)=Bz(90,50);                     %对靠近边界的中央磁场点采样
  Bzxx(tt,2)=Bz(50,90);
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%电磁场的计算部分结束
  figure(1);                                %可视化处理
  clf;   
  mesh(Bz);                                 %磁场的幅值
  axis([0 N 0 N -0.5 0.5]);
  xlabel('i')
  ylabel('j')
  drawnow;
end
figure(2);
plot(Bzxx);
figure(3);
title('磁场幅值分布图');
surface(Bz);
shading interp;
axis square;
toc