%***********************************************************************
%     3-D FDTD code with PML absorbing boundary conditions
%***********************************************************************

%    program auther:   Lijijun 
%                      Department of Physics , Ynagtze University
%                      lijijun@163.net
%                   
%    Date of this version: August 2006                  
%
%***********************************************************************
clear
clc
%***********************************************************************
%     Fundamental constants
%***********************************************************************
nm=1e-9;
cc=2.99792458e8;            %speed of light in free space
muz=4.0*pi*1.0e-7;          %permeability of free space
epsz=1.0/(cc*cc*muz);       %permittivity of free space

lambda=600*nm;
freq=cc/lambda;
omega=2.0*pi*freq;    



%***********************************************************************
%     Grid parameters
%***********************************************************************

ie=30;           %number of grid cells in x-direction
je=40;            %number of grid cells in y-direction
ke=60;             %number of grid cells in z-direction

ib=ie+1;
jb=je+1;
kb=ke+1;

is=15;            %location of  hard source in x axis
js=je/2;          %location of  hard source in y axis
ks=ke/2;          %location of  hard source in z axis

ds=3.0*nm;        %space increment of square lattice
dt=ds/(2.0*cc);   %time step

nmax=300;         %total number of time steps

iebc=8;           %thickness of left and right PML region
jebc=8;           %thickness of front and back PML region
kebc=8;           %thickness of bottom and top PML region
rmax=0.00001;
orderbc=2;
ibbc=iebc+1;
jbbc=jebc+1;
kbbc=kebc+1;

iefbc=ie+2*iebc;
jefbc=je+2*jebc;
kefbc=ke+2*kebc;
ibfbc=iefbc+1;
jbfbc=jefbc+1;
kbfbc=kefbc+1;

%***********************************************************************
%     Material parameters
%***********************************************************************

media=2;

eps=[1.0 1.0];
sig=[0.0 1.0e+7];
mur=[1.0 1.0];
sim=[0.0 0.0];

%***********************************************************************
%     Wave excitation
%***********************************************************************
for n=1:nmax

if n<pi/omega/dt
   U=.5*(1-cos(pi*n/(pi/omega/dt)));
else 
    U=1;
end   %設(shè)置升余弦開(kāi)關(guān)函數(shù),減少波動(dòng)。
source(n)=exp(-i*n*omega*dt)*U;
end 

%***********************************************************************
%     Field arrays
%***********************************************************************

ex=zeros(ie,jb,kb);           %fields in main grid 
ey=zeros(ib,je,kb);
ez=zeros(ib,jb,ke);

hx=zeros(ib,je,ke);
hy=zeros(ie,jb,ke);
hz=zeros(ie,je,kb);


exybcf=zeros(iefbc,jebc,kbfbc);   %fields in front PML region
exzbcf=zeros(iefbc,jebc,kbfbc);
eyzbcf=zeros(ibfbc,jebc,kbfbc);
eyxbcf=zeros(ibfbc,jebc,kbfbc);
ezxbcf=zeros(ibfbc,jebc,kefbc);
ezybcf=zeros(ibfbc,jebc,kefbc);

hxybcf=zeros(ibfbc,jebc,kefbc);
hxzbcf=zeros(ibfbc,jebc,kefbc);
hyzbcf=zeros(iefbc,jebc,kefbc);
hyxbcf=zeros(iefbc,jebc,kefbc);
hzxbcf=zeros(iefbc,jebc,kbfbc);
hzybcf=zeros(iefbc,jebc,kbfbc);


exybcb=zeros(iefbc,jbbc,kbfbc);   %fields in back PML region
exzbcb=zeros(iefbc,jbbc,kbfbc);
eyzbcb=zeros(ibfbc,jebc,kbfbc);
eyxbcb=zeros(ibfbc,jebc,kbfbc);
ezxbcb=zeros(ibfbc,jbbc,kefbc);
ezybcb=zeros(ibfbc,jbbc,kefbc);

hxybcb=zeros(ibfbc,jebc,kefbc);
hxzbcb=zeros(ibfbc,jebc,kefbc);
hyzbcb=zeros(iefbc,jbbc,kefbc);
hyxbcb=zeros(iefbc,jbbc,kefbc);
hzxbcb=zeros(iefbc,jebc,kbfbc);
hzybcb=zeros(iefbc,jebc,kbfbc);


exybcl=zeros(iebc,jb,kbfbc);      %fields in left PML region
exzbcl=zeros(iebc,jb,kbfbc);
eyzbcl=zeros(iebc,je,kbfbc);
eyxbcl=zeros(iebc,je,kbfbc);
ezxbcl=zeros(iebc,jb,kefbc);
ezybcl=zeros(iebc,jb,kefbc);

hxybcl=zeros(iebc,je,kefbc);
hxzbcl=zeros(iebc,je,kefbc);
hyzbcl=zeros(iebc,jb,kefbc);
hyxbcl=zeros(iebc,jb,kefbc);
hzxbcl=zeros(iebc,je,kbfbc);
hzybcl=zeros(iebc,je,kbfbc);


exybcr=zeros(iebc,jb,kbfbc);      %fields in right PML region
exzbcr=zeros(iebc,jb,kbfbc);
eyzbcr=zeros(ibbc,je,kbfbc);
eyxbcr=zeros(ibbc,je,kbfbc);
ezxbcr=zeros(ibbc,jb,kefbc);
ezybcr=zeros(ibbc,jb,kefbc);

hxybcr=zeros(ibbc,je,kefbc);
hxzbcr=zeros(ibbc,je,kefbc);
hyzbcr=zeros(iebc,jb,kefbc);
hyxbcr=zeros(iebc,jb,kefbc);
hzxbcr=zeros(iebc,je,kbfbc);
hzybcr=zeros(iebc,je,kbfbc);


exybcd=zeros(ie,jb,kebc);           %fields in bottom PML region
exzbcd=zeros(ie,jb,kebc);
eyzbcd=zeros(ib,je,kebc);
eyxbcd=zeros(ib,je,kebc);
ezxbcd=zeros(ib,jb,kebc);
ezybcd=zeros(ib,jb,kebc);

hxybcd=zeros(ib,je,kebc);
hxzbcd=zeros(ib,je,kebc);
hyzbcd=zeros(ie,jb,kebc);
hyxbcd=zeros(ie,jb,kebc);
hzybcd=zeros(ie,je,kebc);
hzxbcd=zeros(ie,je,kebc);


exybct=zeros(ie,jb,kbbc);          %fields in top PML region
exzbct=zeros(ie,jb,kbbc);
eyzbct=zeros(ib,je,kbbc); 
eyxbct=zeros(ib,je,kbbc);
ezxbct=zeros(ib,jb,kebc);
ezybct=zeros(ib,jb,kebc);

hxybct=zeros(ib,je,kebc);
hxzbct=zeros(ib,je,kebc);
hyzbct=zeros(ie,jb,kebc);
hyxbct=zeros(ie,jb,kebc);
hzxbct=zeros(ie,je,kbbc);
hzybct=zeros(ie,je,kbbc);


%***********************************************************************
%     Updating coefficients
%***********************************************************************

for i=1:media
  eaf=dt*sig(i)/(2.0*epsz*eps(i));
  ca(i)=(1.0-eaf)/(1.0+eaf);
  cb(i)=dt/epsz/eps(i)/ds/(1.0+eaf);
  haf=dt*sim(i)/(2.0*muz*mur(i));
  da(i)=(1.0-haf)/(1.0+haf);
  db(i)=dt/muz/mur(i)/ds/(1.0+haf);
end

%***********************************************************************
%     Geometry specification (main grid)
%***********************************************************************

%     Initialize entire main grid to free space

caex(1:ie,1:jb,1:kb)=ca(1);     
cbex(1:ie,1:jb,1:kb)=cb(1);

caey(1:ib,1:je,1:kb)=ca(1);
cbey(1:ib,1:je,1:kb)=cb(1);

caez(1:ib,1:jb,1:ke)=ca(1);
cbez(1:ib,1:jb,1:ke)=cb(1);


dahx(1:ib,1:je,1:ke)=da(1);
dbhx(1:ib,1:je,1:ke)=db(1);

dahy(1:ie,1:jb,1:ke)=da(1);
dbhy(1:ie,1:jb,1:ke)=db(1);

dahz(1:ie,1:je,1:kb)=da(1);
dbhz(1:ie,1:je,1:kb)=db(1);

% %     Add metal cylinder
% 
% diam=20;          % diameter of cylinder: 6 cm
% rad=diam/2.0;     % radius of cylinder: 3 cm
% icenter=4*ie/5;   % i-coordinate of cylinder's center
% jcenter=je/2;     % j-coordinate of cylinder's center
% 
% for i=1:ie
% for j=1:je
% for k=10:20
%   dist2=(i+0.5-icenter)^2 + (j-jcenter)^2;
%   if dist2 <= rad^2 
%      caex(i,j,k)=ca(2);
%      cbex(i,j,k)=cb(2);
%   end
%   dist2=(i-icenter)^2 + (j+0.5-jcenter)^2;
%   if dist2 <= rad^2 
%      caey(i,j,k)=ca(2);
%      cbey(i,j,k)=cb(2);
%   end
% end
% end

%***********************************************************************
%     Fill the PML regions
%***********************************************************************

delbc=iebc*ds;
sigmam=-log(rmax/100.0)*epsz*cc*(orderbc+1)/(2*delbc);
bcfactor=eps(1)*sigmam/(ds*(delbc^orderbc)*(orderbc+1));


%     FRONT region 

% front face
caexybcf(1:iefbc,1,1:kbfbc)=1.0;
cbexybcf(1:iefbc,1,1:kbfbc)=0.0;
caexzbcf(1:iefbc,1,1:kbfbc)=1.0;
cbexzbcf(1:iefbc,1,1:kbfbc)=0.0;
caezxbcf(1:ibfbc,1,1:kefbc)=1.0;
cbezxbcf(1:ibfbc,1,1:kefbc)=0.0;
caezybcf(1:ibfbc,1,1:kefbc)=1.0;
cbezybcf(1:ibfbc,1,1:kefbc)=0.0;

dahyxbcf(1:iefbc,1,1:kefbc)=1.0;
dbhyxbcf(1:iefbc,1,1:kefbc)=0.0;
dahyzbcf(1:iefbc,1,1:kefbc)=1.0;
dbhyzbcf(1:iefbc,1,1:kefbc)=0.0;
% left face
caeyxbcf(1,1:jebc,1:kbfbc)=1.0;
cbeyxbcf(1,1:jebc,1:kbfbc)=0.0;
caeyzbcf(1,1:jebc,1:kbfbc)=1.0;
cbeyzbcf(1,1:jebc,1:kbfbc)=0.0;
caezxbcf(1,1:jebc,1:kefbc)=1.0;
cbezxbcf(1,1:jebc,1:kefbc)=0.0;
caezybcf(1,1:jebc,1:kefbc)=1.0;
cbezybcf(1,1:jebc,1:kefbc)=0.0;

dahxybcf(1,1:jebc,1:kebc)=1.0;
dbhxybcf(1,1:jebc,1:kebc)=0.0;
dahxzbcf(1,1:jebc,1:kebc)=1.0;
dbhxzbcf(1,1:jebc,1:kebc)=0.0;
% right face
caeyxbcf(ibfbc,1:jebc,1:kbfbc)=1.0;
cbeyxbcf(ibfbc,1:jebc,1:kbfbc)=0.0;
caeyzbcf(ibfbc,1:jebc,1:kbfbc)=1.0;
cbeyzbcf(ibfbc,1:jebc,1:kbfbc)=0.0;
caezxbcf(ibfbc,1:jebc,1:kefbc)=1.0;
cbezxbcf(ibfbc,1:jebc,1:kefbc)=0.0;
caezybcf(ibfbc,1:jebc,1:kefbc)=1.0;
cbezybcf(ibfbc,1:jebc,1:kefbc)=0.0;

dahxybcf(ibfbc,1:jebc,1:kebc)=1.0;
dbhxybcf(ibfbc,1:jebc,1:kebc)=0.0;
dahxzbcf(ibfbc,1:jebc,1:kebc)=1.0;
dbhxzbcf(ibfbc,1:jebc,1:kebc)=0.0;
% bottom face
caexybcf(1:iefbc,1:jebc,1)=1.0;
cbexybcf(1:iefbc,1:jebc,1)=0.0;
caexzbcf(1:iefbc,1:jebc,1)=1.0;
cbexzbcf(1:iefbc,1:jebc,1)=0.0;
caeyxbcf(1:ibfbc,1:jebc,1)=1.0;
cbeyxbcf(1:ibfbc,1:jebc,1)=0.0;
caeyzbcf(1:ibfbc,1:jebc,1)=1.0;
cbeyzbcf(1:ibfbc,1:jebc,1)=0.0;

dahzxbcf(1:iefbc,1:jebc,1)=1.0;
dbhzxbcf(1:iefbc,1:jebc,1)=0.0;
dahzybcf(1:iefbc,1:jebc,1)=1.0;
dbhzybcf(1:iefbc,1:jebc,1)=0.0;
% top face
caexybcf(1:iefbc,1:jebc,kbfbc)=1.0;
cbexybcf(1:iefbc,1:jebc,kbfbc)=0.0;
caexzbcf(1:iefbc,1:jebc,kbfbc)=1.0;
cbexzbcf(1:iefbc,1:jebc,kbfbc)=0.0;
caeyxbcf(1:ibfbc,1:jebc,kbfbc)=1.0;
cbeyxbcf(1:ibfbc,1:jebc,kbfbc)=0.0;
caeyzbcf(1:ibfbc,1:jebc,kbfbc)=1.0;
cbeyzbcf(1:ibfbc,1:jebc,kbfbc)=0.0;

dahzxbcf(1:iefbc,1:jebc,kbfbc)=1.0;
dbhzxbcf(1:iefbc,1:jebc,kbfbc)=0.0;
dahzybcf(1:iefbc,1:jebc,kbfbc)=1.0;
dbhzybcf(1:iefbc,1:jebc,kbfbc)=0.0;
% $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
for j=2:jebc                               
  y1=(jebc-j+1.5)*ds;
  y2=(jebc-j+0.5)*ds;
  sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
  ca1=exp(-sigmay*dt/(epsz*eps(1)));
  cb1=(1.0-ca1)/(sigmay*ds);
  caexybcf(1:iefbc,j,2:kefbc)=ca1;
  cbexybcf(1:iefbc,j,2:kefbc)=cb1;
  caezybcf(2:iefbc,j,1:kefbc)=ca1;
  cbezybcf(2:iefbc,j,1:kefbc)=cb1;  
  
  caexzbcf(1:iefbc,j,2:kefbc)=ca(1);                  
  cbexzbcf(1:iefbc,j,2:kefbc)=cb(1);
  caezxbcf(2:iefbc,j,1:kefbc)=ca(1);                   
  cbezxbcf(2:iefbc,j,1:kefbc)=cb(1);
  dahyzbcf(1:iefbc,j,1:kefbc)=da(1);
  dbhyzbcf(1:iefbc,j,1:kefbc)=db(1);
  dahyxbcf(1:iefbc,j,1:kefbc)=da(1);
  dbhyxbcf(1:iefbc,j,1:kefbc)=db(1);   
end

sigmay = bcfactor*(0.5*ds)^(orderbc+1);