fdtd3d.m

3-D FDTD code with PML absorbing boundary conditions

title(['real Hz(i=15,j,k), time step = 0']);
xlabel('j coordinate');
ylabel('k coordinate');


% mview(:,:)=real(hz(15,je/2,:));
% tview(:,:)=real(hzxbct(15,je/2,:)+hzybct(15,je/2,:));
% subplot(2,1,1),plot(mview),caxis([-1.0 1.0]);
% ylabel('hz in main')
% subplot(2,1,2),plot(tview),caxis([-1.0 1.0]);
% ylabel('hz in left')

rect=get(gcf,'Position');
rect(1:2)=[0 0];

M=moviein(nmax/2,gcf,rect);


%***********************************************************************
%     BEGIN TIME-STEPPING LOOP
%***********************************************************************

for n=1:nmax

%***********************************************************************
%     Update electric fields (Ex 、Ey、Ez) in main grid
%***********************************************************************

ex(:,2:je,2:ke)=caex(:,2:je,2:ke).*ex(:,2:je,2:ke)+...
           cbex(:,2:je,2:ke).*(hz(:,2:je,2:ke)-hz(:,1:je-1,2:ke)-...
                                    hy(:,2:je,2:ke)+hy(:,2:je,1:ke-1));

ey(2:ie,:,2:ke)=caey(2:ie,:,2:ke).*ey(2:ie,:,2:ke)+...
           cbey(2:ie,:,2:ke).*(hx(2:ie,:,2:ke)-hx(2:ie,:,1:ke-1)-...
                              hz(2:ie,:,2:ke)+hz(1:ie-1,:,2:ke));
       
ez(2:ie,2:je,:)=caez(2:ie,2:je,:).*ez(2:ie,2:je,:)+...
           cbez(2:ie,2:je,:).*(hy(2:ie,2:je,:)-hy(1:ie-1,2:je,:)-...
                               hx(2:ie,2:je,:)+hx(2:ie,1:je-1,:));
%***********************************************************************
%     Update Exy in PML regions
%***********************************************************************

%     FRONT

exybcf(:,2:jebc,2:kefbc)=caexybcf(:,2:jebc,2:kefbc).*exybcf(:,2:jebc,2:kefbc)+...  
  cbexybcf(:,2:jebc,2:kefbc).*(hzxbcf(:,2:jebc,2:kefbc)+hzybcf(:,2:jebc,2:kefbc)-...
                      hzxbcf(:,1:jebc-1,2:kefbc)-hzybcf(:,1:jebc-1,2:kefbc));
   
exybcl(:,1,2:kefbc)=caexybcl(:,1,2:kefbc).*exybcl(:,1,2:kefbc)+...
  cbexybcl(:,1,2:kefbc).*(hzxbcl(:,1,2:kefbc)+hzybcl(:,1,2:kefbc)-...
            hzxbcf(1:iebc,jebc,2:kefbc)-hzybcf(1:iebc,jebc,2:kefbc));                      
                       
ex(1:ie,1,1:kb)=caex(1:ie,1,1:kb).*ex(1:ie,1,1:kb)+...
  cbex(1:ie,1,1:kb).*(hz(1:ie,1,1:kb)-hzxbcf(ibbc:iebc+ie,jebc,kbbc:kebc+kb)-...
                                     hzybcf(ibbc:iebc+ie,jebc,kbbc:kebc+kb));
 
exybcr(:,1,2:kefbc)=caexybcr(:,1,2:kefbc).*exybcr(:,1,2:kefbc)+...
  cbexybcr(:,1,2:kefbc).*(hzxbcr(:,1,2:kefbc)+hzybcr(:,1,2:kefbc)-...
  hzxbcf(ibbc+ie:iefbc,jebc,2:kefbc)-hzybcf(ibbc+ie:iefbc,jebc,2:kefbc));                                          
    
exybcd(:,1,2:kebc)=caexybcd(:,1,2:kebc).*exybcd(:,1,2:kebc)+...
  cbexybcd(:,1,2:kebc).*(hzxbcd(:,1,2:kebc)+hzybcd(:,1,2:kebc)-...
   hzxbcf(ibbc:iebc+ie,jebc,2:kebc)-hzybcf(ibbc:iebc+ie,jebc,2:kebc));   
 
exybct(:,1,2:kebc)=caexybct(:,1,2:kebc).*exybct(:,1,2:kebc)+...
  cbexybct(:,1,2:kebc).*(hzxbct(:,1,2:kebc)+hzybct(:,1,2:kebc)-...
   hzxbcf(ibbc:iebc+ie,jebc,kbbc+kb:kefbc)-hzybcf(ibbc:iebc+ie,jebc,kbbc+kb:kefbc));    
   
                                     
%     BACK

exybcb(:,2:jebc,2:kefbc)=caexybcb(:,2:jebc,2:kefbc).*exybcb(:,2:jebc,2:kefbc)+...
  cbexybcb(:,2:jebc,2:kefbc).*(hzxbcb(:,2:jebc,2:kefbc)+hzybcb(:,2:jebc,2:kefbc)-...
                        hzxbcb(:,1:jebc-1,2:kefbc)-hzybcb(:,1:jebc-1,2:kefbc));
                    
exybcl(:,jb,2:kefbc)=caexybcl(:,jb,2:kefbc).*exybcl(:,jb,2:kefbc)+...
  cbexybcl(:,jb,2:kefbc).*(hzxbcb(1:iebc,1,2:kefbc)+hzybcb(1:iebc,1,2:kefbc)-...
                                       hzxbcl(:,je,2:kefbc)-hzybcl(:,je,2:kefbc));                        
                        
ex(1:ie,jb,1:kb)=caex(1:ie,jb,1:kb).*ex(1:ie,jb,1:kb)+...
  cbex(1:ie,jb,1:kb).*(hzxbcb(ibbc:iebc+ie,1,kbbc:kebc+kb)+...
                       hzybcb(ibbc:iebc+ie,1,kbbc:kebc+kb)-hz(1:ie,je,1:kb));

exybcr(:,jb,2:kefbc)=caexybcr(:,jb,2:kefbc).*exybcr(:,jb,2:kefbc)+...
  cbexybcr(:,jb,2:kefbc).*(hzxbcb(ibbc+ie:iefbc,1,2:kefbc)+...
  hzybcb(ibbc+ie:iefbc,1,2:kefbc)-hzxbcr(:,je,2:kefbc)-hzybcr(:,je,2:kefbc));                  
 
exybcd(:,jb,2:kebc)=caexybcd(:,jb,2:kebc).*exybcd(:,jb,2:kebc)+...
  cbexybcd(:,jb,2:kebc).*(hzxbcb(ibbc:iebc+ie,1,2:kebc)+...
  hzybcb(ibbc:iebc+ie,1,2:kebc)-hzxbcd(:,je,2:kebc)-hzybcd(:,je,2:kebc));  

exybct(:,jb,2:kebc)=caexybct(:,jb,2:kebc).*exybct(:,jb,2:kebc)+...
  cbexybct(:,jb,2:kebc).*(hzxbcb(ibbc:iebc+ie,1,kbbc+kb:kefbc)+...
  hzybcb(ibbc:iebc+ie,1,kbbc+kb:kefbc)-hzxbct(:,je,2:kebc)-hzybct(:,je,2:kebc));

%     LEFT

exybcl(:,2:je,2:kefbc)=caexybcl(:,2:je,2:kefbc).*exybcl(:,2:je,2:kefbc)+...
  cbexybcl(:,2:je,2:kefbc).*(hzxbcl(:,2:je,2:kefbc)+hzybcl(:,2:je,2:kefbc)-...
                             hzxbcl(:,1:je-1,2:kefbc)-hzybcl(:,1:je-1,2:kefbc));

%     RIGHT

exybcr(:,2:je,2:kefbc)=caexybcr(:,2:je,2:kefbc).*exybcr(:,2:je,2:kefbc)+...
  cbexybcr(:,2:je,2:kefbc).*(hzxbcr(:,2:je,2:kefbc)+hzybcr(:,2:je,2:kefbc)-...
                                 hzxbcr(:,1:je-1,2:kefbc)-hzybcr(:,1:je-1,2:kefbc));
 
%   BOTTOM

exybcd(:,2:je,2:kebc)=caexybcd(:,2:je,2:kebc).*exybcd(:,2:je,2:kebc)+...
  cbexybcd(:,2:je,2:kebc).*(hzxbcd(:,2:je,2:kebc)+hzybcd(:,2:je,2:kebc)-...
                                 hzxbcd(:,1:je-1,2:kebc)-hzybcd(:,1:je-1,2:kebc));                           
                             
%   TOP

exybct(:,2:je,2:kebc)=caexybct(:,2:je,2:kebc).*exybct(:,2:je,2:kebc)+...
  cbexybct(:,2:je,2:kebc).*(hzxbct(:,2:je,2:kebc)+hzybct(:,2:je,2:kebc)-...
                                 hzxbct(:,1:je-1,2:kebc)-hzybct(:,1:je-1,2:kebc));
 
%***********************************************************************
%     Update Exz in PML regions
%***********************************************************************


%     FRONT

exzbcf(:,2:jebc,2:kefbc)=caexzbcf(:,2:jebc,2:kefbc).*exzbcf(:,2:jebc,2:kefbc)-...  
  cbexzbcf(:,2:jebc,2:kefbc).*(hyzbcf(:,2:jebc,2:kefbc)+hyxbcf(:,2:jebc,2:kefbc)-...
                      hyzbcf(:,2:jebc,1:kefbc-1)-hyxbcf(:,2:jebc,1:kefbc-1));
                
                      
%     BACK

exzbcb(:,1:jebc,2:kefbc)=caexzbcb(:,1:jebc,2:kefbc).*exzbcb(:,1:jebc,2:kefbc)-...
  cbexzbcb(:,1:jebc,2:kefbc).*(hyzbcb(:,1:jebc,2:kefbc)+hyxbcb(:,1:jebc,2:kefbc)-...
                      hyzbcb(:,1:jebc,1:kefbc-1)-hyxbcb(:,1:jebc,1:kefbc-1));
                    
 
%     LEFT

exzbcl(:,1:jb,2:kefbc)=caexzbcl(:,1:jb,2:kefbc).*exzbcl(:,1:jb,2:kefbc)-...
  cbexzbcl(:,1:jb,2:kefbc).*(hyzbcl(:,1:jb,2:kefbc)+hyxbcl(:,1:jb,2:kefbc)-...
                             hyzbcl(:,1:jb,1:kefbc-1)-hyxbcl(:,1:jb,1:kefbc-1));
                
 
%     RIGHT

exzbcr(:,1:jb,2:kefbc)=caexzbcr(:,1:jb,2:kefbc).*exzbcr(:,1:jb,2:kefbc)-...
  cbexzbcr(:,1:jb,2:kefbc).*(hyzbcr(:,1:jb,2:kefbc)+hyxbcr(:,1:jb,2:kefbc)-...
                             hyzbcr(:,1:jb,1:kefbc-1)-hyxbcr(:,1:jb,1:kefbc-1));
 
 
%   BOTTOM

exzbcd(:,1:jb,2:kebc)=caexzbcd(:,1:jb,2:kebc).*exzbcd(:,1:jb,2:kebc)-...
  cbexzbcd(:,1:jb,2:kebc).*(hyzbcd(:,1:jb,2:kebc)+hyxbcd(:,1:jb,2:kebc)-...
                             hyzbcd(:,1:jb,1:kebc-1)-hyxbcd(:,1:jb,1:kebc-1));
                                 
ex(:,2:je,1)=caex(:,2:je,1).*ex(:,2:je,1)-...
  cbex(:,2:je,1).*(hy(:,2:je,1)-hyzbcd(:,2:je,kebc)-hyxbcd(:,2:je,kebc)); 

%   TOP

exzbct(:,1:je,2:kebc)=caexzbct(:,1:je,2:kebc).*exzbct(:,1:je,2:kebc)-...
  cbexzbct(:,1:je,2:kebc).*(hzxbct(:,1:je,2:kebc)+hzybct(:,1:je,2:kebc)-...
                                 hzxbct(:,1:je,1:kebc-1)-hzybct(:,1:je,1:kebc-1));
                             
ex(:,2:je,kb)=caex(:,2:je,kb).*ex(:,2:je,kb)-...
  cbex(:,2:je,kb).*(hyzbct(:,2:je,1)+hyxbct(:,2:je,1)-hy(:,2:je,1));   
%***********************************************************************
%     Update Eyz in PML regions
%***********************************************************************

%     FRONT

eyzbcf(2:iefbc,:,2:kefbc)=caeyzbcf(2:iefbc,:,2:kefbc).*eyzbcf(2:iefbc,:,2:kefbc)+...
  cbeyzbcf(2:iefbc,:,2:kefbc).*(hxybcf(2:iefbc,:,2:kefbc)+hxzbcf(2:iefbc,:,2:kefbc)-...
                       hxybcf(2:iefbc,:,1:kefbc-1)-hxzbcf(2:iefbc,:,1:kefbc-1));
                     
%     BACK

eyzbcb(2:iefbc,:,2:kefbc)=caeyzbcb(2:iefbc,:,2:kefbc).*eyzbcb(2:iefbc,:,2:kefbc)+...
  cbeyzbcb(2:iefbc,:,2:kefbc).*(hxybcb(2:iefbc,:,2:kefbc)+hxzbcb(2:iefbc,:,2:kefbc)-...
                       hxybcb(2:iefbc,:,1:kefbc-1)-hxzbcb(2:iefbc,:,1:kefbc-1));
 
%     LEFT

eyzbcl(2:iebc,:,2:kefbc)=caeyzbcl(2:iebc,:,2:kefbc).*eyzbcl(2:iebc,:,2:kefbc)+...
  cbeyzbcl(2:iebc,:,2:kefbc).*(hxybcl(2:iebc,:,2:kefbc)+hxzbcl(2:iebc,:,2:kefbc)-...
                       hxybcl(2:iebc,:,1:kefbc-1)-hxzbcl(2:iebc,:,1:kefbc-1));
 
%     RIGHT

eyzbcr(2:iebc,:,2:kefbc)=caeyzbcr(2:iebc,:,2:kefbc).*eyzbcr(2:iebc,:,2:kefbc)+...
  cbeyzbcr(2:iebc,:,2:kefbc).*(hxybcr(2:iebc,:,2:kefbc)+hxzbcr(2:iebc,:,2:kefbc)-...
                       hxybcr(2:iebc,:,1:kefbc-1)-hxzbcr(2:iebc,:,1:kefbc-1));

%   BOTTOM

eyzbcd(:,:,2:kebc)=caeyzbcd(:,:,2:kebc).*eyzbcd(:,:,2:kebc)+...
  cbeyzbcd(:,:,2:kebc).*(hxybcd(:,:,2:kebc)+hxzbcd(:,:,2:kebc)-...
                             hxybcd(:,:,1:kebc-1)-hxzbcd(:,:,1:kebc-1));

ey(2:ie,:,1)=caey(2:ie,:,1).*ey(2:ie,:,1)+...
  cbey(2:ie,:,1).*(hx(2:ie,:,1)-hxybcd(2:ie,:,kebc)-hxzbcd(2:ie,:,kebc));                          
%   TOP

eyzbct(:,:,2:kebc)=caeyzbct(:,:,2:kebc).*eyzbct(:,:,2:kebc)+...
  cbeyzbct(:,:,2:kebc).*(hxybct(:,:,2:kebc)+hxzbct(:,:,2:kebc)-...
                             hxybct(:,:,1:kebc-1)-hxzbct(:,:,1:kebc-1));

 ey(2:ie,:,kb)=caey(2:ie,:,kb).*ey(2:ie,:,kb)+...
  cbey(2:ie,:,kb).*(hxybct(2:ie,:,1)+hxzbct(2:ie,:,1)-hx(2:ie,:,ke));                          
%***********************************************************************
%     Update Eyx in PML regions
%***********************************************************************

%     FRONT

eyxbcf(2:iefbc,:,2:kefbc)=caeyxbcf(2:iefbc,:,2:kefbc).*eyxbcf(2:iefbc,:,2:kefbc)-...
  cbeyxbcf(2:iefbc,:,2:kefbc).*(hzxbcf(2:iefbc,:,2:kefbc)+hzybcf(2:iefbc,:,2:kefbc)-...
                       hzxbcf(1:iefbc-1,:,2:kefbc)-hzybcf(1:iefbc-1,:,2:kefbc));
 
%     BACK

eyxbcb(2:iefbc,:,2:kefbc)=caeyxbcb(2:iefbc,:,2:kefbc).*eyxbcb(2:iefbc,:,2:kefbc)-...
  cbeyxbcb(2:iefbc,:,2:kefbc).*(hzxbcb(2:iefbc,:,2:kefbc)+hzybcb(2:iefbc,:,2:kefbc)-...
                       hzxbcb(1:iefbc-1,:,2:kefbc)-hzybcb(1:iefbc-1,:,2:kefbc));
 
%     LEFT

eyxbcl(2:iebc,:,2:kefbc)=caeyxbcl(2:iebc,:,2:kefbc).*eyxbcl(2:iebc,:,2:kefbc)-...
  cbeyxbcl(2:iebc,:,2:kefbc).*(hzxbcl(2:iebc,:,2:kefbc)+hzybcl(2:iebc,:,2:kefbc)-...
                       hzxbcl(1:iebc-1,:,2:kefbc)-hzybcl(1:iebc-1,:,2:kefbc));
ey(1,:,1:kb)=caey(1,:,1:kb).*ey(1,:,1:kb)+...
  cbey(1,:,1:kb).*(hz(1,:,1:kb)-hzxbcl(iebc,:,kbbc:kebc+kb)-hzybcl(iebc,:,kbbc:kebc+kb));

eyxbcd(1,:,2:kebc)=caeyxbcd(1,:,2:kebc).*eyxbcd(1,:,2:kebc)-...
  cbeyxbcd(1,:,2:kebc).*(hzxbcd(1,:,2:kebc)+hzybcd(1,:,2:kebc)-...
                             hzxbcl(iebc,:,2:kebc)-hzybcl(iebc,:,2:kebc));      

eyxbct(1,:,2:kebc)=caeyxbct(1,:,2:kebc).*eyxbct(1,:,2:kebc)-...
  cbeyxbct(1,:,2:kebc).*(hzxbct(1,:,2:kebc)+hzybct(1,:,2:kebc)-...
  hzxbcl(iebc,:,kbbc+kb:kefbc)-hzybcl(iebc,:,kbbc+kb:kefbc));         
 
%     RIGHT

eyxbcr(2:iebc,:,2:kefbc)=caeyxbcr(2:iebc,:,2:kefbc).*eyxbcr(2:iebc,:,2:kefbc)-...
  cbeyxbcr(2:iebc,:,2:kefbc).*(hzxbcr(2:iebc,:,2:kefbc)+hzybcr(2:iebc,:,2:kefbc)-...
                       hzxbcr(1:iebc-1,:,2:kefbc)-hzybcr(1:iebc-1,:,2:kefbc));
ey(ib,:,1:kb)=caey(ib,:,1:kb).*ey(ib,:,1:kb)+...
  cbey(ib,:,1:kb).*(hzxbcr(1,:,kbbc:kebc+kb)+hzybcr(1,:,kbbc:kebc+kb)-hz(ie,:,1:kb));

eyxbcd(ib,:,2:kebc)=caeyxbcd(ib,:,2:kebc).*eyxbcd(ib,:,2:kebc)-...
  cbeyxbcd(ib,:,2:kebc).*(hzxbcr(1,:,2:kebc)+hzybcr(1,:,2:kebc)-...
                                        hzxbcd(ie,:,2:kebc)-hzybcd(ie,:,2:kebc));   

eyxbct(ib,:,2:kebc)=caeyxbct(ib,:,2:kebc).*eyxbct(ib,:,2:kebc)-...
  cbeyxbct(ib,:,2:kebc).*(hzxbcr(1,:,kbbc+kb:kefbc)+...
  hzybcr(1,:,kbbc+kb:kefbc)-hzxbct(ie,:,2:kebc)-hzybct(ie,:,2:kebc));  

%   BOTTOM

eyxbcd(2:ie,:,2:kebc)=caeyxbcd(2:ie,:,2:kebc).*eyxbcd(2:ie,:,2:kebc)-...
  cbeyxbcd(2:ie,:,2:kebc).*(hzxbcd(2:ie,:,2:kebc)+hzybcd(2:ie,:,2:kebc)-...
                             hzxbcd(1:ie-1,:,2:kebc)-hzybcd(1:ie-1,:,2:kebc));

%   TOP

eyxbct(2:ie,:,2:kebc)=caeyxbct(2:ie,:,2:kebc).*eyxbct(2:ie,:,2:kebc)-...
  cbeyxbct(2:ie,:,2:kebc).*(hzxbct(2:ie,:,2:kebc)+hzybct(2:ie,:,2:kebc)-...
                             hzxbct(1:ie-1,:,2:kebc)-hzybct(1:ie-1,:,2:kebc));
 
%***********************************************************************
%     Update Ezx in PML regions
%***********************************************************************

%     FRONT

ezxbcf(2:iefbc,2:jebc,:)=caezxbcf(2:iefbc,2:jebc,:).*ezxbcf(2:iefbc,2:jebc,:)+...
  cbezxbcf(2:iefbc,2:jebc,:).*(hyzbcf(2:iefbc,2:jebc,:)+hyxbcf(2:iefbc,2:jebc,:)-...
                       hyzbcf(1:iefbc-1,2:jebc,:)-hyxbcf(1:iefbc-1,2:jebc,:));
 
%     BACK

ezxbcb(2:iefbc,2:jebc,:)=caezxbcb(2:iefbc,2:jebc,:).*ezxbcb(2:iefbc,2:jebc,:)+...
  cbezxbcb(2:iefbc,2:jebc,:).*(hyzbcb(2:iefbc,2:jebc,:)+hyxbcb(2:iefbc,2:jebc,:)-...
                       hyzbcb(1:iefbc-1,2:jebc,:)-hyxbcb(1:iefbc-1,2:jebc,:));
 
%     LEFT

ezxbcl(2:iebc,:,:)=caezxbcl(2:iebc,:,:).*ezxbcl(2:iebc,:,:)+...
  cbezxbcl(2:iebc,:,:).*(hyzbcl(2:iebc,:,:)+hyxbcl(2:iebc,:,:)-...
                       hyzbcl(1:iebc-1,:,:)-hyxbcl(1:iebc-1,:,:));
ez(1,:,:)=caez(1,:,:).*ez(1,:,:)+cbez(1,:,:).*(hy(1,:,:)-...
                       hyzbcl(iebc,:,kbbc:kebc+ke)-hyxbcl(iebc,:,kbbc:kebc+ke));                  

ezxbcd(1,:,:)=caezxbcd(1,:,:).*ezxbcd(1,:,:)+...
  cbezxbcd(1,:,:).*(hyzbcd(1,:,:)+hyxbcd(1,:,:)-...
                             hyzbcl(iebc,:,1:kebc)-hyxbcl(iebc,:,1:kebc));  

ezxbct(1,:,:)=caezxbct(1,:,:).*ezxbct(1,:,:)+...
  cbezxbct(1,:,:).*(hyzbct(1,:,:)+hyxbct(1,:,:)-...
  hyzbcl(iebc,:,kebc+kb:kefbc)-hyxbcl(iebc,:,kebc+kb:kefbc));  

%     RIGHT

ezxbcr(2:iebc,:,:)=caezxbcr(2:iebc,:,:).*ezxbcr(2:iebc,:,:)+...
  cbezxbcr(2:iebc,:,:).*(hyzbcr(2:iebc,:,:)+hyxbcr(2:iebc,:,:)-...
                       hyzbcr(1:iebc-1,:,:)-hyxbcr(1:iebc-1,:,:));
ez(ib,:,:)=caez(ib,:,:).*ez(ib,:,:)+cbez(ib,:,:).*(hyzbcr(1,:,kbbc:kebc+ke)+...
                                             hyxbcr(1,:,kbbc:kebc+ke)-hy(ie,:,:));         

ezxbcd(ib,:,:)=caezxbcd(ib,:,:).*ezxbcd(ib,:,:)+...
  cbezxbcd(ib,:,:).*(hyzbcr(1,:,1:kebc)+hyxbcr(1,:,1:kebc)-...
                                               hyzbcd(ie,:,:)+hyxbcd(ie,:,:));                                         

ezxbct(ib,:,:)=caezxbct(ib,:,:).*ezxbct(ib,:,:)+...
  cbezxbct(ib,:,:).*(hyzbcr(1,:,kbbc+ke:kefbc)+...
  hyxbcr(1,:,kbbc+ke:kefbc)-hyzbct(ie,:,:)-hyxbct(ie,:,:));
                                           
%   BOTTOM

ezxbcd(2:ie,:,:)=caezxbcd(2:ie,:,:).*ezxbcd(2:ie,:,:)+...
  cbezxbcd(2:ie,:,:).*(hyzbcd(2:ie,:,:)+hyxbcd(2:ie,:,:)-...
                       hyzbcd(1:ie-1,:,:)-hyxbcd(1:ie-1,:,:));

%   TOP
ezxbct(2:ie,:,:)=caezxbct(2:ie,:,:).*ezxbct(2:ie,:,:)+...
  cbezxbct(2:ie,:,:).*(hyzbct(2:ie,:,:)+hyxbct(2:ie,:,:)-...
                       hyzbct(1:ie-1,:,:)-hyxbct(1:ie-1,:,:));

%***********************************************************************
%     Update Ezy in PML regions
%***********************************************************************

%     FRONT

ezybcf(2:iefbc,2:jebc,:)=caezybcf(2:iefbc,2:jebc,:).*ezybcf(2:iefbc,2:jebc,:)-...
  cbezybcf(2:iefbc,2:jebc,:).*(hxybcf(2:iefbc,2:jebc,:)+hxzbcf(2:iefbc,2:jebc,:)-...
                       hxybcf(2:iefbc,1:jebc-1,:)-hxzbcf(2:iefbc,1:jebc-1,:));

ezybcl(2:iebc,1,:)=caezybcl(2:iebc,1,:).*ezybcl(2:iebc,1,:)-...
  cbezybcl(2:iebc,1,:).*(hxybcl(2:iebc,1,:)+hxzbcl(2:iebc,1,:)-...
                       hxybcf(2:iebc,jebc,:)-hxzbcf(2:iebc,jebc,:));                   
                       
ez(:,1,:)=caez(:,1,:).*ez(:,1,:)+cbez(:,1,:).*(hx(:,1,:)-...
   hxybcf(ibbc:iebc+ib,jebc,kbbc:kebc+ke)-hxzbcf(ibbc:iebc+ib,jebc,kbbc:kebc+ke));

ezybcr(2:iebc,1,:)=caezybcr(2:iebc,1,:).*ezybcr(2:iebc,1,:)-...
  cbezybcr(2:iebc,1,:).*(hxybcr(2:iebc,1,:)+hxzbcr(2:iebc,1,:)-...
                       hxybcf(ibbc+ib:iefbc,jebc,:)-hxzbcf(ibbc+ib:iefbc,jebc,:));

ezybcd(:,1,:)=caezybcd(:,1,:).*ezybcd(:,1,:)-...
  cbezybcd(:,1,:).*(hxybcd(:,1,:)+hxzbcd(:,1,:)-...