upml_2d_tmmode.m

二维TM模电磁波FDTD仿真

%top
%Ez
for ii=1:2*iebc+ix
     for jj=iebc+jy+1:iebc+jy+iebc
         eff=M_epsr(ii,jj);   %未平均化的eps
%          Omax=(m+1)/sqrt(eff)/150/pi/dx;
         if (ii<=iebc)
           
             x1  =(iebc-ii)*dx;
             x2  =(iebc-ii+1)*dx;
             Ox  =Omax*Consx*(x2^(m+1)-x1^(m+1));   %半网格处的值
         else
             if ii>iebc+ix
                 x1=(ii-iebc-ix)*dx;
                 x2=(ii-iebc-ix-1)*dx;
                 Ox  =Omax*Consx*(x1^(m+1)-x2^(m+1)); 
             else
                 Ox=0;
             end
         end
%          Omax=(m+1)/sqrt(eff)/150/pi/dy;
         y1=(jj-iebc-jy)*dy;
         y2=(jj-iebc-jy-1)*dy;
         Oy =Omax*Consy*(y1^(m+1)-y2^(m+1));   %半网格处的值
         DCP(ii,jj)=(2*e-Oy*dt)/(2*e+Oy*dt);
         DCQ(ii,jj)=2*e*dt/(2*e+Oy*dt);
         ECP(ii,jj)=(2*e-Ox*dt)/(2*e+Ox*dt);
         ECQ(ii,jj)=2*e/(2*e+Ox*dt)/eff/e;
     end
 end
 
 %Hy top
 %对左右边界上的Hy系数赋0
 for ii=2:2*iebc+ix
     for jj=iebc+jy+1:2*iebc+jy
         eff=(M_epsr(ii,jj)+M_epsr(ii-1,jj))/2;     %应用X方向平均后的eps
         y1=(jj-iebc-jy)*dy;
         y2=(jj-iebc-jy-1)*dy;
         Oy=Omax*Consy*(y1^(m+1)-y2^(m+1));     %应用半网格处的Oy
         %分为三个区域并对交接处特殊处理Ox
%          Omax=(m+1)/sqrt(eff)/150/pi/dx;
         if ii==iebc+1
         Ox  = Omax*Consx*(0.5*dx)^(m+1);     %边界上特殊处理的Ox
         else
             if ii==iebc+ix+1
             Ox=Omax*Consx*(0.5*dx)^(m+1);       %边界上特殊处理的Ox
             else
                 if  ii<iebc+1
                     x1=(iebc-ii+0.5)*dx;
                     x2=(iebc-ii+1.5)*dx;
                     Ox= Omax*Consx*(x2^(m+1)-x1^(m+1));                       %应用整网格上的Ox
                 else
                     if ii>iebc+ix+1
                         x1=(ii-iebc-ix-1.5)*dx;
                         x2=(ii-iebc-ix-0.5)*dx;
                         Ox=Omax*Consx*(x2^(m+1)-x1^(m+1));                 %应用整网格上的Ox
                     else
                         Ox=0;
                     end
                 end
             end
         end
         ByCA(ii,jj)=1;
         ByCB(ii,jj)=dt;
         HyCA(ii,jj)=(2*e-Ox*dt)/(2*e+Ox*dt);
         HyCB1(ii,jj)=(2*e+Oy*dt)/(2*e+Ox*dt)/muz;
         HyCB2(ii,jj)=(2*e-Oy*dt)/(2*e+Ox*dt)/muz;
     end
 end
 
 %Hx top
 %对两个角区域的Ox特殊处理,同时对边界处的Oy特殊处理
 for ii=1:2*iebc+ix
     for jj=iebc+jy+1:iebc+jy+iebc
         eff=(M_epsr(ii,jj)+M_epsr(ii,jj-1))/2;    %在y方向平均的eps
%          Omax=(m+1)/sqrt(eff)/150/pi/dy;
     if jj==iebc+jy+1
         Oy=Omax*Consy*(0.5*dy)^(m+1);       % 边界上的Oy特殊处理
     else
         y1=dy*(jj-iebc-jy-0.5);
         y2=dy*(jj-iebc-jy-1.5);
         Oy=Omax*Consy*(y1^(m+1)-y2^(m+1));      %整网格处的Oy
     end
     
%      Omax=(m+1)/sqrt(eff)/150/pi/dx;
     if ii<=iebc
         x1=(iebc-ii)*dx;
         x2=(iebc-ii+1)*dx;
         Ox=Omax*Consx*(x2^(m+1)-x1^(m+1));      %半网格处的 Ox
     else
         if ii>iebc+ix
             x1=(ii-iebc-ix)*dx;
             x2=(ii-iebc-ix-1)*dx;
             Ox=Omax*Consx*(x1^(m+1)-x2^(m+1));   %半网格处的Ox
         else
             Ox=0;
         end
     end
     BxCA(ii,jj)=1;
     BxCB(ii,jj)=dt;
     HxCA(ii,jj)=(2*e-Oy*dt)/(2*e+Oy*dt);
     HxCB1(ii,jj)=(2*e+Ox*dt)/(2*e+Oy*dt)/muz;
     HxCB2(ii,jj)=(2*e-Ox*dt)/(2*e+Oy*dt)/muz;
     end
  end
 
  %计算区域的系数处理
  %Ez CS
  for ii=iebc+1:iebc+ix
      for jj=iebc+1:iebc+jy
          Ox=0;
          Oy=0;
          eff=M_epsr(ii,jj);
          DCP(ii,jj)=1;
          DCQ(ii,jj)=dt;
          ECP(ii,jj)=1;
          ECQ(ii,jj)=1/eff/e;
      end
  end
  
  %Hy CS
  for ii=iebc+1:iebc+ix
      for jj=iebc+1:iebc+jy
          eff=(M_epsr(ii,jj)+M_epsr(ii-1,jj))/2;
          if ii==iebc+1
%               Omax=(m+1)/sqrt(eff)/150/pi/dx;
              Ox=Omax*Consx*(0.5*dx)^(m+1);
          else
              Ox=0;
          end
          Oy=0;
          ByCA(ii,jj)=1;
          ByCB(ii,jj)=dt;
          HyCA(ii,jj)=(2*e-Ox*dt)/(2*e+Ox*dt);
          HyCB1(ii,jj)=(2*e+Oy*dt)/(2*e+Ox*dt)/muz;
          HyCB2(ii,jj)=(2*e-Oy*dt)/(2*e+Ox*dt)/muz;
      end
  end
  
  %Hx CS
  for ii=iebc+1:iebc+ix
      for jj=iebc+1:iebc+jy
          eff=(M_epsr(ii,jj)+M_epsr(ii,jj-1))/2;
          if jj==iebc+1 % <--+++++++++++++++++++++++ if ii==iebc+1
%               Omax=(m+1)/sqrt(eff)/150/pi/dy;
              Oy=Omax*Consx*(0.5*dx)^(m+1);
          else
              Oy=0;
          end
          Ox=0;
          BxCA(ii,jj)=1;
          BxCB(ii,jj)=dt;
          HxCA(ii,jj)=(2*e-Oy*dt)/(2*e+Oy*dt);
          HxCB1(ii,jj)=(2*e+Ox*dt)/(2*e+Oy*dt)/muz;
          HxCB2(ii,jj)=(2*e-Ox*dt)/(2*e+Oy*dt)/muz;
      end
  end

  %Movie initiate            
   %*********************************************************************** 
%     Movie initialization 
%*********************************************************************** 

% subplot(3,1,1),pcolor(Hx((iebc+1):(iebc+ix+1),(iebc+1):(iebc+1+jy))'); 
% shading flat; 
% caxis([-80.0 80.0]); 
% axis([1 ie 1 jb]); 
% colorbar; 
% axis image; 
% axis off; 
% title(['Hx at time step = 0']); 
% 
% subplot(3,1,2),pcolor(Hy((iebc+1):(iebc+ix+1),(iebc+1):(iebc+1+jy))'); 
% shading flat; 
% caxis([-80.0 80.0]); 
% axis([1 ib 1 je]); 
% colorbar; 
% axis image; 
% axis off; 
% title(['Hy at time step = 0']); 
% 
% subplot(3,1,3),pcolor(Ez((iebc+1):(iebc+ix+1),(iebc+1):(iebc+1+jy))'); 
% shading flat; 
% caxis([-0.2 0.2]); 
% axis([1 ie 1 je]); 
% colorbar; 
% axis image; 
% axis off; 
% title(['Ez at time step = 0']); 
% 
% rect=get(gcf,'Position'); 
% rect(1:2)=[0 0]; 
% 
% M=moviein(nmax/4,gcf,rect);        

%ramp CW source implement parameters


  for n=1:nmax
      dstore(1:ie,1:je)=Dz(1:ie,1:je);
      Dz(1:ie,1:je)=DCP(1:ie,1:je).*Dz(1:ie,1:je)-DCQ(1:ie,1:je).*...
           ((Hx(1:ie,2:jb)-Hx(1:ie,1:je))  / dy -...
             (Hy(2:ib,1:je)-Hy(1:ie,1:je)) / dx );
      Ez(1:ie,1:je) = ECP(1:ie,1:je).* Ez(1:ie,1:je) + ECQ(1:ie,1:je) .*...
           (Dz(1:ie,1:je) -  dstore(1:ie,1:je));
%       Hz(iebc+40,20)=source(n);
%       source=sin(omega*(n)*dt); 

      start_x=5;start_y=5;length=10;
      sxsta=iebc+start_x; sxend=iebc+start_x+length-1; systa=iebc+start_y+5;
      Ez(iebc+start_x:iebc+start_x+length-1,iebc+start_y+5)=source(n);
%       Ez(iebc+start_x+10,iebc+start_y:iebc+start_y+width)=source(n);
%       Hz(iebc+start_x:iebc+start_x+length,iebc+start_y+5)=
% Hz(79+iebc:101+iebc,iebc+1)=source(n)
      dstore(1:ib,1:je)=By(1:ib,1:je);
      By(2:ie,1:je)=ByCA(2:ie,1:je).* By(2:ie,1:je) + ByCB(2:ie,1:je) .*... 
           ((Ez(2:ie,1:je)-Ez(1:(ie-1),1:je)) / dx);
      Hy(2:ie,1:je)= HyCA(2:ie,1:je) .* Hy(2:ie,1:je) + (HyCB1(2:ie,1:je) .* By(2:ie,1:je) - HyCB2(2:ie,1:je) .* dstore(2:ie,1:je));
      
      dstore(1:ib,1:jb)=Bx(1:ib,1:jb);
      Bx(1:ie,2:je) = BxCA(1:ie,2:je).* Bx(1:ie,2:je) - BxCB(1:ie,2:je) .*...
           ((Ez(1:ie,2:je)-Ez(1:ie,1:(je-1))) / dy);
      Hx(1:ie,2:je) =HxCA(1:ie,2:je) .* Hx(1:ie,2:je) + (HxCB1(1:ie,2:je) .* Bx(1:ie,2:je) - HxCB2(1:ie,2:je) .* dstore(1:ie,2:je));

%       probe(n)= Ez(iebc+start_x+65,iebc+start_y+width/2);
%       probe(n)= Ez(iebc+start_x+5,iebc+start_y+width-40);
%*****************plot field distribution of E/H
%*********************************************************************** 
%     Visualize fields 
%*********************************************************************** 
% if mod(n,100)==0
%     save n Ez;
% end
% if mod(n,4)==0; 
% 
% timestep=int2str(n); 
% % surf(Hz');
% subplot(3,1,1),pcolor(Hx((iebc+1):(iebc+ix+1),(iebc+1):(iebc+1+jy))'); 
% % subplot(3,1,1),pcolor(Ex'); 
% shading flat; 
% caxis([-80.0 80.0]); 
% axis([1 ib 1 jb]); 
% colorbar; 
% axis image; 
% axis off; 
% title(['Hx at time step = ',timestep]); 
% 
% subplot(3,1,2),pcolor(Hy((iebc+1):(iebc+ix+1),(iebc+1):(iebc+1+jy))'); 
% % subplot(3,1,2),pcolor(Ey'); 
% shading flat; 
% caxis([-80.0 80.0]); 
% axis([1 ib 1 jb]); 
% colorbar; 
% axis image; 
% axis off; 
% title(['Hy at time step = ',timestep]); 
% % 
% subplot(3,1,3),pcolor(Ez((iebc+1):(iebc+ix+1),(iebc+1):(iebc+1+jy))'); 
% % subplot(3,1,3),pcolor(Hz'); 
% shading flat; 
% caxis([-0.2 0.2]); 
% axis([1 ib 1 jb]); 
% colorbar; 
% axis image; 
% axis off; 
% title(['Ez at time step = ',timestep]); 
% 
% nn=n/4; 
% M(:,nn)=getframe(gcf,rect); 
% % 
%  end
 probe1(n)=Ez(sxsta,systa);
 probe2(n)=Ez(sxsta+10,systa+10);
 
 
 %*****************other signal processing way
end

 clear n;
 n=1:1:nmax;
 figure(1);
 plot(n,probe1);
 figure(2);
 plot(n,probe2);