%     where i varies over the PML region.
%  
%***********************************************************************

rmax=exp(-16);  %desired reflection error, designated as R(0) in Equation 7.62 

orderbc=4;      %order of the polynomial grading, designated as m in Equation 7.60a,b

%   x-varying material properties
delbc=upml*delta;
sigmam=-log(rmax)*(orderbc+1.0)/(2.0*eta*delbc); 
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1.0));
kmax=1;
kfactor=(kmax-1.0)/delta/(orderbc+1.0)/delbc^orderbc;

for i=1:upml
    
    % Coefficients for field components in the center of the grid cell
    x1=(upml-i+1)*delta;
    x2=(upml-i)*delta;
    sigma=sigfactor*(x1^(orderbc+1)-x2^(orderbc+1));
    ki=1+kfactor*(x1^(orderbc+1)-x2^(orderbc+1));
    facm=(2*epsr*epsz*ki-sigma*dt);
    facp=(2*epsr*epsz*ki+sigma*dt);

    C5ex(i,:,:)=facp;
    C5ex(ie_tot-i+1,:,:)=facp;
    C6ex(i,:,:)=facm;
    C6ex(ie_tot-i+1,:,:)=facm;
    D1hz(i,:,:)=facm/facp;
    D1hz(ie_tot-i+1,:,:)=facm/facp;
    D2hz(i,:,:)=2.0*epsr*epsz*dt/facp;
    D2hz(ie_tot-i+1,:,:)=2.0*epsr*epsz*dt/facp;
    D3hy(i,:,:)=facm/facp;
    D3hy(ie_tot-i+1,:,:)=facm/facp;
    D4hy(i,:,:)=1.0/facp/mur/muz;
    D4hy(ie_tot-i+1,:,:)=1.0/facp/mur/muz;

    % Coefficients for field components on the grid cell boundary
    x1=(upml-i+1.5)*delta;
    x2=(upml-i+0.5)*delta;
    sigma=sigfactor*(x1^(orderbc+1)-x2^(orderbc+1));
    ki=1.0+kfactor*(x1^(orderbc+1)-x2^(orderbc+1));
    facm=(2.0*epsr*epsz*ki-sigma*dt);
    facp=(2.0*epsr*epsz*ki+sigma*dt);

    C1ez(i,:,:)=facm/facp;
    C1ez(ih_tot-i+1,:,:)=facm/facp;
    C2ez(i,:,:)=2.0*epsr*epsz*dt/facp;
    C2ez(ih_tot-i+1,:,:)=2.0*epsr*epsz*dt/facp;
    C3ey(i,:,:)=facm/facp;
    C3ey(ih_tot-i+1,:,:)=facm/facp;
    C4ey(i,:,:)=1.0/facp/epsr/epsz;
    C4ey(ih_tot-i+1,:,:)=1.0/facp/epsr/epsz;
    D5hx(i,:,:)=facp;
    D5hx(ih_tot-i+1,:,:)=facp;
    D6hx(i,:,:)=facm;
    D6hx(ih_tot-i+1,:,:)=facm;
    
end

%   PEC walls
C1ez(1,:,:)=-1.0;
C1ez(ih_tot,:,:)=-1.0;
C2ez(1,:,:)=0.0;
C2ez(ih_tot,:,:)=0.0;
C3ey(1,:,:)=-1.0;
C3ey(ih_tot,:,:)=-1.0;
C4ey(1,:,:)=0.0;
C4ey(ih_tot,:,:)=0.0;

%   y-varying material properties
delbc=upml*delta;
sigmam=-log(rmax)*epsr*epsz*cc*(orderbc+1.0)/(2.0*delbc); 
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1.0));
kmax=1.0;
kfactor=(kmax-1.0)/delta/(orderbc+1.0)/delbc^orderbc;

for j=1:upml
    
    % Coefficients for field components in the center of the grid cell
    y1=(upml-j+1)*delta;
    y2=(upml-j)*delta;
    sigma=sigfactor*(y1^(orderbc+1)-y2^(orderbc+1));
    ki=1+kfactor*(y1^(orderbc+1)-y2^(orderbc+1));
    facm=(2*epsr*epsz*ki-sigma*dt);
    facp=(2*epsr*epsz*ki+sigma*dt);
    
    C5ey(:,j,:)=facp;
    C5ey(:,je_tot-j+1,:)=facp;
    C6ey(:,j,:)=facm;
    C6ey(:,je_tot-j+1,:)=facm;
    D1hx(:,j,:)=facm/facp;
    D1hx(:,je_tot-j+1,:)=facm/facp;
    D2hx(:,j,:)=2*epsr*epsz*dt/facp;
    D2hx(:,je_tot-j+1,:)=2*epsr*epsz*dt/facp;
    D3hz(:,j,:)=facm/facp;
    D3hz(:,je_tot-j+1,:)=facm/facp;
    D4hz(:,j,:)=1/facp/mur/muz;
    D4hz(:,je_tot-j+1,:)=1/facp/mur/muz;
    
    % Coefficients for field components on the grid cell boundary
    y1=(upml-j+1.5)*delta;
    y2=(upml-j+0.5)*delta;
    sigma=sigfactor*(y1^(orderbc+1)-y2^(orderbc+1));
    ki=1+kfactor*(y1^(orderbc+1)-y2^(orderbc+1));
    facm=(2*epsr*epsz*ki-sigma*dt);
    facp=(2*epsr*epsz*ki+sigma*dt);    
     
    C1ex(:,j,:)=facm/facp;
    C1ex(:,jh_tot-j+1,:)=facm/facp;
    C2ex(:,j,:)=2*epsr*epsz*dt/facp;
    C2ex(:,jh_tot-j+1,:)=2*epsr*epsz*dt/facp;
    C3ez(:,j,:)=facm/facp;
    C3ez(:,jh_tot-j+1,:)=facm/facp;
    C4ez(:,j,:)=1/facp/epsr/epsz;
    C4ez(:,jh_tot-j+1,:)=1/facp/epsr/epsz;   
    D5hy(:,j,:)=facp;
    D5hy(:,jh_tot-j+1,:)=facp;
    D6hy(:,j,:)=facm;
    D6hy(:,jh_tot-j+1,:)=facm;

end

%   PEC walls
C1ex(:,1,:)=-1;
C1ex(:,jh_tot,:)=-1;
C2ex(:,1,:)=0;
C2ex(:,jh_tot,:)=0;
C3ez(:,1,:)=-1;
C3ez(:,jh_tot,:)=-1;
C4ez(:,1,:)=0;
C4ez(:,jh_tot,:)=0;   

%   z-varying material properties
delbc=upml*delta;
sigmam=-log(rmax)*epsr*epsz*cc*(orderbc+1)/(2*delbc); 
sigfactor=sigmam/(delta*(delbc^orderbc)*(orderbc+1));
kmax=1;
kfactor=(kmax-1)/delta/(orderbc+1)/delbc^orderbc;

for k=1:upml

    % Coefficients for field components in the center of the grid cell
    z1=(upml-k+1)*delta;
    z2=(upml-k)*delta;
    sigma=sigfactor*(z1^(orderbc+1)-z2^(orderbc+1));
    ki=1+kfactor*(z1^(orderbc+1)-z2^(orderbc+1));
    facm=(2*epsr*epsz*ki-sigma*dt);
    facp=(2*epsr*epsz*ki+sigma*dt);
    
    C5ez(:,:,k)=facp;
    C5ez(:,:,ke_tot-k+1)=facp;
    C6ez(:,:,k)=facm;
    C6ez(:,:,ke_tot-k+1)=facm;
    D1hy(:,:,k)=facm/facp;
    D1hy(:,:,ke_tot-k+1)=facm/facp;
    D2hy(:,:,k)=2*epsr*epsz*dt/facp;
    D2hy(:,:,ke_tot-k+1)=2*epsr*epsz*dt/facp;
    D3hx(:,:,k)=facm/facp;
    D3hx(:,:,ke_tot-k+1)=facm/facp;
    D4hx(:,:,k)=1/facp/mur/muz;
    D4hx(:,:,ke_tot-k+1)=1/facp/mur/muz;
    
    % Coefficients for field components on the grid cell boundary
    z1=(upml-k+1.5)*delta;
    z2=(upml-k+0.5)*delta;
    sigma=sigfactor*(z1^(orderbc+1)-z2^(orderbc+1));
    ki=1+kfactor*(z1^(orderbc+1)-z2^(orderbc+1));
    facm=(2*epsr*epsz*ki-sigma*dt);
    facp=(2*epsr*epsz*ki+sigma*dt);
    
    C1ey(:,:,k)=facm/facp;
    C1ey(:,:,kh_tot-k+1)=facm/facp;
    C2ey(:,:,k)=2*epsr*epsz*dt/facp;
    C2ey(:,:,kh_tot-k+1)=2*epsr*epsz*dt/facp;
    C3ex(:,:,k)=facm/facp;
    C3ex(:,:,kh_tot-k+1)=facm/facp;
    C4ex(:,:,k)=1/facp/epsr/epsz;
    C4ex(:,:,kh_tot-k+1)=1/facp/epsr/epsz;
    D5hz(:,:,k)=facp;
    D5hz(:,:,kh_tot-k+1)=facp;
    D6hz(:,:,k)=facm;
    D6hz(:,:,kh_tot-k+1)=facm;

end

%   PEC walls
C1ey(:,:,1)=-1;
C1ey(:,:,kh_tot)=-1;
C2ey(:,:,1)=0;
C2ey(:,:,kh_tot)=0;
C3ex(:,:,1)=-1;
C3ex(:,:,kh_tot)=-1;
C4ex(:,:,1)=0;
C4ex(:,:,kh_tot)=0;

%figure
%set(gcf,'DoubleBuffer','on')

%***********************************************************************
%     Begin time stepping loop
%***********************************************************************

for n=1:nmax
    
    % Update magnetic field
    bstore=bx;
    bx(2:ie_tot,:,:)=D1hx(2:ie_tot,:,:).*  bx(2:ie_tot,:,:)-...
                     D2hx(2:ie_tot,:,:).*((ez(2:ie_tot,2:jh_tot,:)-ez(2:ie_tot,1:je_tot,:))-...
                                          (ey(2:ie_tot,:,2:kh_tot)-ey(2:ie_tot,:,1:ke_tot)))./delta;
    hx(2:ie_tot,:,:)= D3hx(2:ie_tot,:,:).*hx(2:ie_tot,:,:)+...
                      D4hx(2:ie_tot,:,:).*(D5hx(2:ie_tot,:,:).*bx(2:ie_tot,:,:)-...
                                           D6hx(2:ie_tot,:,:).*bstore(2:ie_tot,:,:));
    bstore=by;
    by(:,2:je_tot,:)=D1hy(:,2:je_tot,:).*  by(:,2:je_tot,:)-...
                     D2hy(:,2:je_tot,:).*((ex(:,2:je_tot,2:kh_tot)-ex(:,2:je_tot,1:ke_tot))-...
                                          (ez(2:ih_tot,2:je_tot,:)-ez(1:ie_tot,2:je_tot,:)))./delta;
    hy(:,2:je_tot,:)= D3hy(:,2:je_tot,:).*hy(:,2:je_tot,:)+...
                      D4hy(:,2:je_tot,:).*(D5hy(:,2:je_tot,:).*by(:,2:je_tot,:)-...
                                           D6hy(:,2:je_tot,:).*bstore(:,2:je_tot,:));
    bstore=bz;
    bz(:,:,2:ke_tot)=D1hz(:,:,2:ke_tot).*  bz(:,:,2:ke_tot)-...
                     D2hz(:,:,2:ke_tot).*((ey(2:ih_tot,:,2:ke_tot)-ey(1:ie_tot,:,2:ke_tot))-...
                                          (ex(:,2:jh_tot,2:ke_tot)-ex(:,1:je_tot,2:ke_tot)))./delta;
    hz(:,:,2:ke_tot)= D3hz(:,:,2:ke_tot).*hz(:,:,2:ke_tot)+...
                      D4hz(:,:,2:ke_tot).*(D5hz(:,:,2:ke_tot).*bz(:,:,2:ke_tot)-...
                                           D6hz(:,:,2:ke_tot).*bstore(:,:,2:ke_tot));
    
    % Update electric field
    dstore=dx;
    dx(:,2:je_tot,2:ke_tot)=C1ex(:,2:je_tot,2:ke_tot).*  dx(:,2:je_tot,2:ke_tot)+...
                            C2ex(:,2:je_tot,2:ke_tot).*((hz(:,2:je_tot,2:ke_tot)-hz(:,1:je_tot-1,2:ke_tot))-...
                                                        (hy(:,2:je_tot,2:ke_tot)-hy(:,2:je_tot,1:ke_tot-1)))./delta;
    ex(:,2:je_tot,2:ke_tot)=C3ex(:,2:je_tot,2:ke_tot).*ex(:,2:je_tot,2:ke_tot)+...
                            C4ex(:,2:je_tot,2:ke_tot).*(C5ex(:,2:je_tot,2:ke_tot).*dx(:,2:je_tot,2:ke_tot)-...
                                                        C6ex(:,2:je_tot,2:ke_tot).*dstore(:,2:je_tot,2:ke_tot));
    dstore=dy;
    dy(2:ie_tot,:,2:ke_tot)=C1ey(2:ie_tot,:,2:ke_tot).*  dy(2:ie_tot,:,2:ke_tot)+...
                            C2ey(2:ie_tot,:,2:ke_tot).*((hx(2:ie_tot,:,2:ke_tot)-hx(2:ie_tot,:,1:ke_tot-1))-...
                                                        (hz(2:ie_tot,:,2:ke_tot)-hz(1:ie_tot-1,:,2:ke_tot)))./delta;
    ey(2:ie_tot,:,2:ke_tot)=C3ey(2:ie_tot,:,2:ke_tot).*ey(2:ie_tot,:,2:ke_tot)+...
                            C4ey(2:ie_tot,:,2:ke_tot).*(C5ey(2:ie_tot,:,2:ke_tot).*dy(2:ie_tot,:,2:ke_tot)-...
                                                        C6ey(2:ie_tot,:,2:ke_tot).*dstore(2:ie_tot,:,2:ke_tot));
    dstore=dz;
    dz(2:ie_tot,2:je_tot,:)=C1ez(2:ie_tot,2:je_tot,:).*  dz(2:ie_tot,2:je_tot,:)+...
                            C2ez(2:ie_tot,2:je_tot,:).*((hy(2:ie_tot,2:je_tot,:)-hy(1:ie_tot-1,2:je_tot,:))-...
                                                        (hx(2:ie_tot,2:je_tot,:)-hx(2:ie_tot,1:je_tot-1,:)))./delta;
    dz(is,js,ks:ks+1)=dz(is,js,ks:ks+1)+J0*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));
    ez(2:ie_tot,2:je_tot,:)=C3ez(2:ie_tot,2:je_tot,:).*ez(2:ie_tot,2:je_tot,:)+...
                            C4ez(2:ie_tot,2:je_tot,:).*(C5ez(2:ie_tot,2:je_tot,:).*dz(2:ie_tot,2:je_tot,:)-...
                                                        C6ez(2:ie_tot,2:je_tot,:).*dstore(2:ie_tot,2:je_tot,:));

    %***********************************************************************
    %     Visualize fields
    %***********************************************************************

    timestep=int2str(n);
    tview(:,:)=squeeze(ez(ih_bc:upml+ie,jh_bc:upml+je,ks));
    sview(:,:)=squeeze(ez(ih_bc:upml+ie,js,kh_bc:upml+ke));
    
    subplot('position',[0.15 0.57 0.7 0.35])
    %imagesc(tview'); 
    pcolor(tview');
    shading interp;
    caxis([-0.2 0.2]);
    colorbar;
    axis image; axis xy;
    title(['E_z(i,j,k=k_s_o_u_r_c_e), time step = ',timestep]);
    xlabel('i coordinate');
    ylabel('j coordinate');
    
    subplot('position',[0.15 0.08 0.7 0.35])
    %imagesc(tview');
    pcolor(tview');
    shading interp;
    caxis([-0.2 0.2]);
    colorbar;
    axis image; axis xy;
    title(['E_z(i,j=j_s_o_u_r_c_e,k), time step = ',timestep]);
    xlabel('i coordinate');
    ylabel('k coordinate');
    
    pause(0.05)
    
end