cb1=(1-ca1)/(sigmamax*dx);
caey(1,1:je)=ca1;
cbey(1,1:je)=cb1;
caeybcf(iebc+1,1:jebc)=ca1;
cbeybcf(iebc+1,1:jebc)=cb1;
caeybcb(iebc+1,1:jebc)=ca1;
cbeybcb(iebc+1,1:jebc)=cb1;

for i=1:iebc
  x1=(iebc-i+1)*dx;
  x2=(iebc-i)*dx;
  sigmamax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
  sigmamaxs=sigmamax*(muz/(epsz*eps(1)));
  da1=exp(-sigmamaxs*dt/muz);
  db1=(1-da1)/(sigmamaxs*dx);
  dahzxbcl(i,1:je)=da1;
  dbhzxbcl(i,1:je)=db1;
  dahzxbcf(i,1:jebc)=da1;
  dbhzxbcf(i,1:jebc)=db1;
  dahzxbcb(i,1:jebc)=da1;
  dbhzxbcb(i,1:jebc)=db1;
  caexbcl(i,2:je)=ca(1);
  cbexbcl(i,2:je)=cb(1);
  dahzybcl(i,1:je)=da(1);
  dbhzybcl(i,1:je)=db(1);
end

%     RIGHT region 

caeybcr(ibbc,1:je)=1.0;
cbeybcr(ibbc,1:je)=0.0;
for i=2:iebc
  x1=(i-0.5)*dx;
  x2=(i-1.5)*dx;
  sigmamax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
  ca1=exp(-sigmamax*dt/(epsz*eps(1)));
  cb1=(1-ca1)/(sigmamax*dx);
  caeybcr(i,1:je)=ca1;
  cbeybcr(i,1:je)=cb1;
  caeybcf(i+iebc+ie,1:jebc)=ca1;
  cbeybcf(i+iebc+ie,1:jebc)=cb1;
  caeybcb(i+iebc+ie,1:jebc)=ca1;
  cbeybcb(i+iebc+ie,1:jebc)=cb1;
end
sigmamax=bcfactor*(0.5*dx)^(orderbc+1);
ca1=exp(-sigmamax*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmamax*dx);
caey(ib,1:je)=ca1;
cbey(ib,1:je)=cb1;
caeybcf(iebc+ib,1:jebc)=ca1;
cbeybcf(iebc+ib,1:jebc)=cb1;
caeybcb(iebc+ib,1:jebc)=ca1;
cbeybcb(iebc+ib,1:jebc)=cb1;

for i=1:iebc
  x1=i*dx;
  x2=(i-1)*dx;
  sigmamax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
  sigmamaxs=sigmamax*(muz/(epsz*eps(1)));
  da1=exp(-sigmamaxs*dt/muz);
  db1=(1-da1)/(sigmamaxs*dx);
  dahzxbcr(i,1:je) = da1;
  dbhzxbcr(i,1:je) = db1;
  dahzxbcf(i+ie+iebc,1:jebc)=da1;
  dbhzxbcf(i+ie+iebc,1:jebc)=db1;
  dahzxbcb(i+ie+iebc,1:jebc)=da1;
  dbhzxbcb(i+ie+iebc,1:jebc)=db1;
  caexbcr(i,2:je)=ca(1);
  cbexbcr(i,2:je)=cb(1);
  dahzybcr(i,1:je)=da(1);
  dbhzybcr(i,1:je)=db(1);
end

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

subplot(3,1,1),pcolor(ex');
shading flat;
caxis([-80.0 80.0]);
axis([1 ie 1 jb]);
colorbar;
axis image;
axis off;
title(['Ex at time step = 0']);

subplot(3,1,2),pcolor(ey');
shading flat;
caxis([-80.0 80.0]);
axis([1 ib 1 je]);
colorbar;
axis image;
axis off;
title(['Ey at time step = 0']);

subplot(3,1,3),pcolor(hz');
shading flat;
caxis([-0.2 0.2]);
axis([1 ie 1 je]);
colorbar;
axis image;
axis off;
title(['Hz at time step = 0']);

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

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

%**************************************************************************
%     BEGIN TIME-STEPPING LOOP
%**************************************************************************
for n=1:nmax
    %**********************************************************************
    %     Update electric fields (EX and EY) in main grid
    %**********************************************************************
    ex(:,2:je)=caex(:,2:je).*ex(:,2:je)+...
        cbex(:,2:je).*(hz(:,2:je)-hz(:,1:je-1));
    
    ey(2:ie,:)=caey(2:ie,:).*ey(2:ie,:)+...
        cbey(2:ie,:).*(hz(1:ie-1,:)-hz(2:ie,:));
    
    %**********************************************************************
    %     Update EX in PML regions
    %**********************************************************************
    %     FRONT
    exbcf(:,2:jebc)=caexbcf(:,2:jebc).*exbcf(:,2:jebc)-...
            cbexbcf(:,2:jebc).*(hzxbcf(:,1:jebc-1)+hzybcf(:,1:jebc-1)-...
            hzxbcf(:,2:jebc)-hzybcf(:,2:jebc));
    ex(1:ie,1)=caex(1:ie,1).*ex(1:ie,1)-...
            cbex(1:ie,1).*(hzxbcf(ibbc:iebc+ie,jebc)+...
            hzybcf(ibbc:iebc+ie,jebc)-hz(1:ie,1));
        %     BACK
    exbcb(:,2:jebc-1)=caexbcb(:,2:jebc-1).*exbcb(:,2:jebc-1)-...
            cbexbcb(:,2:jebc-1).*(hzxbcb(:,1:jebc-2)+hzybcb(:,1:jebc-2)-...
            hzxbcb(:,2:jebc-1)-hzybcb(:,2:jebc-1));
    ex(1:ie,jb)=caex(1:ie,jb).*ex(1:ie,jb)-...
            cbex(1:ie,jb).*(hz(1:ie,jb-1)-hzxbcb(ibbc:iebc+ie,1)-...
            hzybcb(ibbc:iebc+ie,1));
    %     LEFT
    exbcl(:,2:je)=caexbcl(:,2:je).*exbcl(:,2:je)-...
            cbexbcl(:,2:je).*(hzxbcl(:,1:je-1)+hzybcl(:,1:je-1)-...
            hzxbcl(:,2:je)-hzybcl(:,2:je));
    exbcl(:,1)=caexbcl(:,1).*exbcl(:,1)-...
            cbexbcl(:,1).*(hzxbcf(1:iebc,jebc)+hzybcf(1:iebc,jebc)-...
            hzxbcl(:,1)-hzybcl(:,1));
    exbcl(:,jb)=caexbcl(:,jb).*exbcl(:,jb)-...
            cbexbcl(:,jb).*(hzxbcl(:,je)+hzybcl(:,je)-...
            hzxbcb(1:iebc,1)-hzybcb(1:iebc,1));
    %     RIGHT
    exbcr(:,2:je)=caexbcr(:,2:je).*exbcr(:,2:je)-...
            cbexbcr(:,2:je).*(hzxbcr(:,1:je-1)+hzybcr(:,1:je-1)-...
            hzxbcr(:,2:je)-hzybcr(:,2:je));
    exbcr(:,1)=caexbcr(:,1).*exbcr(:,1)-...
            cbexbcr(:,1).*(hzxbcf(1+iebc+ie:iefbc,jebc)+...
            hzybcf(1+iebc+ie:iefbc,jebc)-...
            hzxbcr(:,1)-hzybcr(:,1));
    exbcr(:,jb)=caexbcr(:,jb).*exbcr(:,jb)-...
            cbexbcr(:,jb).*(hzxbcr(:,je)+hzybcr(:,je)-...
            hzxbcb(1+iebc+ie:iefbc,1)-...
            hzybcb(1+iebc+ie:iefbc,1));
        
    %**********************************************************************
    %     Update EY in PML regions
    %**********************************************************************
    %     FRONT
    eybcf(2:iefbc,:)=caeybcf(2:iefbc,:).*eybcf(2:iefbc,:)-...
            cbeybcf(2:iefbc,:).*(hzxbcf(2:iefbc,:)+hzybcf(2:iefbc,:)-...
            hzxbcf(1:iefbc-1,:)-hzybcf(1:iefbc-1,:));
        %     BACK
    eybcb(2:iefbc,:)=caeybcb(2:iefbc,:).*eybcb(2:iefbc,:)-...
            cbeybcb(2:iefbc,:).*(hzxbcb(2:iefbc,:)+hzybcb(2:iefbc,:)-...
            hzxbcb(1:iefbc-1,:)-hzybcb(1:iefbc-1,:));
    %     LEFT
    eybcl(2:iebc,:)=caeybcl(2:iebc,:).*eybcl(2:iebc,:)-...
            cbeybcl(2:iebc,:).*(hzxbcl(2:iebc,:)+hzybcl(2:iebc,:)-...
            hzxbcl(1:iebc-1,:)-hzybcl(1:iebc-1,:));
    ey(1,:)=caey(1,:).*ey(1,:)-...
            cbey(1,:).*(hz(1,:)-hzxbcl(iebc,:)-hzybcl(iebc,:));
    %     RIGHT
    eybcr(2:iebc,:)=caeybcr(2:iebc,:).*eybcr(2:iebc,:)-...
            cbeybcr(2:iebc,:).*(hzxbcr(2:iebc,:)+hzybcr(2:iebc,:)-...
            hzxbcr(1:iebc-1,:)-hzybcr(1:iebc-1,:));
    ey(ib,:)=caey(ib,:).*ey(ib,:)-...
            cbey(ib,:).*(hzxbcr(1,:)+hzybcr(1,:)- hz(ie,:));
    %**********************************************************************
    %     Update magnetic fields (HZ) in main grid
    %**********************************************************************
    hz(1:ie,1:je)=dahz(1:ie,1:je).*hz(1:ie,1:je)+...
            dbhz(1:ie,1:je).*(ex(1:ie,2:jb)-ex(1:ie,1:je)+...
            ey(1:ie,1:je)-ey(2:ib,1:je));
    hz(is,js)=source(n);
    %**********************************************************************
    %     Update HZX in PML regions
    %**********************************************************************
    %     FRONT
    hzxbcf(1:iefbc,:)=dahzxbcf(1:iefbc,:).*hzxbcf(1:iefbc,:)-...
            dbhzxbcf(1:iefbc,:).*(eybcf(2:ibfbc,:)-eybcf(1:iefbc,:));
    %     BACK
    hzxbcb(1:iefbc,:)=dahzxbcb(1:iefbc,:).*hzxbcb(1:iefbc,:)-...
            dbhzxbcb(1:iefbc,:).*(eybcb(2:ibfbc,:)-eybcb(1:iefbc,:));
    %     LEFT
    hzxbcl(1:iebc-1,:)=dahzxbcl(1:iebc-1,:).*hzxbcl(1:iebc-1,:)-...
            dbhzxbcl(1:iebc-1,:).*(eybcl(2:iebc,:)-eybcl(1:iebc-1,:));
    hzxbcl(iebc,:)=dahzxbcl(iebc,:).*hzxbcl(iebc,:)-...
            dbhzxbcl(iebc,:).*(ey(1,:)-eybcl(iebc,:));
    %     RIGHT
    hzxbcr(2:iebc,:)=dahzxbcr(2:iebc,:).*hzxbcr(2:iebc,:)-...
            dbhzxbcr(2:iebc,:).*(eybcr(3:ibbc,:)-eybcr(2:iebc,:));
    hzxbcr(1,:)=dahzxbcr(1,:).*hzxbcr(1,:)-...
            dbhzxbcr(1,:).*(eybcr(2,:)-ey(ib,:));
    %**********************************************************************
    %     Update HZY in PML regions
    %**********************************************************************
    %     FRONT
    hzybcf(:,1:jebc-1)=dahzybcf(:,1:jebc-1).*hzybcf(:,1:jebc-1)-...
            dbhzybcf(:,1:jebc-1).*(exbcf(:,1:jebc-1)-exbcf(:,2:jebc));
    hzybcf(1:iebc,jebc)=dahzybcf(1:iebc,jebc).*hzybcf(1:iebc,jebc)-...
            dbhzybcf(1:iebc,jebc).*(exbcf(1:iebc,jebc)-exbcl(1:iebc,1));
    hzybcf(iebc+1:iebc+ie,jebc)=...
            dahzybcf(iebc+1:iebc+ie,jebc).*hzybcf(iebc+1:iebc+ie,jebc)-...
            dbhzybcf(iebc+1:iebc+ie,jebc).*(exbcf(iebc+1:iebc+ie,jebc)-...
            ex(1:ie,1));
    hzybcf(iebc+ie+1:iefbc,jebc)=...
            dahzybcf(iebc+ie+1:iefbc,jebc).*hzybcf(iebc+ie+1:iefbc,jebc)-...
            dbhzybcf(iebc+ie+1:iefbc,jebc).*(exbcf(iebc+ie+1:iefbc,jebc)-...
            exbcr(1:iebc,1));
    %     BACK
    hzybcb(1:iefbc,2:jebc)=dahzybcb(1:iefbc,2:jebc).*hzybcb(1:iefbc,2:jebc)-...
            dbhzybcb(1:iefbc,2:jebc).*(exbcb(1:iefbc,2:jebc)-exbcb(1:iefbc,3:jbbc));
    hzybcb(1:iebc,1)=dahzybcb(1:iebc,1).*hzybcb(1:iebc,1)-...
            dbhzybcb(1:iebc,1).*(exbcl(1:iebc,jb)-exbcb(1:iebc,2));
    hzybcb(iebc+1:iebc+ie,1)=...
            dahzybcb(iebc+1:iebc+ie,1).*hzybcb(iebc+1:iebc+ie,1)-...
            dbhzybcb(iebc+1:iebc+ie,1).*(ex(1:ie,jb)-exbcb(iebc+1:iebc+ie,2));
    hzybcb(iebc+ie+1:iefbc,1)=...
            dahzybcb(iebc+ie+1:iefbc,1).*hzybcb(iebc+ie+1:iefbc,1)-...
            dbhzybcb(iebc+ie+1:iefbc,1).*(exbcr(1:iebc,jb)-...
            exbcb(iebc+ie+1:iefbc,2));
    %     LEFT
    hzybcl(:,1:je)=dahzybcl(:,1:je).*hzybcl(:,1:je)-...
           dbhzybcl(:,1:je).*(exbcl(:,1:je)-exbcl(:,2:jb));
       
    %     RIGHT
       
    hzybcr(:,1:je)=dahzybcr(:,1:je).*hzybcr(:,1:je)-...
           dbhzybcr(:,1:je).*(exbcr(:,1:je)-exbcr(:,2:jb));
       
    %**********************************************************************
    %     Visualize fields
    %**********************************************************************
    if mod(n,4)==0;
        timestep=int2str(n);
        
        subplot(3,1,1),pcolor(ex');
        shading flat;
        caxis([-80.0 80.0]);
        axis([1 ie 1 jb]);
        colorbar;
        axis image;
        axis off;
        title(['Ex at time step = ',timestep]);
        
        subplot(3,1,2),pcolor(ey');
        shading flat;
        caxis([-80.0 80.0]);
        axis([1 ib 1 je]);
        colorbar;
        axis image;
        axis off;
        title(['Ey at time step = ',timestep]);
        
        subplot(3,1,3),pcolor(hz');
        shading flat;
        caxis([-0.2 0.2]);
        axis([1 ie 1 je]);
        colorbar;
        axis image;
        axis off;
        title(['Hz at time step = ',timestep]);
        nn=n/4;
        M(:,nn)=getframe(gcf,rect);
    end;
    %**********************************************************************
    %     END TIME-STEPPING LOOP
    %**********************************************************************
end
movie(gcf,M,0,10,rect);