fdtd3d.m

3-D FDTD code with PML absorbing boundary conditions

ca1=exp(-sigmay*dt/(epsz*eps(1)));
cb1=(1.0-ca1)/(sigmay*ds);
caex(1:ie,1,1:kb)=ca1;
cbex(1:ie,1,1:kb)=cb1;
caez(1:ib,1,1:ke)=ca1;
cbez(1:ib,1,1:ke)=cb1;
caexybcd(1:ie,1,1:kebc)=ca1;
cbexybcd(1:ie,1,1:kebc)=cb1;
caezybcd(1:ib,1,1:kebc)=ca1;
cbezybcd(1:ib,1,1:kebc)=cb1;
caexybct(1:ie,1,1:kebc)=ca1;
cbexybct(1:ie,1,1:kebc)=cb1;
caezybct(1:ib,1,1:kebc)=ca1;
cbezybct(1:ib,1,1:kebc)=cb1;
caexybcl(1:iebc,1,2:kefbc)=ca1;
cbexybcl(1:iebc,1,2:kefbc)=cb1;
caezybcl(1:iebc,1,1:kefbc)=ca1;
cbezybcl(1:iebc,1,1:kefbc)=cb1;
caexybcr(1:iebc,1,2:kefbc)=ca1;
cbexybcr(1:iebc,1,2:kefbc)=cb1;
caezybcr(1:iebc,1,1:kefbc)=ca1;
cbezybcr(1:iebc,1,1:kefbc)=cb1;


for j=1:jebc                              
  y1=(jebc-j+1)*ds;
  y2=(jebc-j)*ds;
  sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
  sigmays=sigmay*(muz/(epsz*eps(1)));
  da1=exp(-sigmays*dt/muz);
  db1=(1.0-da1)/(sigmays*ds);
  dahxybcf(2:iefbc,j,1:kefbc)=da1;
  dbhxybcf(2:iefbc,j,1:kefbc)=db1;
  dahzybcf(1:iefbc,j,2:kefbc)=da1;
  dbhzybcf(1:iefbc,j,2:kefbc)=db1;

  caeyxbcf(2:iefbc,j,2:kefbc)=ca(1);                   
  cbeyxbcf(2:iefbc,j,2:kefbc)=cb(1);
  caeyzbcf(2:iefbc,j,2:kefbc)=ca(1);                   
  cbeyzbcf(2:iefbc,j,2:kefbc)=cb(1);
  dahxzbcf(1:ibfbc,j,1:kefbc)=da(1);
  dbhxzbcf(1:ibfbc,j,1:kefbc)=db(1);   
  dahzxbcf(1:iefbc,j,1:kbfbc)=da(1);
  dbhzxbcf(1:iefbc,j,1:kbfbc)=db(1);
end



%     BACK region 
% back face
caexybcb(1:iefbc,jbbc,1:kbfbc)=1.0;
cbexybcb(1:iefbc,jbbc,1:kbfbc)=0.0;
caexzbcb(1:iefbc,jbbc,1:kbfbc)=1.0;
cbexzbcb(1:iefbc,jbbc,1:kbfbc)=0.0;
caezxbcb(1:ibfbc,jbbc,1:kefbc)=1.0;
cbezxbcb(1:ibfbc,jbbc,1:kefbc)=0.0;
caezybcb(1:ibfbc,jbbc,1:kefbc)=1.0;
cbezybcb(1:ibfbc,jbbc,1:kefbc)=0.0;

dahyxbcb(1:iefbc,jbbc,1:kefbc)=1.0;
dbhyxbcb(1:iefbc,jbbc,1:kefbc)=0.0;
dahyzbcb(1:iefbc,jbbc,1:kefbc)=1.0;
dbhyzbcb(1:iefbc,jbbc,1:kefbc)=0.0;
% left face
caeyxbcb(1,1:jebc,1:kbfbc)=1.0;
cbeyxbcb(1,1:jebc,1:kbfbc)=0.0;
caeyzbcb(1,1:jebc,1:kbfbc)=1.0;
cbeyzbcb(1,1:jebc,1:kbfbc)=0.0;
caezxbcb(1,1:jbbc,1:kefbc)=1.0;
cbezxbcb(1,1:jbbc,1:kefbc)=0.0;
caezybcb(1,1:jbbc,1:kefbc)=1.0;
cbezybcb(1,1:jbbc,1:kefbc)=0.0;

dahxybcb(1,1:jebc,1:kebc)=1.0;
dbhxybcb(1,1:jebc,1:kebc)=0.0;
dahxzbcb(1,1:jebc,1:kebc)=1.0;
dbhxzbcb(1,1:jebc,1:kebc)=0.0;
% right face
caeyxbcb(ibfbc,1:jebc,1:kbfbc)=1.0;
cbeyxbcb(ibfbc,1:jebc,1:kbfbc)=0.0;
caeyzbcb(ibfbc,1:jebc,1:kbfbc)=1.0;
cbeyzbcb(ibfbc,1:jebc,1:kbfbc)=0.0;
caezxbcb(ibfbc,1:jbbc,1:kefbc)=1.0;
cbezxbcb(ibfbc,1:jbbc,1:kefbc)=0.0;
caezybcb(ibfbc,1:jbbc,1:kefbc)=1.0;
cbezybcb(ibfbc,1:jbbc,1:kefbc)=0.0;

dahxybcb(ibfbc,1:jebc,1:kebc)=1.0;
dbhxybcb(ibfbc,1:jebc,1:kebc)=0.0;
dahxzbcb(ibfbc,1:jebc,1:kebc)=1.0;
dbhxzbcb(ibfbc,1:jebc,1:kebc)=0.0;
% bottom face
caexybcb(1:iefbc,1:jbbc,1)=1.0;
cbexybcb(1:iefbc,1:jbbc,1)=0.0;
caexzbcb(1:iefbc,1:jbbc,1)=1.0;
cbexzbcb(1:iefbc,1:jbbc,1)=0.0;
caeyxbcb(1:ibfbc,1:jebc,1)=1.0;
cbeyxbcb(1:ibfbc,1:jebc,1)=0.0;
caeyzbcb(1:ibfbc,1:jebc,1)=1.0;
cbeyzbcb(1:ibfbc,1:jebc,1)=0.0;

dahzxbcb(1:iefbc,1:jebc,1)=1.0;
dbhzxbcb(1:iefbc,1:jebc,1)=0.0;
dahzybcb(1:iefbc,1:jebc,1)=1.0;
dbhzybcb(1:iefbc,1:jebc,1)=0.0;
% top face
caexybcb(1:iefbc,1:jbbc,kbfbc)=1.0;
cbexybcb(1:iefbc,1:jbbc,kbfbc)=0.0;
caexzbcb(1:iefbc,1:jbbc,kbfbc)=1.0;
cbexzbcb(1:iefbc,1:jbbc,kbfbc)=0.0;
caeyxbcb(1:ibfbc,1:jebc,kbfbc)=1.0;
cbeyxbcb(1:ibfbc,1:jebc,kbfbc)=0.0;
caeyzbcb(1:ibfbc,1:jebc,kbfbc)=1.0;
cbeyzbcb(1:ibfbc,1:jebc,kbfbc)=0.0;

dahzxbcb(1:iefbc,1:jebc,kbfbc)=1.0;
dbhzxbcb(1:iefbc,1:jebc,kbfbc)=0.0;
dahzybcb(1:iefbc,1:jebc,kbfbc)=1.0;
dbhzybcb(1:iefbc,1:jebc,kbfbc)=0.0;


for j=2:jebc
  y1=(j-0.5)*ds;
  y2=(j-1.5)*ds;
  sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
  ca1=exp(-sigmay*dt/(epsz*eps(1)));
  cb1=(1.0-ca1)/(sigmay*ds);
  caexybcb(1:iefbc,j,2:kefbc)=ca1;
  cbexybcb(1:iefbc,j,2:kefbc)=cb1;
  caezybcb(2:iefbc,j,1:kefbc)=ca1;
  cbezybcb(2:iefbc,j,1:kefbc)=cb1;
  
  caexzbcb(2:iefbc,j,2:kefbc)=ca(1);
  cbexzbcb(2:iefbc,j,2:kefbc)=cb(1); 
  caezxbcb(2:iefbc,j,1:kefbc)=ca(1);
  cbezxbcb(2:iefbc,j,1:kefbc)=cb(1);  
end

sigmay = bcfactor*(0.5*ds)^(orderbc+1);
ca1=exp(-sigmay*dt/(epsz*eps(1)));
cb1=(1.0-ca1)/(sigmay*ds);
caex(1:ie,jb,1:kb)=ca1;
cbex(1:ie,jb,1:kb)=cb1;
caez(1:ib,jb,1:ke)=ca1;
cbez(1:ib,jb,1:ke)=cb1;
caexybcd(1:ie,jb,1:kebc)=ca1;
cbexybcd(1:ie,jb,1:kebc)=cb1;
caezybcd(1:ib,jb,1:kebc)=ca1;
cbezybcd(1:ib,jb,1:kebc)=cb1;
caexybct(1:ie,jb,1:kebc)=ca1;
cbexybct(1:ie,jb,1:kebc)=cb1;
caezybct(1:ib,jb,1:kebc)=ca1;
cbezybct(1:ib,jb,1:kebc)=cb1;
caexybcl(1:iebc,jb,2:kefbc)=ca1;
cbexybcl(1:iebc,jb,2:kefbc)=cb1;
caezybcl(2:iebc,jb,1:kefbc)=ca1;
cbezybcl(2:iebc,jb,1:kefbc)=cb1;
caexybcr(1:iebc,jb,2:kefbc)=ca1;
cbexybcr(1:iebc,jb,2:kefbc)=cb1;
caezybcr(1:iebc,jb,1:kefbc)=ca1;
cbezybcr(1:iebc,jb,1:kefbc)=cb1;

for j=1:jebc
  y1=j*ds;
  y2=(j-1)*ds;
  sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
  sigmays=sigmay*(muz/(epsz*eps(1)));
  da1=exp(-sigmays*dt/muz);
  db1=(1.0-da1)/(sigmays*ds);
  dahxybcb(1:ibfbc,j,1:kefbc)=da1;
  dbhxybcb(1:ibfbc,j,1:kefbc)=db1;
  dahzybcb(1:iefbc,j,1:kbfbc)=da1;
  dbhzybcb(1:iefbc,j,1:kbfbc)=db1;    
   
  caeyxbcb(2:iefbc,j,2:kefbc)=ca(1);
  cbeyxbcb(2:iefbc,j,2:kefbc)=cb(1);
  caeyzbcb(2:iefbc,j,2:kefbc)=ca(1);
  cbeyzbcb(2:iefbc,j,2:kefbc)=cb(1);   
  dahxzbcb(1:ibfbc,j,1:kefbc)=da(1);
  dbhxzbcb(1:ibfbc,j,1:kefbc)=db(1);
  dahyzbcb(1:iefbc,j,1:kefbc)=da(1);
  dbhyzbcb(1:iefbc,j,1:kefbc)=db(1); 
  dahyxbcb(1:iefbc,j,1:kefbc)=da(1);
  dbhyxbcb(1:iefbc,j,1:kefbc)=db(1); 
  dahzxbcb(1:iefbc,j,1:kbfbc)=da(1);
  dbhzxbcb(1:iefbc,j,1:kbfbc)=db(1);
end



%     LEFT region 
% left face
caeyzbcl(1,1:je,1:kbfbc)=1.0;
cbeyzbcl(1,1:je,1:kbfbc)=0.0;
caeyxbcl(1,1:je,1:kbfbc)=1.0;
cbeyxbcl(1,1:je,1:kbfbc)=0.0;
caezxbcl(1,1:jb,1:kefbc)=1.0;
cbezxbcl(1,1:jb,1:kefbc)=0.0;
caezybcl(1,1:jb,1:kefbc)=1.0;
cbezybcl(1,1:jb,1:kefbc)=0.0;

dahxybcl(1,1:je,1:kefbc)=1.0;
dbhxybcl(1,1:je,1:kefbc)=0.0;
dahxzbcl(1,1:je,1:kefbc)=1.0;
dbhxzbcl(1,1:je,1:kefbc)=0.0;
% top face
caexybcl(1:iebc,1:jb,kbfbc)=1.0;
cbexybcl(1:iebc,1:jb,kbfbc)=0.0;
caexzbcl(1:iebc,1:jb,kbfbc)=1.0;
cbexzbcl(1:iebc,1:jb,kbfbc)=0.0;
caeyxbcl(1:iebc,1:je,kbfbc)=1.0;
cbeyxbcl(1:iebc,1:je,kbfbc)=0.0;
caeyzbcl(1:iebc,1:je,kbfbc)=1.0;
cbeyzbcl(1:iebc,1:je,kbfbc)=0.0;

dahzxbcl(1:iebc,1:je,kbfbc)=1.0;
dbhzxbcl(1:iebc,1:je,kbfbc)=0.0;
dahzybcl(1:iebc,1:je,kbfbc)=1.0;
dbhzybcl(1:iebc,1:je,kbfbc)=0.0;
% bottom face
caexybcl(1:iebc,1:jb,1)=1.0;
cbexybcl(1:iebc,1:jb,1)=0.0;
caexzbcl(1:iebc,1:jb,1)=1.0;
cbexzbcl(1:iebc,1:jb,1)=0.0;
caeyxbcl(1:iebc,1:je,1)=1.0;
cbeyxbcl(1:iebc,1:je,1)=0.0;
caeyzbcl(1:iebc,1:je,1)=1.0;
cbeyzbcl(1:iebc,1:je,1)=0.0;

dahzxbcl(1:iebc,1:je,1)=1.0;
dbhzxbcl(1:iebc,1:je,1)=0.0;
dahzybcl(1:iebc,1:je,1)=1.0;
dbhzybcl(1:iebc,1:je,1)=0.0;

for i=2:iebc
  x1=(iebc-i+1.5)*ds;
  x2=(iebc-i+0.5)*ds;
  sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
  ca1=exp(-sigmax*dt/(epsz*eps(1)));
  cb1=(1.0-ca1)/(sigmax*ds);
  caeyxbcl(i,1:je,2:kefbc)=ca1;
  cbeyxbcl(i,1:je,2:kefbc)=cb1;
  caezxbcl(i,1:jb,1:kefbc)=ca1;
  cbezxbcl(i,1:jb,1:kefbc)=cb1; 
  caeyzbcl(i,1:je,2:kefbc)=ca(1);
  cbeyzbcl(i,1:je,2:kefbc)=cb(1); 
  caezybcl(i,2:je,1:kefbc)=ca(1);
  cbezybcl(i,2:je,1:kefbc)=cb(1);
  caeyxbcf(i,1:jebc,2:kefbc)=ca1;
  cbeyxbcf(i,1:jebc,2:kefbc)=cb1;
  caezxbcf(i,2:jebc,1:kefbc)=ca1;
  cbezxbcf(i,2:jebc,1:kefbc)=cb1; 
  caeyxbcb(i,1:jebc,2:kefbc)=ca1;
  cbeyxbcb(i,1:jebc,2:kefbc)=cb1;
  caezxbcb(i,1:jebc,1:kefbc)=ca1;
  cbezxbcb(i,1:jebc,1:kefbc)=cb1;  
  dahxybcl(i,1:je,1:kefbc)=da(1);
  dbhxybcl(i,1:je,1:kefbc)=db(1);
  dahxzbcl(i,1:je,1:kefbc)=da(1);
  dbhxzbcl(i,1:je,1:kefbc)=db(1);  
end
 
sigmax=bcfactor*(0.5*ds)^(orderbc+1);
ca1=exp(-sigmax*dt/(epsz*eps(1)));
cb1=(1.0-ca1)/(sigmax*ds);
caey(1,1:je,1:kb)=ca1;
cbey(1,1:je,1:kb)=cb1;
caez(1,1:jb,1:ke)=ca1;
cbez(1,1:jb,1:ke)=cb1;
caeyxbcd(1,1:je,2:kebc)=ca1;
cbeyxbcd(1,1:je,2:kebc)=cb1;
caezxbcd(1,1:jb,1:kebc)=ca1;
cbezxbcd(1,1:jb,1:kebc)=cb1;
caeyxbct(1,1:je,1:kebc)=ca1;
cbeyxbct(1,1:je,1:kebc)=cb1;
caezxbct(1,1:jb,1:kebc)=ca1;
cbezxbct(1,1:jb,1:kebc)=cb1;
caeyxbcf(ibbc,1:jebc,2:kefbc)=ca1;
cbeyxbcf(ibbc,1:jebc,2:kefbc)=cb1;
caezxbcf(ibbc,1:jebc,1:kefbc)=ca1;
cbezxbcf(ibbc,1:jebc,1:kefbc)=cb1;
caeyxbcb(ibbc,1:jebc,2:kefbc)=ca1;
cbeyxbcb(ibbc,1:jebc,2:kefbc)=cb1;
caezxbcb(ibbc,2:jebc,1:kefbc)=ca1;
cbezxbcb(ibbc,2:jebc,1:kefbc)=cb1;

 
for i=1:iebc
  x1=(iebc-i+1)*ds;
  x2=(iebc-i)*ds;
  sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
  sigmaxs=sigmax*(muz/(epsz*eps(1)));
  da1=exp(-sigmaxs*dt/muz);
  db1=(1.0-da1)/(sigmaxs*ds);
  dahyxbcl(i,1:jb,1:kefbc)=da1;
  dbhyxbcl(i,1:jb,1:kefbc)=db1;
  dahzxbcl(i,1:je,2:kefbc)=da1;
  dbhzxbcl(i,1:je,2:kefbc)=db1;  
  dahyxbcf(i,1:jebc,1:kefbc)=da1;
  dbhyxbcf(i,1:jebc,1:kefbc)=db1;  
  dahzxbcf(i,1:jebc,2:kefbc)=da1;
  dbhzxbcf(i,1:jebc,2:kefbc)=db1;
  dahyxbcb(i,1:jebc,1:kefbc)=da1;
  dbhyxbcb(i,1:jebc,1:kefbc)=db1;  
  dahzxbcb(i,1:jebc,2:kefbc)=da1;
  dbhzxbcb(i,1:jebc,2:kefbc)=db1;  

  caexybcl(i,2:je,2:kefbc)=ca(1);
  cbexybcl(i,2:je,2:kefbc)=cb(1);
  caexzbcl(i,1:jb,2:kefbc)=ca(1);
  cbexzbcl(i,1:jb,2:kefbc)=cb(1); 
  dahyzbcl(i,1:jb,1:kefbc)=da(1);
  dbhyzbcl(i,1:jb,1:kefbc)=db(1); 
  dahzybcl(i,1:je,2:kefbc)=da(1);
  dbhzybcl(i,1:je,2:kefbc)=db(1);
end


%     RIGHT region 
% right face
caeyxbcr(ibbc,1:je,1:kbfbc)=1.0;
cbeyxbcr(ibbc,1:je,1:kbfbc)=0.0;
caeyzbcr(ibbc,1:je,1:kbfbc)=1.0;
cbeyzbcr(ibbc,1:je,1:kbfbc)=0.0;
caezxbcr(ibbc,1:jb,1:kefbc)=1.0;
cbezxbcr(ibbc,1:jb,1:kefbc)=0.0;
caezybcr(ibbc,1:jb,1:kefbc)=1.0;
cbezybcr(ibbc,1:jb,1:kefbc)=0.0;

dahxybcr(ibbc,1:je,1:kefbc)=1.0;
dbhxybcr(ibbc,1:je,1:kefbc)=0.0;
dahxzbcr(ibbc,1:je,1:kefbc)=1.0;
dbhxzbcr(ibbc,1:je,1:kefbc)=0.0;
% top face
caexybcr(1:iebc,1:jb,kbfbc)=1.0;
cbexybcr(1:iebc,1:jb,kbfbc)=0.0;
caexzbcr(1:iebc,1:jb,kbfbc)=1.0;
cbexzbcr(1:iebc,1:jb,kbfbc)=0.0;
caeyxbcr(1:ibbc,1:je,kbfbc)=1.0;
cbeyxbcr(1:ibbc,1:je,kbfbc)=0.0;
caeyzbcr(1:ibbc,1:je,kbfbc)=1.0;
cbeyzbcr(1:ibbc,1:je,kbfbc)=0.0;

dahzxbcr(1:iebc,1:je,kbfbc)=1.0;
dbhzxbcr(1:iebc,1:je,kbfbc)=0.0;
dahzybcr(1:iebc,1:je,kbfbc)=1.0;
dbhzybcr(1:iebc,1:je,kbfbc)=0.0;
%bottom face
caexybcr(1:iebc,1:jb,1)=1.0;
cbexybcr(1:iebc,1:jb,1)=0.0;
caexzbcr(1:iebc,1:jb,1)=1.0;
cbexzbcr(1:iebc,1:jb,1)=0.0;
caeyxbcr(1:ibbc,1:je,1)=1.0;
cbeyxbcr(1:ibbc,1:je,1)=0.0;
caeyzbcr(1:ibbc,1:je,1)=1.0;
cbeyzbcr(1:ibbc,1:je,1)=0.0;