fdtd3d.m

3-D FDTD code with PML absorbing boundary conditions

  hxybcf(ibbc:iebc+ib,jebc,1:kebc)-hxzbcf(ibbc:iebc+ib,jebc,1:kebc));   
                   
                       
ezybct(:,1,:)=caezybct(:,1,:).*ezybct(:,1,:)-...
  cbezybct(:,1,:).*(hxybct(:,1,:)+hxzbct(:,1,:)-...
  hxybcf(ibbc:iebc+ib,jebc,kebc+kb:kefbc)-hxzbcf(ibbc:iebc+ib,jebc,kebc+kb:kefbc));   
                       
%     BACK

ezybcb(2:iefbc,2:jebc,:)=caezybcb(2:iefbc,2:jebc,:).*ezybcb(2:iefbc,2:jebc,:)-...
  cbezybcb(2:iefbc,2:jebc,:).*(hxybcb(2:iefbc,2:jebc,:)+hxzbcb(2:iefbc,2:jebc,:)-...
                       hxybcb(2:iefbc,1:jebc-1,:)-hxzbcb(2:iefbc,1:jebc-1,:));

ezybcl(2:iebc,jb,:)=caezybcl(2:iebc,jb,:).*ezybcl(2:iebc,jb,:)-...
  cbezybcl(2:iebc,jb,:).*(hxybcb(2:iebc,1,:)+hxzbcb(2:iebc,1,:)-...
                                      hxybcl(2:iebc,je,:)-hxzbcl(2:iebc,je,:));                     
                       
ez(:,jb,:)=caez(:,jb,:).*ez(:,jb,:)+cbez(:,jb,:).*(hxybcb(ibbc:iebc+ib,1,kbbc:kebc+ke)+...
                                          hxzbcb(ibbc:iebc+ib,1,kbbc:kebc+ke)-hx(:,je,:));

ezybcr(2:iebc,jb,:)=caezybcr(2:iebc,jb,:).*ezybcr(2:iebc,jb,:)-...
  cbezybcr(2:iebc,jb,:).*(hxybcb(ibbc+ib:iefbc,1,:)+hxzbcb(ibbc+ib:iefbc,1,:)-...
                       hxybcr(2:iebc,je,:)-hxzbcr(2:iebc,je,:)); 
                   
ezybcd(:,jb,:)=caezybcd(:,jb,:).*ezybcd(:,jb,:)-...
  cbezybcd(:,jb,:).*(hxybcb(ibbc:iebc+ib,1,1:kebc)+hxzbcb(ibbc:iebc+ib,1,1:kebc)-...
                                                       hxybcd(:,je,:)-hxzbcd(:,je,:));   
 
ezybct(:,jb,:)=caezybct(:,jb,:).*ezybct(:,jb,:)-...
  cbezybct(:,jb,:).*(hxybcb(ibbc:iebc+ib,1,kebc+kb:kefbc)+...
  hxzbcb(ibbc:iebc+ib,1,kebc+kb:kefbc)-hxybct(:,je,:)-hxzbct(:,je,:));                                                        
                       
%     LEFT

ezybcl(2:iebc,2:je,:)=caezybcl(2:iebc,2:je,:).*ezybcl(2:iebc,2:je,:)-...
  cbezybcl(2:iebc,2:je,:).*(hxybcl(2:iebc,2:je,:)+hxzbcl(2:iebc,2:je,:)-...
                       hxybcl(2:iebc,1:je-1,:)-hxzbcl(2:iebc,1:je-1,:));
                    
%     RIGHT

ezybcr(2:iebc,2:je,:)=caezybcr(2:iebc,2:je,:).*ezybcr(2:iebc,2:je,:)-...
  cbezybcr(2:iebc,2:je,:).*(hxybcr(2:iebc,2:je,:)+hxzbcr(2:iebc,2:je,:)-...
                       hxybcr(2:iebc,1:je-1,:)-hxzbcr(2:iebc,1:je-1,:));
               
%   BOTTOM

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

%   TOP

ezybct(:,2:je,:)=caezybct(:,2:je,:).*ezybct(:,2:je,:)-...
  cbezybct(:,2:je,:).*(hxybct(:,2:je,:)+hxzbct(:,2:je,:)-...
                       hxybct(:,1:je-1,:)-hxzbct(:,1:je-1,:));
 
%***********************************************************************
%     Update magnetic fields (Hx、Hy、Hz) in main grid
%***********************************************************************

hx(:,:,:)=dahx(:,:,:).*hx(:,:,:)-... 
              dbhx(:,:,:).*(ez(:,2:jb,:)-ez(:,1:je,:)+...
                                ey(:,:,1:ke)-ey(:,:,2:kb));

hy(:,:,:)=dahy(:,:,:).*hy(:,:,:)-... 
              dbhy(:,:,:).*(ex(:,:,2:kb)-ex(:,:,1:ke)+...
                                ez(1:ie,:,:)-ez(2:ib,:,:));

hz(:,:,:)=dahz(:,:,:).*hz(:,:,:)-... 
              dbhz(:,:,:).*(ey(2:ib,:,:)-ey(1:ie,:,:)+...
                                 ex(:,1:je,:)-ex(:,2:jb,:));

hz(is,js,ks)=source(n)+hz(is,js,ks);
% hz(is,js,ks)=1;

%***********************************************************************
%     Update Hxy in PML regions
%***********************************************************************

%     FRONT

hxybcf(2:iefbc,1:jebc-1,:)=dahxybcf(2:iefbc,1:jebc-1,:).*hxybcf(2:iefbc,1:jebc-1,:)-...
  dbhxybcf(2:iefbc,1:jebc-1,:).*(ezxbcf(2:iefbc,2:jebc,:)+ezybcf(2:iefbc,2:jebc,:)-...
                              ezxbcf(2:iefbc,1:jebc-1,:)-ezybcf(2:iefbc,1:jebc-1,:));
 
hxybcf(2:iebc,jebc,:)=dahxybcf(2:iebc,jebc,:).*hxybcf(2:iebc,jebc,:)-...
  dbhxybcf(2:iebc,jebc,:).*(ezxbcl(2:iebc,1,:)+ezybcl(2:iebc,1,:)-...
                              ezxbcf(2:iebc,jebc,:)-ezybcf(2:iebc,jebc,:));

hxybcf(ibbc:iebc+ib,jebc,1:kebc)=dahxybcf(ibbc:iebc+ib,jebc,1:kebc).*...
    hxybcf(ibbc:iebc+ib,jebc,1:kebc)-dbhxybcf(ibbc:iebc+ib,jebc,1:kebc).*...
   (ezxbcd(:,1,1:kebc)+ezybcd(:,1,1:kebc)-ezxbcf(ibbc:iebc+ib,jebc,1:kebc)-...
                                              ezybcf(ibbc:iebc+ib,jebc,1:kebc));
                          
hxybcf(ibbc:iebc+ib,jebc,kbbc:kebc+ke)=dahxybcf(ibbc:iebc+ib,jebc,kbbc:kebc+ke).*...
    hxybcf(ibbc:iebc+ib,jebc,kbbc:kebc+ke)-dbhxybcf(ibbc:iebc+ib,jebc,kbbc:kebc+ke).*...
    (ez(:,1,:)-ezxbcf(ibbc:iebc+ib,jebc,kbbc:kebc+ke)-ezybcf(ibbc:iebc+ib,jebc,kbbc:kebc+ke));                              
 
hxybcf(ibbc:iebc+ib,jebc,kebc+kb:kefbc)=dahxybcf(ibbc:iebc+ib,jebc,kebc+kb:kefbc).*...
    hxybcf(ibbc:iebc+ib,jebc,kebc+kb:kefbc)-dbhxybcf(ibbc:iebc+ib,jebc,kebc+kb:kefbc).*...
   (ezxbct(:,1,1:kebc)+ezybct(:,1,1:kebc)-ezxbcf(ibbc:iebc+ib,jebc,kebc+kb:kefbc)-...
                                              ezybcf(ibbc:iebc+ib,jebc,kebc+kb:kefbc));                                                     
                                                     
hxybcf(ibbc+ib:iefbc,jebc,:)=dahxybcf(ibbc+ib:iefbc,jebc,:).*hxybcf(ibbc+ib:iefbc,jebc,:)-...
                 dbhxybcf(ibbc+ib:iefbc,jebc,:).*(ezxbcr(2:iebc,1,:)+ezybcr(2:iebc,1,:)-...
                                ezxbcf(ibbc+ib:iefbc,jebc,:)-ezybcf(ibbc+ib:iefbc,jebc,:));                                                     
                              
%     BACK

hxybcb(2:iefbc,2:jebc,:)=dahxybcb(2:iefbc,2:jebc,:).*hxybcb(2:iefbc,2:jebc,:)-...
  dbhxybcb(2:iefbc,2:jebc,:).*(ezxbcb(2:iefbc,3:jbbc,:)+ezybcb(2:iefbc,3:jbbc,:)-...
                              ezxbcb(2:iefbc,2:jebc,:)-ezybcb(2:iefbc,2:jebc,:));
 
hxybcb(2:iebc,1,:)=dahxybcb(2:iebc,1,:).*hxybcb(2:iebc,1,:)-...
  dbhxybcb(2:iebc,1,:).*(ezxbcb(2:iebc,2,:)+ezybcb(2:iebc,2,:)-...
                                      ezxbcl(2:iebc,jb,:)-ezybcl(2:iebc,jb,:));

hxybcb(ibbc:iebc+ib,1,1:kebc)=dahxybcb(ibbc:iebc+ib,1,1:kebc).*...
    hxybcb(ibbc:iebc+ib,1,1:kebc)-dbhxybcb(ibbc:iebc+ib,1,1:kebc).*...
   (ezxbcb(ibbc:iebc+ib,2,1:kebc)+ezybcb(ibbc:iebc+ib,2,1:kebc)-...
                                ezxbcd(:,jb,1:kebc)-ezybcd(:,jb,1:kebc));
                          
hxybcb(ibbc:iebc+ib,1,kbbc:kebc+ke)=dahxybcb(ibbc:iebc+ib,1,kbbc:kebc+ke).*...
    hxybcb(ibbc:iebc+ib,1,kbbc:kebc+ke)-dbhxybcb(ibbc:iebc+ib,1,kbbc:kebc+ke).*...
    (ezxbcb(ibbc:iebc+ib,2,kbbc:kebc+ke)+ezybcb(ibbc:iebc+ib,2,kbbc:kebc+ke)-ez(:,jb,:)); 
 
hxybcb(ibbc:iebc+ib,1,kebc+kb:kefbc)=dahxybcb(ibbc:iebc+ib,1,kebc+kb:kefbc).*...
    hxybcb(ibbc:iebc+ib,1,kebc+kb:kefbc)-dbhxybcb(ibbc:iebc+ib,1,kebc+kb:kefbc).*...
   (ezxbcb(ibbc:iebc+ib,2,kebc+kb:kefbc)+ezybcb(ibbc:iebc+ib,2,kebc+kb:kefbc)-...
                                              ezxbct(:,jb,1:kebc)-ezybct(:,jb,1:kebc));                                                     
                                                     
hxybcb(ibbc+ib:iefbc,1,:)=dahxybcb(ibbc+ib:iefbc,1,:).*hxybcb(ibbc+ib:iefbc,1,:)-...
                 dbhxybcb(ibbc+ib:iefbc,1,:).*(ezxbcb(ibbc+ib:iefbc,2,:)+...
                 ezybcb(ibbc+ib:iefbc,2,:)-ezxbcr(2:iebc,jb,:)+ezybcr(2:iebc,jb,:)); 

%     LEFT
 
hxybcl(2:iebc,:,:)=dahxybcl(2:iebc,:,:).*hxybcl(2:iebc,:,:)-...
  dbhxybcl(2:iebc,:,:).*(ezxbcl(2:iebc,2:jb,:)+ezybcl(2:iebc,2:jb,:)-...
                              ezxbcl(2:iebc,1:je,:)-ezybcl(2:iebc,1:je,:));

%     RIGHT

hxybcr(2:iebc,:,:)=dahxybcr(2:iebc,:,:).*hxybcr(2:iebc,:,:)-...
  dbhxybcr(2:iebc,:,:).*(ezxbcr(2:iebc,2:jb,:)+ezybcr(2:iebc,2:jb,:)-...
                              ezxbcr(2:iebc,1:je,:)-ezybcr(2:iebc,1:je,:));

%    BOTTOM

hxybcd(:,:,:)=dahxybcd(:,:,:).*hxybcd(:,:,:)-...
  dbhxybcd(:,:,:).*(ezxbcd(:,2:jb,:)+ezybcd(:,2:jb,:)-...
                              ezxbcd(:,1:je,:)-ezybcd(:,1:je,:));

                                                
%    TOP

hxybct(:,:,:)=dahxybct(:,:,:).*hxybct(:,:,:)-...
  dbhxybct(:,:,:).*(ezxbct(:,2:jb,:)+ezybct(:,2:jb,:)-...
                              ezxbct(:,1:je,:)-ezybct(:,1:je,:));

%***********************************************************************
%     Update Hxz in PML regions
%***********************************************************************

%     FRONT

hxzbcf(2:iefbc,:,:)=dahxzbcf(2:iefbc,:,:).*hxzbcf(2:iefbc,:,:)+...
  dbhxzbcf(2:iefbc,:,:).*(eyzbcf(2:iefbc,:,2:kbfbc)+eyxbcf(2:iefbc,:,2:kbfbc)-...
                              eyzbcf(2:iefbc,:,1:kefbc)-eyxbcf(2:iefbc,:,1:kefbc));                       
                              
%     BACK

hxzbcb(2:iefbc,:,:)=dahxzbcb(2:iefbc,:,:).*hxzbcb(2:iefbc,:,:)+...
  dbhxzbcb(2:iefbc,:,:).*(eyzbcb(2:iefbc,:,2:kbfbc)+eyxbcb(2:iefbc,:,2:kbfbc)-...
                              eyzbcb(2:iefbc,:,1:kefbc)-eyxbcb(2:iefbc,:,1:kefbc));

%     LEFT
 
hxzbcl(2:iebc,:,:)=dahxzbcl(2:iebc,:,:).*hxzbcl(2:iebc,:,:)+...
  dbhxzbcl(2:iebc,:,:).*(eyzbcl(2:iebc,:,2:kbfbc)+eyxbcl(2:iebc,:,2:kbfbc)-...
                              eyzbcl(2:iebc,:,1:kefbc)-eyxbcl(2:iebc,:,1:kefbc));

%     RIGHT

hxzbcr(2:iebc,:,:)=dahxzbcr(2:iebc,:,:).*hxzbcr(2:iebc,:,:)+...
  dbhxzbcr(2:iebc,:,:).*(eyzbcr(2:iebc,:,2:kbfbc)+eyxbcr(2:iebc,:,2:kbfbc)-...
                              eyzbcr(2:iebc,:,1:kefbc)-eyxbcr(2:iebc,:,1:kefbc));


%    BOTTOM

hxzbcd(:,:,1:kebc-1)=dahxzbcd(:,:,1:kebc-1).*hxzbcd(:,:,1:kebc-1)+...
  dbhxzbcd(:,:,1:kebc-1).*(eyzbcd(:,:,2:kebc)+eyxbcd(:,:,2:kebc)-...
                              eyzbcd(:,:,1:kebc-1)-eyxbcd(:,:,1:kebc-1));
                          
hxzbcd(:,:,kebc)=dahxzbcd(:,:,kebc).*hxzbcd(:,:,kebc)+...
  dbhxzbcd(:,:,kebc).*(ey(:,:,1)-eyzbcd(:,:,kebc)-eyxbcd(:,:,kebc));                              
                                                
%    TOP

hxzbct(:,:,2:kebc)=dahxzbct(:,:,2:kebc).*hxzbct(:,:,2:kebc)+...
  dbhxzbct(:,:,2:kebc).*(eyzbct(:,:,3:kbbc)+eyxbct(:,:,3:kbbc)-...
                              eyzbct(:,:,2:kebc)-eyxbct(:,:,2:kebc));
                          
hxzbct(:,:,1)=dahxzbct(:,:,1).*hxzbct(:,:,1)+...
  dbhxzbct(:,:,1).*(eyzbct(:,:,2)+eyxbct(:,:,2)-ey(:,:,kb));

%***********************************************************************
%     Update Hyz in PML regions
%***********************************************************************
 
%     FRONT

hyzbcf(:,2:jebc,:)=dahyzbcf(:,2:jebc,:).*hyzbcf(:,2:jebc,:)-...
  dbhyzbcf(:,2:jebc,:).*(exybcf(:,2:jebc,2:kbfbc)+exzbcf(:,2:jebc,2:kbfbc)-...
                              exybcf(:,2:jebc,1:kefbc)-exzbcf(:,2:jebc,1:kefbc));
                          
                              
%     BACK

hyzbcb(:,2:jebc,:)=dahyzbcb(:,2:jebc,:).*hyzbcb(:,2:jebc,:)-...
  dbhyzbcb(:,2:jebc,:).*(exybcb(:,2:jebc,2:kbfbc)+exzbcb(:,2:jebc,2:kbfbc)-...
                              exybcb(:,2:jebc,1:kefbc)-exzbcb(:,2:jebc,1:kefbc));

%     LEFT
 
hyzbcl(:,:,:)=dahyzbcl(:,:,:).*hyzbcl(:,:,:)-...
  dbhyzbcl(:,:,:).*(exybcl(:,:,2:kbfbc)+exzbcl(:,:,2:kbfbc)-...
                              exybcl(:,:,1:kefbc)-exzbcl(:,:,1:kefbc));


%     RIGHT

hyzbcr(:,:,:)=dahyzbcr(:,:,:).*hyzbcr(:,:,:)-...
  dbhyzbcr(:,:,:).*(exybcr(:,:,2:kbfbc)+exzbcr(:,:,2:kbfbc)-...
                              exybcr(:,:,1:kefbc)-exzbcr(:,:,1:kefbc));


%    BOTTOM

hyzbcd(:,:,1:kebc-1)=dahyzbcd(:,:,1:kebc-1).*hyzbcd(:,:,1:kebc-1)-...
  dbhyzbcd(:,:,1:kebc-1).*(exybcd(:,:,2:kebc)+exzbcd(:,:,2:kebc)-...
                              exybcd(:,:,1:kebc-1)-exzbcd(:,:,1:kebc-1));

hyzbcd(:,:,kebc)=dahyzbcd(:,:,kebc).*hyzbcd(:,:,kebc)-...
  dbhyzbcd(:,:,kebc).*(ex(:,:,1)-exybcd(:,:,kebc)-exzbcd(:,:,kebc));                           
                                                
%    TOP

hyzbct(:,:,2:kebc)=dahyzbct(:,:,2:kebc).*hyzbct(:,:,2:kebc)-...
  dbhyzbct(:,:,2:kebc).*(exybct(:,:,3:kbbc)+exzbct(:,:,3:kbbc)-...
                              exybct(:,:,2:kebc)-exzbct(:,:,2:kebc));

hyzbct(:,:,1)=dahyzbct(:,:,1).*hyzbct(:,:,1)-...
  dbhyzbct(:,:,1).*(exybct(:,:,1)+exzbct(:,:,1)-ex(:,:,kb));  


%***********************************************************************
%     Update Hyx in PML regions
%***********************************************************************

%     FRONT

hyxbcf(:,2:jebc,:)=dahyxbcf(:,2:jebc,:).*hyxbcf(:,2:jebc,:)+...
  dbhyxbcf(:,2:jebc,:).*(ezxbcf(2:ibfbc,2:jebc,:)+ezybcf(2:ibfbc,2:jebc,:)-...
                              ezxbcf(1:iefbc,2:jebc,:)-ezybcf(1:iefbc,2:jebc,:));
                          
                              
%     BACK

hyxbcb(:,2:jebc,:)=dahyxbcb(:,2:jebc,:).*hyxbcb(:,2:jebc,:)+...
  dbhyxbcb(:,2:jebc,:).*(ezxbcb(2:ibfbc,2:jebc,:)+ezybcb(2:ibfbc,2:jebc,:)-...
                              ezxbcb(1:iefbc,2:jebc,:)-ezybcb(1:iefbc,2:jebc,:));

%     LEFT
 
hyxbcl(1:iebc-1,:,:)=dahyxbcl(1:iebc-1,:,:).*hyxbcl(1:iebc-1,:,:)+...
  dbhyxbcl(1:iebc-1,:,:).*(ezxbcl(2:iebc,:,:)+ezybcl(2:iebc,:,:)-...
                              ezxbcl(1:iebc-1,:,:)-ezybcl(1:iebc-1,:,:));
                          
hyxbcl(iebc,:,1:kebc)=dahyxbcl(iebc,:,1:kebc).*hyxbcl(iebc,:,1:kebc)+...
  dbhyxbcl(iebc,:,1:kebc).*(ezxbcd(1,:,1:kebc)+ezybcd(1,:,1:kebc)-...
                              ezxbcl(iebc,:,1:kebc)-ezybcl(iebc,:,1:kebc)); 
                          
hyxbcl(iebc,:,kbbc:kebc+ke)=dahyxbcl(iebc,:,kbbc:kebc+ke).*hyxbcl(iebc,:,kbbc:kebc+ke)+...
  dbhyxbcl(iebc,:,kbbc:kebc+ke).*(ez(1,:,1:ke)-ezxbcl(iebc,:,kbbc:kebc+ke)-...
                                                        ezybcl(iebc,:,kbbc:kebc+ke));                               
                              
hyxbcl(iebc,:,kebc+kb:kefbc)=dahyxbcl(iebc,:,kebc+kb:kefbc).*hyxbcl(iebc,:,kebc+kb:kefbc)+...
  dbhyxbcl(iebc,:,kebc+kb:kefbc).*(ezxbct(1,:,1:kebc)+ezybct(1,:,1:kebc)-...
                              ezxbcl(iebc,:,kebc+kb:kefbc)-ezybcl(iebc,:,kebc+kb:kefbc));                             

%     RIGHT

hyxbcr(2:iebc,:,:)=dahyxbcr(2:iebc,:,:).*hyxbcr(2:iebc,:,:)+...
  dbhyxbcr(2:iebc,:,:).*(ezxbcr(3:ibbc,:,:)+ezybcr(3:ibbc,:,:)-...
                              ezxbcr(2:iebc,:,:)-ezybcr(2:iebc,:,:));
                          
hyxbcr(1,:,1:kebc)=dahyxbcr(1,:,1:kebc).*hyxbcr(1,:,1:kebc)+...
  dbhyxbcr(1,:,1:kebc).*(ezxbcr(2,:,1:kebc)+ezybcr(2,:,1:kebc)-...
                                    ezxbcd(ib,:,1:kebc)-ezybcd(ib,:,1:kebc)); 
                          
hyxbcr(1,:,kbbc:kebc+ke)=dahyxbcr(iebc,:,kbbc:kebc+ke).*hyxbcr(iebc,:,kbbc:kebc+ke)+...
  dbhyxbcr(1,:,kbbc:kebc+ke).*(ezxbcr(2,:,kbbc:kebc+ke)+ezybcr(2,:,kbbc:kebc+ke)-...
                                                                          ez(ib,:,1:ke));                               
                              
hyxbcr(1,:,kebc+kb:kefbc)=dahyxbcr(1,:,kebc+kb:kefbc).*hyxbcr(1,:,kebc+kb:kefbc)+...
  dbhyxbcr(1,:,kebc+kb:kefbc).*(ezxbcr(2,:,kebc+kb:kefbc)+ezybcr(2,:,kebc+kb:kefbc)-...
                                                ezxbct(ib,:,1:kebc)-ezybct(ib,:,1:kebc));  

%    BOTTOM

hyxbcd(:,:,:)=dahyxbcd(:,:,:).*hyxbcd(:,:,:)+...
  dbhyxbcd(:,:,:).*(ezxbcd(2:ib,:,:)+ezybcd(2:ib,:,:)-...
                              ezxbcd(1:ie,:,:)-ezybcd(1:ie,:,:));                      
                                                
%    TOP

hyxbct(:,:,:)=dahyxbct(:,:,:).*hyxbct(:,:,:)+...
  dbhyxbct(:,:,:).*(ezxbct(2:ib,:,:)+ezybct(2:ib,:,:)-...
                              ezxbct(1:ie,:,:)-ezybct(1:ie,:,:)); 


%***********************************************************************
%     Update Hzx in PML regions
%***********************************************************************

%     FRONT

hzxbcf(:,:,2:kefbc)=dahzxbcf(:,:,2:kefbc).*hzxbcf(:,:,2:kefbc)-...
  dbhzxbcf(:,:,2:kefbc).*(eyzbcf(2:ibfbc,:,2:kefbc)+eyxbcf(2:ibfbc,:,2:kefbc)-...
                              eyzbcf(1:iefbc,:,2:kefbc)-eyxbcf(1:iefbc,:,2:kefbc));
 
%     BACK

hzxbcb(:,:,2:kefbc)=dahzxbcb(:,:,2:kefbc).*hzxbcb(:,:,2:kefbc)-...
  dbhzxbcb(:,:,2:kefbc).*(eyzbcb(2:ibfbc,:,2:kefbc)+eyxbcb(2:ibfbc,:,2:kefbc)-...
                              eyzbcb(1:iefbc,:,2:kefbc)-eyxbcb(1:iefbc,:,2:kefbc));
 

%     LEFT
 
hzxbcl(1:iebc-1,:,2:kefbc)=dahzxbcl(1:iebc-1,:,2:kefbc).*hzxbcl(1:iebc-1,:,2:kefbc)-...
  dbhzxbcl(1:iebc-1,:,2:kefbc).*(eyzbcl(2:iebc,:,2:kefbc)+eyxbcl(2:iebc,:,2:kefbc)-...
                                eyzbcl(1:iebc-1,:,2:kefbc)-eyxbcl(1:iebc-1,:,2:kefbc));
                            
hzxbcl(iebc,:,2:kebc)=dahzxbcl(iebc,:,2:kebc).*hzxbcl(iebc,:,2:kebc)-...
  dbhzxbcl(iebc,:,2:kebc).*(eyzbcd(1,:,2:kebc)+eyxbcd(1,:,2:kebc)-...
                                eyzbcl(iebc,:,2:kebc)-eyxbcl(iebc,:,2:kebc));
                                
hzxbcl(iebc,:,kbbc:kebc+kb)=dahzxbcl(iebc,:,kbbc:kebc+kb).*hzxbcl(iebc,:,kbbc:kebc+kb)-...
  dbhzxbcl(iebc,:,kbbc:kebc+kb).*(ey(1,:,1:kb)-...
  eyzbcl(iebc,:,kbbc:kebc+kb)-eyxbcl(iebc,:,kbbc:kebc+kb));

hzxbcl(iebc,:,kbbc+kb:kefbc)=dahzxbcl(i