www.pudn.com > fdtd_planewave.rar > fdtd_2d_TM_PML_planewave.m

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