fdtd_2d_tm_pml_planewave.m

二维fdtd程序

clc
clear

N=100;
L=2000;

c=3*10^8;
f=1.5*10^9;
lamda=c/f;

ddx=lamda/20;
dt=ddx/(2*c);

epsz=1/(4*pi*9*10^9);
epsilon=1; 
sigma=0;  

spread=12;             
t0=40;                
i_source=N/2;                
j_source=N/2;

ez_inc=zeros(1,N+1);
hx_inc=zeros(1,N);
low1=0;
low2=0;
high1=0;
high2=0;

% 总场边界
ia=N/4;
ib=3*N/4;
ja=ia;
jb=ib;

npml=15;
halfgrid=20;

dz=zeros(N+1,N+1);
ez=zeros(N+1,N+1);
iz=zeros(N+1,N+1);
hx=zeros(N+1,N);
hy=zeros(N,N+1);
ihx=zeros(N+1,N);
ihy=zeros(N,N+1);


% PML 初始设置

for i=1:N+1
    gi2(i)=1;
    gi3(i)=1;
    fi1(i)=0;
end

for i=1:N
    fi2(i)=1;
    fi3(i)=1;
end

for j=1:N+1
    gj2(j)=1;
    gj3(j)=1;
    fj1(j)=0;
end

for j=1:N
    fj2(j)=1;
    fj3(j)=1;
end

for i=1:N                                  
    for j=1:N;
        ga(i,j)=1/(epsilon+sigma*dt/epsz);  
        gb(i,j)=sigma*dt/epsz;              
    end;
end;

% PML 阻抗渐变设置
for i=1:npml+1
    xnum=npml+1-i;
    xxn=xnum/npml;
    xn=0.33*(xxn^3);
    
    gi2(i)=1/(1+xn);
    gi2(N+2-i)=1/(1+xn);
    gi3(i)=(1-xn)/(1+xn);
    gi3(N+2-i)=(1-xn)/(1+xn);
    
    xxn=(xnum-0.5)/npml;
    xn=0.25*(xxn^3);
    fi1(i)=xn;
    fi1(N+2-i)=xn;
end
    
 for i=1:npml
     xnum=npml+1-i;
     xxn=(xnum-0.5)/npml;
     xn=0.25*(xxn^3);
     fi2(i)=1/(1+xn);
     fi2(N+1-i)=1/(1+xn);
     fi3(i)=(1-xn)/(1+xn);
     fi3(N+1-i)=(1-xn)/(1+xn);
 end
 
 for j=1:npml+1
    xnum=npml+1-j;
    xxn=xnum/npml;
    xn=0.33*(xxn^3);
    
    gj2(j)=1/(1+xn);
    gj2(N+2-j)=1/(1+xn);
    gj3(j)=(1-xn)/(1+xn);
    gj3(N+2-j)=(1-xn)/(1+xn);
    
    xxn=(xnum-0.5)/npml;
    xn=0.25*(xxn^3);
    fj1(j)=xn;
    fj1(N+2-j)=xn;
end
    
 for j=1:npml
     xnum=npml+1-j;
     xxn=(xnum-0.5)/npml;
     xn=0.25*(xxn^3);
     fj2(j)=1/(1+xn);
     fj2(N+1-j)=1/(1+xn);
     fj3(j)=(1-xn)/(1+xn);
     fj3(N+1-j)=(1-xn)/(1+xn);
 end 
 % 程序主体
 for t=1:L
   %  TM波Y方向传播
    for j=2:N
         ez_inc(j)=ez_inc(j)+0.5*(hx_inc(j-1)-hx_inc(j));
    end
    ez_inc(1)=low2;
    low2=low1;
    low1=ez_inc(2);
    
    ez_inc(N+1)=high2;
    high2=high1;
    high1=ez_inc(N);
    for i=2:N
         for j=2:N
            dz(i,j)=gi3(i)*gj3(j)*dz(i,j)+gi2(i)*gj2(j)*0.5*( hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1) );
         end
      end      
      % 源的设置
      pulse=exp((-0.5)*( (t0-t)/spread ).^2);
     % pulse=sin(2*pi*f*t*dt);
      %dz(i_source,j_source)=dz(i_source,j_source)+pulse;
      ez_inc(5)=pulse;
      
      for i=ia:ib
          dz(i,ja)=dz(i,ja)+0.5*hx_inc(ja-1);
          dz(i,jb)=dz(i,jb)-0.5*hx_inc(jb);
      end     
    %  Ez的迭代  
      for i=2:N
          for j=2:N
              ez(i,j)=ga(i,j)*( dz(i,j)-iz(i,j) );
              iz(i,j)=iz(i,j)+gb(i,j)*ez(i,j) ;
          end
      end      
     % 截断边界
      for j=1:N+1
          ez(1,j)=0;
          ez(N+1,j)=0;
      end      
      for i=1:N+1
          ez(i,1)=0;
          ez(i,N+1)=0;
      end      
      %  金属边界条件
     for i=N/2-halfgrid:N/2+halfgrid-1;
         for j=N/2-halfgrid:N/2+halfgrid-1;
           ez(i,j)=0;
       end
     end     
     for j=1:N
         hx_inc(j)=hx_inc(j)+0.5*(ez_inc(j)-ez_inc(j+1));
     end       
      % Hx的迭代
    for i=1:N+1
          for j=1:N
              curl_e=ez(i,j)-ez(i,j+1);
              ihx(i,j)=ihx(i,j)+fi1(i)*curl_e;
              hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*0.5*(curl_e+ihx(i,j));
         end
    end      
     for i=ia:ib
         hx(i,ja-1)=hx(i,ja-1)+0.5*ez_inc(ja);
         hx(i,jb)=hx(i,jb)-0.5*ez_inc(jb);
     end  
     % Hy的迭代
      for i=1:N    
        for j=1:N+1
            curl_e=ez(i+1,j)-ez(i,j);
            ihy(i,j)=ihy(i,j)+fj1(j)*curl_e;
            hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5*(curl_e+ihy(i,j));
        end
    end
    
    for j=ja:jb
        hy(ia-1,j)=hy(ia-1,j)-0.5*ez_inc(j);
        hy(ib,j)=hy(ib,j)+0.5*ez_inc(j);
    end
    % 结果显示
    mesh(ez)
    %view(90,90)
    tm=['T=',num2str(t)]
    text(10,100,0.5,tm)
    axis([1 101 1 101 -1 1]);
    drawnow;
end