% 2-D TE mode FDTD with UPML-ABC
clear all;
nlayerx=8;  % number of PML in x-axis
nlayery=8;
% some physical constant
mu0=4.0*pi*1.0e-7; % permeability in free space
c=2.99792458e8;    % speed of light in free space
eps0=1.0/(c^2*mu0);%permittivity in free space

%basic parameters for the simulated area
ngmx=80;  % number of grid cells in x-axis in the main FDTD area
ngmy=80;
dx=0.02e-6; dy=0.02e-6; % step of grid cells in x; y-axis
nnmx=ngmx+1; %number of nodes in x-axis in the main FDTD area
nnmy=ngmy+1;
ngx=ngmx+2*nlayerx; % number of grid cells in x-axis for the whole simulated area
ngy=ngmy+2*nlayery;
nnx=ngx+1;  % number of nodes in x-axis for the whole simulated area
nny=ngy+1;
nt=3000; %number of time steps
dt=dx/(2*c); % step of time
epsr=ones(nnx,nny); % relative permittivity in simulated area  
mur=1;               % relative permeability in simulated area

% allocate memory for coefficient matrix
dxdx=1; % coefficient of dx for computing dx
dxhz=dt;% coefficient of hz for computing dx
exex=zeros(nnx,nny);
exdx1=zeros(nnx,nny);
exdx2=zeros(nnx,nny);
dydy=1;
dyhz=dt;
eyey=zeros(nnx,nny);
eydy1=zeros(nnx,nny);
eydy2=zeros(nnx,nny);
bzbz=zeros(nnx,nny);
bzexy=zeros(nnx,nny);
hzhz=zeros(nnx,nny);
hzbz=zeros(nnx,nny);

% allocate memory for field component
Dx=zeros(nnx,nny);
Ex=zeros(nnx,nny);
Dy=zeros(nnx,nny);
Ey=zeros(nnx,nny);
Bz=zeros(nnx,nny);
Hz=zeros(nnx,nny);
% instruction: we can write size of Ex Ey and Hz in a uniform formula (nnx,nny).The component outside ...
%  our simulated area are zeros as we have gaven,that is Ex(nnx,:)==0,Ey(:,nny)==0,Hz(nnx,:)==0,Hz(:,nny)...
%  ==0 and we need not operate on them 

% wave excitation
f0=193*10^12;
lambda=c/f0;
omega=2.0*pi*f0;
T=1/f0;
t0=3*T;
n0=t0/dt;
source=zeros(1,nt);
for n=1:7.0*(T/dt)
    source(n)=sin(omega*(n-n0)*dt)*exp(-((n-n0)^2/(T/dt)^2)); % modulated sine wave
end


rmax=10^(-8); % the maximal reflect coefficient
m=4; % the order of UPML
Om=-log(rmax/100)*eps0*c*(m+1)/(2*dx*nlayerx);  % unified to epsr(已對epsr歸一化）
factorx=1/(dx*(nlayerx*dx)^m*(m+1));
factory=1/(dy*(nlayery*dy)^m*(m+1));

%************************************************************************** decide the coefficient matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% the front area (under PML)
% coefficient for computing Ex in front
for i=1:ngx
    for j=2:nlayery
        epsreff=(epsr(i,j)+epsr(i+1,j))/2; % effective epsr
        if i<=nlayerx
            x1=(nlayerx-i)*dx;
            x2=(nlayerx+1-i)*dx;
            Ox=epsreff*Om*factorx*(x2^(m+1)-x1^(m+1));
        else if i>=nlayerx+ngmx+1
                 x1=(i-nlayerx-ngmx-1)*dx;
                 x2=(i-nlayerx-ngmx)*dx;
                 Ox=epsreff*Om*factorx*(x2^(m+1)-x1^(m+1));
            else
                Ox=0;
            end
        end
        y1=(nlayery-j+0.5)*dy;
        y2=(nlayery-j+1.5)*dy;
        Oy=epsreff*Om*factory*(y2^(m+1)-y1^(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 front
for i=2:nnx-1 %it's PEC boundary for i==1 & i==nny ,all coefficients are zero,need not treating
    for j=1:nlayery
        epsreff=(epsr(i,j)+epsr(i,j+1))/2; % effective epsr
        if i<=nlayerx
            x1=(nlayerx+0.5-i)*dx;
            x2=(nlayerx+1.5-i)*dx;
            Ox=epsreff*Om*factorx*(x2^(m+1)-x1^(m+1));
        else if i>=nlayerx+ngmx+2
                 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));
            else if i==nlayerx+1|i==nlayerx+ngmx+1  % special treating for Ox in two boundaries in front PML 
                    Ox=epsreff*Om*factorx*((0.5*dx)^(m+1));% substitute r=0 with r=0.5*dx to make sigma decreased  more slowly  
                else
                   Ox=0;
                end
            end
        end
        y1=(nlayery-j)*dy;
        y2=(nlayery+1-j)*dy;
        Oy=epsreff*Om*factory*(y2^(m+1)-y1^(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 front
for i=1:ngx
    for j=1:nlayery
        epsreff=0.25*(epsr(i,j)+epsr(i,j+1)+epsr(i+1,j)+epsr(i+1,j+1));
        if i<=nlayerx
            x1=(nlayerx-i)*dx;
            x2=(nlayerx+1-i)*dx;
            Ox=epsreff*Om*factorx*(x2^(m+1)-x1^(m+1));
        else if i>=nlayerx+ngmx+1
                 x1=(i-nlayerx-ngmx-1)*dx;
                 x2=(i-nlayerx-ngmx)*dx;
                 Ox=epsreff*Om*factorx*(x2^(m+1)-x1^(m+1));
            else 
                Ox=0;
            end
        end
        y1=(nlayery-j)*dy;
        y2=(nlayery+1-j)*dy;
        Oy=epsreff*Om*factory*(y2^(m+1)-y1^(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 back area (upper PML)            
 % coefficient for computing Ex in back
for i=1:ngx
    for j=nlayery+ngmy+1:ngy % PEC boundary when j==nny
        epsreff=(epsr(i,j)+epsr(i+1,j))/2; % effective epsr
        if i<=nlayerx
            x1=(nlayerx-i)*dx;
            x2=(nlayerx+1-i)*dx;
            Ox=epsreff*Om*factorx*(x2^(m+1)-x1^(m+1));
        else if i>=nlayerx+ngmx+1
                 x1=(i-nlayerx-ngmx-1)*dx;
                 x2=(i-nlayerx-ngmx)*dx;
                 Ox=epsreff*Om*factorx*(x2^(m+1)-x1^(m+1));
            else
                Ox=0;
            end
        end
        y1=(j-nlayery-ngmy-1.5)*dy;
        y2=(j-nlayery-ngmy-0.5)*dy;
        if j==nlayery+ngmy+1
            Oy=epsreff*Om*factory*((0.5*dy)^(m+1));  % special treating in the under boundary of back PML area
        else
            Oy=epsreff*Om*factory*(y2^(m+1)-y1^(m+1));
        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 back area
for i=2:nnx-1 %it's PEC boundary for i==1 & i==nnx ,all coefficients are zero,need not treating
    for j=nlayery+ngmy+1:ngy
        epsreff=(epsr(i,j)+epsr(i,j+1))/2; % effective epsr
        if i<=nlayerx
            x1=(nlayerx+0.5-i)*dx;
            x2=(nlayerx+1.5-i)*dx;
            Ox=epsreff*Om*factorx*(x2^(m+1)-x1^(m+1));
        else if i>=nlayerx+ngmx+2
                 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));
            else if i==nlayerx+1|i==nlayerx+ngmx+1  % special treating for Ox in two boundaries in back PML
                    Ox=epsreff*Om*factorx*((0.5*dx)^(m+1));
                else
                   Ox=0;
                end
            end
        end
        y1=(j-nlayery-ngmy-1)*dy;
        y2=(j-nlayery-ngmy)*dy;
        Oy=epsreff*Om*factory*(y2^(m+1)-y1^(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 back
for i=1:ngx
    for j=nlayery+ngmy+1:ngy
        epsreff=0.25*(epsr(i,j)+epsr(i,j+1)+epsr(i+1,j)+epsr(i+1,j+1));
        if i<=nlayerx
            x1=(nlayerx-i)*dx;
            x2=(nlayerx+1-i)*dx;
            Ox=epsreff*Om*factorx*(x2^(m+1)-x1^(m+1));
        else if i>=nlayerx+ngmx+1
                 x1=(i-nlayerx-ngmx-1)*dx;
                 x2=(i-nlayerx-ngmx)*dx;
                 Ox=epsreff*Om*factorx*(x2^(m+1)-x1^(m+1));
            else 
                Ox=0;
            end
        end
        y1=(j-nlayery-ngmy-1)*dy;
        y2=(j-nlayery-ngmy)*dy;
        Oy=epsreff*Om*factory*(y2^(m+1)-y1^(m+1));
        bzbz(i,j)=(1/dt-Oy/2/eps0)/(1/dt+Oy/2/eps0);