d2pmlxishou.m

分裂场PML吸收层

clc,clear
%%%% 二维FDTD TE wave pml吸收边界
ke=80;%元包数
kp=8; %pml 层数
ex(ke+2*kp+1,ke+2*kp+1)=0;%实际电场217
ey(ke+2*kp+1,ke+2*kp+1)=0;%实际电场217
hz(ke+2*kp+1,ke+2*kp+1)=0;%实际磁场216，即两端为电场边界
hzx(ke+2*kp+1,ke+2*kp+1)=0;
hzy(ke+2*kp+1,ke+2*kp+1)=0;
kc=(ke+2*kp)/2;%激励源中心 hz(kc,kc)

e0=8.85*10^-12; %   真空介电常数
mu0=4*3.14159*10^-7; % 真空磁导率
c0=1/sqrt(e0*mu0);% 真空中光速
f=1e9; %  频率
lambda=c0/f;  %波长
dx=lambda/10;% dx 
dy=lambda/10;% dy 
dt=(dx/2)/c0; %dt

sigxm=0.33;  % 电导率最大值   
sigmxm=sigxm*(mu0/e0);%导磁率最大值 
nn=3;    %pml层中电导率.导磁率指数 

%%%%%%%% 开始主程序 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=0;
%Nsteps=input('Nsteps=')%循环次数
Nsteps=1000;
for n=1:Nsteps
    T=T+1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hz=hzx+hzy;
hz(kc,kc)=sin(2*pi*f*T*dt);

%%%%%%%%%%%%%%%%%%%%%%%%计算电场%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%%%%%%%%%%% ex%%%%%%%%%%%%%%%%%%
    for ki=1:(ke+2*kp)
        for kj=(kp+1):(ke+kp+1)
        ex(ki,kj)=ex(ki,kj)+0.5*(hz(ki,kj)-hz(ki,kj-1));
        end
    end
  %%%%%% 下边pml层电场ex%%%%%%%%%%%%%%%%%%%%%
 
  for ki=1:(ke+2*kp)
      for kj=2:kp
        sigx=sigxm*((kp+1-kj)/kp)^nn; 
        ex(ki,kj)=ex(ki,kj)*exp(-sigx*dt/e0)+(1-exp(-sigx*dt/e0))/(dy*sigx)*sqrt(e0/mu0)*(hz(ki,kj)-hz(ki,kj-1));
        end
  end
 
 %%%%%% 上边pml层电场ex%%%%%%%%%%%%%%%%%%%%%
    for ki=1:(ke+2*kp)
        for kj=(ke+kp+2):(ke+2*kp)
        sigx=sigxm*((kj-ke-kp-1)/kp)^nn;
        ex(ki,kj)=ex(ki,kj)*exp(-sigx*dt/e0)+(1-exp(-sigx*dt/e0))/(dy*sigx)*sqrt(e0/mu0)*(hz(ki,kj)-hz(ki,kj-1));
        end
    end
    
%%%%%%%%%%%%%% ey%%%%%%%%%%%%%%%%%%
    for ki=(kp+1):(ke+kp+1)
        for kj=1:(ke+2*kp)
            ey(ki,kj)=ey(ki,kj)-0.5*(hz(ki,kj)-hz(ki-1,kj));
        end
    end
 %%%%%%%%%%%%%%%左边pml层中ey%%%%%%%%%%%%%%%%%%%
  for ki=2:kp
      for kj=1:(ke+2*kp)
        sigx=sigxm*((kp+1-ki)/kp)^nn; 
        ey(ki,kj)=ey(ki,kj)*exp(-sigx*dt/e0)-(1-exp(-sigx*dt/e0))/(dx*sigx)*sqrt(e0/mu0)*(hz(ki,kj)-hz(ki-1,kj));
        end
    end
     %%%%%% 右边pml层电场ey%%%%%%%%%%%%%%%%%%%%%
    for ki=(ke+kp+2):(ke+2*kp)
        for kj=1:(ke+2*kp)
        sigx=sigxm*((ki-ke-kp-1)/kp)^nn;
        ey(ki,kj)=ey(ki,kj)*exp(-sigx*dt/e0)-(1-exp(-sigx*dt/e0))/(dx*sigx)*sqrt(e0/mu0)*(hz(ki,kj)-hz(ki-1,kj));
        end
    end

 
%%%%%%%%%%%%%%%电场循环结束%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  

%%%%%%%%%%%%%%%磁场循环
 %%%%%%%%%%%%%%下边pml中hzy,含sigmy  
 for ki=1:(ke+2*kp)
    for kj=1:kp
      sigmx=sigmxm*((kp-kj+0.5)/kp)^nn;
      hzy(ki,kj)=hzy(ki,kj)*exp(-sigmx*dt/mu0)+(1-exp(-sigmx*dt/mu0))/(dx*sigmx)*sqrt(mu0/e0)*(ex(ki,kj+1)-ex(ki,kj)); 
    end
end
%%%%中间hzy(不含sigmy)   
 for ki=1:(ke+2*kp)
     for kj=(kp+1):(kp+ke)
         hzy(ki,kj)=hzy(ki,kj)+0.5*(ex(ki,kj+1)-ex(ki,kj));
     end
 end
 %%%%%%%%%%%%%%%pml磁场%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
%%%%%%%%%%%%%%上边pml中hzy,含sigmx        
 for ki=1:(ke+2*kp)
     for kj=(kp+ke+1):(ke+2*kp)
      sigmx=sigmxm*((kj-kp-ke-0.5)/kp)^nn;  
      hzy(ki,kj)=hzy(ki,kj)*exp(-sigmx*dt/mu0)+(1-exp(-sigmx*dt/mu0))/(dx*sigmx)*sqrt(mu0/e0)*(ex(ki,kj+1)-ex(ki,kj)); 
     end
 end
 

%%%%%%%%%%%%%%%左边pml中hzx,含sigmx

for ki=1:kp
    for kj=1:(ke+2*kp)
      sigmx=sigmxm*((kp-ki+0.5)/kp)^nn;
      hzx(ki,kj)=hzx(ki,kj)*exp(-sigmx*dt/mu0)-(1-exp(-sigmx*dt/mu0))/(dx*sigmx)*sqrt(mu0/e0)*(ey(ki+1,kj)-ey(ki,kj)); 
    end
end
%%%%中间hzx(不含sigmx) 
  for ki=(kp+1):(kp+ke)
     for kj=1:(ke+2*kp)
         hzx(ki,kj)=hzx(ki,kj)-0.5*(ey(ki+1,kj)-ey(ki,kj));
     end
  end  
%%%%%%%%%%%%%%右边pml中hzx,含sigmx        
 for ki=(kp+ke+1):(ke+2*kp)
     for kj=1:(ke+2*kp)
      sigmx=sigmxm*((ki-kp-ke-0.5)/kp)^nn;  
      hzx(ki,kj)=hzx(ki,kj)*exp(-sigmx*dt/mu0)-(1-exp(-sigmx*dt/mu0))/(dx*sigmx)*sqrt(mu0/e0)*(ey(ki+1,kj)-ey(ki,kj)); 
     end
 end
 
 


 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%plot(hz(:,kc))
%axis([1  ke+2*kp -1 1])
end
%pcolor(hz)
contour(hz,5)
axis([kp kp+ke kp kp+ke ],'square')


shading interp
%print -djpeg D2pml相位图
print -djpeg D2pml等高线图