bzexy(i,j)=1/(1/dt+Oy/2/eps0);
        hzhz(i,j)=(1-dt*Ox/2/eps0)/(1+dt*Ox/2/eps0);
        hzbz(i,j)=1/mu0/mur/(1+dt*Ox/2/eps0);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% the left area (left PML)  
% coefficient for computing Ex in left
for i=1:nlayerx
    for j=nlayery+1:nlayery+ngmy % the boundary should be treated in the area where the field component is placed really  eg. j==nlayery+1
        epsreff=(epsr(i,j)+epsr(i+1,j))/2; % effective epsr
        if j==nlayery+1
            Oy=epsreff*Om*factory*(0.5*dy)^(m+1);  % special treating for Oy in left-front boundary
        else
            Oy=0;
        end
        x1=(nlayerx-i)*dx;
        x2=(nlayerx+1-i)*dx;
        Ox=epsreff*Om*factorx*(x2^(m+1)-x1^(m+1));
        exex(i,j)=(eps0/dt-Oy/2)/(eps0/dt+Oy/2);
        exdx1(i,j)=(eps0/dt+Ox/2)/(eps0*epsreff*(eps0/dt+Oy/2));
        exdx2(i,j)=(eps0/dt-Ox/2)/(eps0*epsreff*(eps0/dt+Oy/2));
    end
end       

%% coefficient for computing Ey in left 
for i=2:nlayerx %it's PEC boundary for i==1,all coefficients are zero,need not treating
    for j=nlayery+1:nlayery+ngmy
        epsreff=(epsr(i,j)+epsr(i,j+1))/2; % effective epsr
        Oy=0;
        x1=(nlayerx-i+0.5)*dx;
        x2=(nlayerx-i+1.5)*dx;
        Ox=epsreff*Om*factorx*(x2^(m+1)-x1^(m+1));
        eyey(i,j)=(eps0/dt-Ox/2)/(eps0/dt+Ox/2);
        eydy1(i,j)=(eps0/dt+Oy/2)/(eps0*epsreff)/(eps0/dt+Ox/2);
        eydy2(i,j)=(eps0/dt-Oy/2)/(eps0*epsreff)/(eps0/dt+Ox/2);
    end
end        

%% coefficient for computing Hz in left
for i=1:nlayerx
    for j=nlayery+1:nlayery+ngmy
        epsreff=0.25*(epsr(i,j)+epsr(i,j+1)+epsr(i+1,j)+epsr(i+1,j+1));
        Oy=0;
        x1=(nlayerx-i)*dx;
        x2=(nlayerx+1-i)*dx;
        Ox=epsreff*Om*factorx*(x2^(m+1)-x1^(m+1));
        bzbz(i,j)=(1/dt-Oy/2/eps0)/(1/dt+Oy/2/eps0);
        bzexy(i,j)=1/(1/dt+Oy/2/eps0);
        hzhz(i,j)=(1-dt*Ox/2/eps0)/(1+dt*Ox/2/eps0);
        hzbz(i,j)=1/mu0/mur/(1+dt*Ox/2/eps0);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% the right area (right PML)  
% coefficient for computing Ex in right
for i=nlayerx+ngmx+1:ngx
    for j=nlayery+1:nlayery+ngmy % attention: the right-front & right-back boundary already computed 
        epsreff=(epsr(i,j)+epsr(i+1,j))/2; % effective epsr
        if j==nlayery+1
             Oy=epsreff*Om*factory*(0.5*dy)^(m+1);  % special treating for Oy in right-front boundary
        else
             Oy=0;
        end
        x1=(i-nlayerx-ngmx-1)*dx;
        x2=(i-nlayerx-ngmx)*dx;
        Ox=epsreff*Om*factorx*(x2^(m+1)-x1^(m+1));
        exex(i,j)=(eps0/dt-Oy/2)/(eps0/dt+Oy/2);
        exdx1(i,j)=(eps0/dt+Ox/2)/(eps0*epsreff*(eps0/dt+Oy/2));
        exdx2(i,j)=(eps0/dt-Ox/2)/(eps0*epsreff*(eps0/dt+Oy/2));
    end
end       

%% coefficient for computing Ey in right 
for i=nlayerx+ngmx+1:ngx %it's PEC boundary for i==nnx,all coefficients are zero,need not treating
    for j=nlayery+1:nlayery+ngmy
        epsreff=(epsr(i,j)+epsr(i,j+1))/2; % effective epsr
        Oy=0;
        if i==nlayerx+ngmx+1  
            Ox=epsreff*Om*factorx*((0.5*dx)^(m+1)); % special treating right-main boundary
        else
            x1=(i-nlayerx-ngmx-2+0.5)*dx;
            x2=(i-nlayerx-ngmx-2+1.5)*dx;
            Ox=epsreff*Om*factorx*(x2^(m+1)-x1^(m+1));
        end
        eyey(i,j)=(eps0/dt-Ox/2)/(eps0/dt+Ox/2);
        eydy1(i,j)=(eps0/dt+Oy/2)/(eps0*epsreff)/(eps0/dt+Ox/2);
        eydy2(i,j)=(eps0/dt-Oy/2)/(eps0*epsreff)/(eps0/dt+Ox/2);
    end
end        

%% coefficient for computing Hz in right
for i=nlayerx+ngmx+1:ngx
    for j=nlayery+1:nlayery+ngmy
        epsreff=0.25*(epsr(i,j)+epsr(i,j+1)+epsr(i+1,j)+epsr(i+1,j+1));
        Oy=0;
        x1=(i-nlayerx-ngmx-1)*dx;
        x2=(i-nlayerx-ngmx)*dx;
        Ox=epsreff*Om*factorx*(x2^(m+1)-x1^(m+1));
        bzbz(i,j)=(1/dt-Oy/2/eps0)/(1/dt+Oy/2/eps0);
        bzexy(i,j)=1/(1/dt+Oy/2/eps0);
        hzhz(i,j)=(1-dt*Ox/2/eps0)/(1+dt*Ox/2/eps0);
        hzbz(i,j)=1/mu0/mur/(1+dt*Ox/2/eps0);
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% the main area ( normal FDTD area)  

% coefficient for computing Ex in main
for i=nlayerx+1:nlayerx+ngmx
    for j=nlayery+1:nlayery+ngmy 
        epsreff=(epsr(i,j)+epsr(i+1,j))/2; % effective epsr
        Ox=0;
        if j==nlayery+1
           Oy=epsreff*Om*factory*(0.5*dy)^(m+1);  % special treating for Oy in main-front boundary    
        else
           Oy=0;
        end
        exex(i,j)=(eps0/dt-Oy/2)/(eps0/dt+Oy/2);
        exdx1(i,j)=(eps0/dt+Ox/2)/(eps0*epsreff*(eps0/dt+Oy/2));
        exdx2(i,j)=(eps0/dt-Ox/2)/(eps0*epsreff*(eps0/dt+Oy/2));
    end
end 

%% coefficient for computing Ey in main 
for i=nlayerx+1:nlayerx+ngmx
    for j=nlayery+1:nlayery+ngmy 
        epsreff=(epsr(i,j)+epsr(i,j+1))/2; % effective epsr
        if i==nlayerx+1
           Ox=epsreff*Om*factorx*(0.5*dx)^(m+1);  % special treating for Oy in main-left boundary    
        else
           Ox=0;
        end
        Oy=0;
        eyey(i,j)=(eps0/dt-Ox/2)/(eps0/dt+Ox/2);
        eydy1(i,j)=(eps0/dt+Oy/2)/(eps0*epsreff)/(eps0/dt+Ox/2);
        eydy2(i,j)=(eps0/dt-Oy/2)/(eps0*epsreff)/(eps0/dt+Ox/2);
    end
end        

%% coefficient for computing Hz in main
for i=nlayerx+1:nlayerx+ngmx
    for j=nlayery+1:nlayery+ngmy 
        epsreff=0.25*(epsr(i,j)+epsr(i,j+1)+epsr(i+1,j)+epsr(i+1,j+1));
        Ox=0;
        Oy=0;     
        bzbz(i,j)=(1/dt-Oy/2/eps0)/(1/dt+Oy/2/eps0);
        bzexy(i,j)=1/(1/dt+Oy/2/eps0);
        hzhz(i,j)=(1-dt*Ox/2/eps0)/(1+dt*Ox/2/eps0);
        hzbz(i,j)=1/mu0/mur/(1+dt*Ox/2/eps0);
    end
end


%********************* compute E & H field component
for n=1:nt
    Dxstore(1:ngx,1:ngy)=Dx(1:ngx,1:ngy); % save  value of previous time to be used later
    Dx(1:ngx,2:ngy)=dxdx*Dx(1:ngx,2:ngy)+dxhz*(Hz(1:ngx,2:ngy)-Hz(1:ngx,1:ngy-1))/dy;
    Ex(1:ngx,2:ngy)=exex(1:ngx,2:ngy).*Ex(1:ngx,2:ngy)+exdx1(1:ngx,2:ngy).*Dx(1:ngx,2:ngy)-...
        exdx2(1:ngx,2:ngy).*Dxstore(1:ngx,2:ngy);
    Dystore(1:ngx,1:ngy)=Dy(1:ngx,1:ngy);
    Dy(2:ngx,1:ngy)=dydy*Dy(2:ngx,1:ngy)-dyhz*(Hz(2:ngx,1:ngy)-Hz(1:ngx-1,1:ngy))/dx;
    Ey(2:ngx,1:ngy)=eyey(2:ngx,1:ngy).*Ey(2:ngx,1:ngy)+eydy1(2:ngx,1:ngy).*Dy(2:ngx,1:ngy)-...
        eydy2(2:ngx,1:ngy).*Dystore(2:ngx,1:ngy);
    
   % add source
    xssta=nlayerx+8;
    xsend=nlayerx+8;
    yssta=nlayery+10;
    ysend=nlayery+25;    
    Ey(xssta,yssta:ysend)=source(n);
    
    Bzstore(1:ngx,1:ngy)=Bz(1:ngx,1:ngy);
    Bz(1:ngx,1:ngy)=bzbz(1:ngx,1:ngy).*Bz(1:ngx,1:ngy)+bzexy(1:ngx,1:ngy).*((Ex(1:ngx,2:nny)-...
        Ex(1:ngx,1:ngy))/dy-(Ey(2:nnx,1:ngy)-Ey(1:ngx,1:ngy))/dx);
    Hz(1:ngx,1:ngy)=hzhz(1:ngx,1:ngy).*Hz(1:ngx,1:ngy)+hzbz(1:ngx,1:ngy).*(Bz(1:ngx,1:ngy)-...
        Bzstore(1:ngx,1:ngy)); 
  
      %*********************** observe position
    xobs=xssta+10;
    yobs=yssta+10;
    I1(n)=Hz(xobs,yobs);
    I2(n)=Hz(xobs+10,yobs+10);
    I3(n)=Hz(xobs+20,yobs+20);
    Vin(n)=-sum(Ey(xssta,2:ngy)*dy);
    Vtotal1(n)=-sum(Ey(xobs,2:ngy)*dy);
    Vtotal2(n)=-sum(Ey(xobs+15,2:ngy)*dy);
    Vre1(n)=Vtotal1(n)-Vin(n);
    Vre2(n)=Vtotal2(n)-Vin(n);
    R1(n)=-20*log(abs(Vre1(n))/abs(Vin(n)));
    R2(n)=-20*log(abs(Vre2(n))/abs(Vin(n)));
end
clear n;
% t=1*dt*1e9:1*dt*1e9:nt*dt*1e9; % convert time steps into real time (unit: ns)
n=1:nt;
figure(1);
plot(n,I1);
figure(2);
plot(n,I2);
figure(3);
plot(n,I3);
% figure(4);
% plot(n,R1-R2);


% Fmax=40e9;
% n=1:nt;
% k=1;
%   Nf=30;
%   df=Fmax/Nf;
%   while k<=Nf
%       Vinf(k)=sum(Vin.*exp(-i*2*pi*(k*df)*(dt*n)));
%       Vrf(k)=sum(Vr.*exp(-i*2*pi*(k*df)*(dt*n)));
%       k=k+1;
%   end
%   Vinf=dt*Vinf;
%   Vrf=dt*Vrf;
%   R=Vrf./Vinf;