fdtd3d.m

3-D FDTD code with PML absorbing boundary conditions

dahzxbcr(1:iebc,1:je,1)=1.0;
dbhzxbcr(1:iebc,1:je,1)=0.0;
dahzybcr(1:iebc,1:je,1)=1.0;
dbhzybcr(1:iebc,1:je,1)=0.0;


for i=2:iebc
  x1=(i-0.5)*ds;
  x2=(i-1.5)*ds;
  sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
  ca1=exp(-sigmax*dt/(epsz*eps(1)));
  cb1=(1.0-ca1)/(sigmax*ds); 
  caeyxbcr(i,1:je,2:kefbc)=ca1;
  cbeyxbcr(i,1:je,2:kefbc)=cb1;  
  caezxbcr(i,1:jb,1:kefbc)=ca1;
  cbezxbcr(i,1:jb,1:kefbc)=cb1; 
  caeyzbcr(i,1:je,2:kefbc)=ca(1);
  cbeyzbcr(i,1:je,2:kefbc)=cb(1); 
  caezybcr(i,2:je,1:kefbc)=ca(1);
  cbezybcr(i,2:je,1:kefbc)=cb(1);

  caeyxbcf(i+iebc+ie,1:jebc,2:kefbc)=ca1;
  cbeyxbcf(i+iebc+ie,1:jebc,2:kefbc)=cb1;
  caezxbcf(i+iebc+ie,2:jebc,1:kefbc)=ca1;
  cbezxbcf(i+iebc+ie,2:jebc,1:kefbc)=cb1;  
  caeyxbcb(i+iebc+ie,1:jebc,2:kefbc)=ca1;
  cbeyxbcb(i+iebc+ie,1:jebc,2:kefbc)=cb1;
  caezxbcb(i+iebc+ie,2:jebc,1:kefbc)=ca1;
  cbezxbcb(i+iebc+ie,2:jebc,1:kefbc)=cb1;  
  dahxybcr(i,1:je,1:kefbc)=da(1);
  dbhxybcr(i,1:je,1:kefbc)=db(1);  
  dahxzbcr(i,1:je,1:kefbc)=da(1);
  dbhxzbcr(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(ib,1:je,1:kb)=ca1;
cbey(ib,1:je,1:kb)=cb1;
caez(ib,1:jb,1:ke)=ca1;
cbez(ib,1:jb,1:ke)=cb1;
caeyxbcd(ib,1:je,2:kebc)=ca1;
cbeyxbcd(ib,1:je,2:kebc)=cb1;
caezxbcd(ib,1:jb,1:kebc)=ca1;
cbezxbcd(ib,1:jb,1:kebc)=cb1;
caeyxbct(ib,1:je,1:kebc)=ca1;
cbeyxbct(ib,1:je,1:kebc)=cb1;
caezxbct(ib,1:jb,1:kebc)=ca1;
cbezxbct(ib,1:jb,1:kebc)=cb1;
caeyxbcf(iebc+ib,1:jebc,2:kefbc)=ca1;
cbeyxbcf(iebc+ib,1:jebc,2:kefbc)=cb1;
caezxbcf(iebc+ib,1:jebc,1:kefbc)=ca1;
cbezxbcf(iebc+ib,1:jebc,1:kefbc)=cb1;
caeyxbcb(iebc+ib,1:jebc,2:kefbc)=ca1;
cbeyxbcb(iebc+ib,1:jebc,2:kefbc)=cb1;
caezxbcb(iebc+ib,2:jebc,1:kefbc)=ca1;
cbezxbcb(iebc+ib,2:jebc,1:kefbc)=cb1;


for i=1:iebc
  x1=i*ds;
  x2=(i-1)*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);  
  dahyxbcr(i,1:jb,1:kefbc)=da1;
  dbhyxbcr(i,1:jb,1:kefbc)=db1;  
  dahzxbcr(i,1:je,2:kefbc)=da1;
  dbhzxbcr(i,1:je,2:kefbc)=db1;  
  dahyxbcf(i+ie+iebc,1:jebc,1:kefbc)=da1;
  dbhyxbcf(i+ie+iebc,1:jebc,1:kefbc)=db1;  
  dahzxbcf(i+ie+iebc,1:jebc,2:kefbc)=da1;
  dbhzxbcf(i+ie+iebc,1:jebc,2:kefbc)=db1;  
  dahyxbcb(i+ie+iebc,1:jebc,1:kefbc)=da1;
  dbhyxbcb(i+ie+iebc,1:jebc,1:kefbc)=db1;  
  dahzxbcb(i+ie+iebc,1:jebc,2:kefbc)=da1;
  dbhzxbcb(i+ie+iebc,1:jebc,2:kefbc)=db1;   

  caeyzbcr(i,1:je,2:kefbc)=ca(1);
  cbeyzbcr(i,1:je,2:kefbc)=cb(1); 
  caezybcr(i,2:je,1:kefbc)=ca(1);
  cbezybcr(i,2:je,1:kefbc)=cb(1);    
  caexybcr(i,2:je,2:kefbc)=ca(1);
  cbexybcr(i,2:je,2:kefbc)=cb(1);
  caexzbcr(i,1:jb,2:kefbc)=ca(1);
  cbexzbcr(i,1:jb,2:kefbc)=cb(1); 
  dahyzbcr(i,1:jb,1:kefbc)=da(1);
  dbhyzbcr(i,1:jb,1:kefbc)=db(1);  
  dahzybcr(i,1:je,2:kefbc)=da(1);
  dbhzybcr(i,1:je,2:kefbc)=db(1);
end


%     BOTTOM region

caexybcd(1:ie,1:jb,1)=1.0;
cbexybcd(1:ie,1:jb,1)=0.0;
caexzbcd(1:ie,1:jb,1)=1.0;
cbexzbcd(1:ie,1:jb,1)=0.0;
caeyxbcd(1:ib,1:je,1)=1.0;
cbeyxbcd(1:ib,1:je,1)=0.0;
caeyzbcd(1:ib,1:je,1)=1.0;
cbeyzbcd(1:ib,1:je,1)=0.0;

dahzxbcd(1:ie,1:je,1)=1.0;
dbhzxbcd(1:ie,1:je,1)=0.0;
dahzybcd(1:ie,1:je,1)=1.0;
dbhzybcd(1:ie,1:je,1)=0.0;


for k=2:kebc                               % 与sigmaz有关的量
  z1=(kebc-k+1.5)*ds;
  z2=(kebc-k+0.5)*ds;
  sigmaz=bcfactor*(z1^(orderbc+1)-z2^(orderbc+1));
  ca1=exp(-sigmaz*dt/(epsz*eps(1)));
  cb1=(1.0-ca1)/(sigmaz*ds);  
  caexzbcd(1:ie,1:jb,k)=ca1;
  cbexzbcd(1:ie,1:jb,k)=cb1;  
  caeyzbcd(1:ib,1:je,k)=ca1;
  cbeyzbcd(1:ib,1:je,k)=cb1;  
  
  caexzbcf(1:iefbc,1:jebc,k)=ca1;
  cbexzbcf(1:iefbc,1:jebc,k)=cb1;  
  caeyzbcf(2:iefbc,1:jebc,k)=ca1;
  cbeyzbcf(2:iefbc,1:jebc,k)=cb1;  
  caexzbcb(1:iefbc,1:jebc,k)=ca1;
  cbexzbcb(1:iefbc,1:jebc,k)=cb1; 
  caeyzbcb(2:iefbc,1:jebc,k)=ca1;
  cbeyzbcb(2:iefbc,1:jebc,k)=cb1;  
  caexzbcl(1:iebc,1:jb,k)=ca1;
  cbexzbcl(1:iebc,1:jb,k)=cb1;  
  caeyzbcl(1:iebc,1:je,k)=ca1;
  cbeyzbcl(1:iebc,1:je,k)=cb1;  
  caexzbcr(1:iebc,1:jb,k)=ca1;
  cbexzbcr(1:iebc,1:jb,k)=cb1;  
  caeyzbcr(2:iebc,1:je,k)=ca1;
  cbeyzbcr(2:iebc,1:je,k)=cb1; 
  
  caexybcd(1:ie,2:je,k)=ca(1);
  cbexybcd(1:ie,2:je,k)=cb(1);  
  caeyxbcd(2:ie,1:je,k)=ca(1);                   
  cbeyxbcd(2:ie,1:je,k)=cb(1);  
  dahzxbcd(1:ie,1:je,k)=da(1);
  dbhzxbcd(1:ie,1:je,k)=db(1);  
  dahzybcd(1:ie,1:je,k)=da(1);
  dbhzybcd(1:ie,1:je,k)=db(1);   
end

sigmaz = bcfactor*(0.5*ds)^(orderbc+1); 
ca1=exp(-sigmaz*dt/(epsz*eps(1)));
cb1=(1.0-ca1)/(sigmaz*ds);
caex(1:ie,1:jb,1)=ca1;
cbex(1:ie,1:jb,1)=cb1;
caey(1:ib,1:je,1)=ca1;
cbey(1:ib,1:je,1)=cb1;
caexzbcf(1:iefbc,1:jebc,kbbc)=ca1;
cbexzbcf(1:iefbc,1:jebc,kbbc)=cb1;
caeyzbcf(2:iefbc,1:jebc,kbbc)=ca1;
cbeyzbcf(2:iefbc,1:jebc,kbbc)=cb1;
caexzbcb(1:iefbc,2:jebc,kbbc)=ca1;
cbexzbcb(1:iefbc,2:jebc,kbbc)=cb1;
caeyzbcb(2:iefbc,1:jebc,kbbc)=ca1;
cbeyzbcb(2:iefbc,1:jebc,kbbc)=cb1;
caexzbcl(1:iebc,1:jb,kbbc)=ca1;
cbexzbcl(1:iebc,1:jb,kbbc)=cb1;
caeyzbcl(1:iebc,1:je,kbbc)=ca1;
cbeyzbcl(1:iebc,1:je,kbbc)=cb1;
caexzbcr(1:iebc,1:jb,kbbc)=ca1;
cbexzbcr(1:iebc,1:jb,kbbc)=cb1;
caeyzbcr(1:iebc,1:je,kbbc)=ca1;
cbeyzbcr(1:iebc,1:je,kbbc)=cb1;

for k=1:kebc                               % 与sigmazs有关的量
  y1=(kebc-k+1)*ds;
  y2=(kebc-k)*ds;
  sigmaz=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
  sigmazs=sigmaz*(muz/(epsz*eps(1)));
  da1=exp(-sigmazs*dt/muz);
  db1=(1.0-da1)/(sigmazs*ds);  
  dahxzbcd(1:ib,1:je,k)=da1;
  dbhxzbcd(1:ib,1:je,k)=db1;
  dahyzbcd(1:ie,1:jb,k)=da1;
  dbhyzbcd(1:ie,1:jb,k)=db1;   
  dahxzbcf(2:iefbc,1:jebc,k)=da1;
  dbhxzbcf(2:iefbc,1:jebc,k)=db1;   
  dahyzbcf(1:iefbc,1:jebc,k)=da1;
  dbhyzbcf(1:iefbc,1:jebc,k)=db1;   
  dahxzbcb(2:iefbc,1:jebc,k)=da1;
  dbhxzbcb(2:iefbc,1:jebc,k)=db1;   
  dahyzbcb(1:iefbc,1:jebc,k)=da1;
  dbhyzbcb(1:iefbc,1:jebc,k)=db1;   
  dahxzbcl(2:iebc,1:je,k)=da1;
  dbhxzbcl(2:iebc,1:je,k)=db1;   
  dahyzbcl(1:iebc,1:jb,k)=da1;
  dbhyzbcl(1:iebc,1:jb,k)=db1;   
  dahxzbcr(1:iebc,1:je,k)=da1;
  dbhxzbcr(1:iebc,1:je,k)=db1;   
  dahyzbcr(1:iebc,1:jb,k)=da1;
  dbhyzbcr(1:iebc,1:jb,k)=db1;   

  caezxbcd(2:ie,1:jb,k)=ca(1);              %与sigmaz、sigmazs无关的量      
  cbezxbcd(2:ie,1:jb,k)=cb(1);  
  caezybcd(1:ib,2:je,k)=ca(1);                   
  cbezybcd(1:ib,2:je,k)=cb(1);   
  dahxybcd(1:ib,1:je,k)=da(1);  
  dbhxybcd(1:ib,1:je,k)=db(1);  
  dahyxbcd(1:ie,1:jb,k)=da(1);
  dbhyxbcd(1:ie,1:jb,k)=db(1); 
end


%     TOP region

caexybct(1:ie,1:jb,kbbc)=1.0;
cbexybct(1:ie,1:jb,kbbc)=0.0;
caexzbct(1:ie,1:jb,kbbc)=1.0;
cbexzbct(1:ie,1:jb,kbbc)=0.0;
caeyxbct(1:ib,1:je,kbbc)=1.0;
cbeyxbct(1:ib,1:je,kbbc)=0.0;
caeyzbct(1:ib,1:je,kbbc)=1.0;
cbeyzbct(1:ib,1:je,kbbc)=0.0;

dahzxbct(1:ie,1:je,kbbc)=1.0;
dbhzxbct(1:ie,1:je,kbbc)=0.0;
dahzybct(1:ie,1:je,kbbc)=1.0;
dbhzybct(1:ie,1:je,kbbc)=0.0;

for k=2:kebc                               % 与sigmaz有关的量
  z1=(k-0.5)*ds;
  z2=(k-1.5)*ds;  
  sigmaz=bcfactor*(z1^(orderbc+1)-z2^(orderbc+1));
  ca1=exp(-sigmaz*dt/(epsz*eps(1)));
  cb1=(1.0-ca1)/(sigmaz*ds);  
  caexzbct(1:ie,1:jb,k)=ca1;
  cbexzbct(1:ie,1:jb,k)=cb1;  
  caeyzbct(1:ib,1:je,k)=ca1;
  cbeyzbct(1:ib,1:je,k)=cb1;  
  
  caexzbcf(1:iefbc,2:jebc,kebc+ke+k)=ca1;
  cbexzbcf(1:iefbc,2:jebc,kebc+ke+k)=cb1;  
  caeyzbcf(2:iefbc,1:jebc,kebc+ke+k)=ca1;
  cbeyzbcf(2:iefbc,1:jebc,kebc+ke+k)=cb1;  
  caexzbcb(1:iefbc,1:jebc,kebc+ke+k)=ca1;
  cbexzbcb(1:iefbc,1:jebc,kebc+ke+k)=cb1; 
  caeyzbcb(2:iefbc,1:jebc,kebc+ke+k)=ca1;
  cbeyzbcb(2:iefbc,1:jebc,kebc+ke+k)=cb1;  
  caexzbcl(1:iebc,1:jb,kebc+ke+k)=ca1;
  cbexzbcl(1:iebc,1:jb,kebc+ke+k)=cb1;  
  caeyzbcl(2:iebc,1:je,kebc+ke+k)=ca1;
  cbeyzbcl(2:iebc,1:je,kebc+ke+k)=cb1;  
  caexzbcr(1:iebc,1:jb,kebc+ke+k)=ca1;
  cbexzbcr(1:iebc,1:jb,kebc+ke+k)=cb1;  
  caeyzbcr(2:iebc,1:je,kebc+ke+k)=ca1;
  cbeyzbcr(2:iebc,1:je,kebc+ke+k)=cb1;   
end


sigmaz = bcfactor*(0.5*ds)^(orderbc+1); 
ca1=exp(-sigmaz*dt/(epsz*eps(1)));
cb1=(1.0-ca1)/(sigmaz*ds);
caex(1:ie,1:jb,kb)=ca1;
cbex(1:ie,1:jb,kb)=cb1;
caey(1:ib,1:je,kb)=ca1;
cbey(1:ib,1:je,kb)=cb1;

caexzbcf(1:iefbc,1:jebc,kbbc+ke)=ca1;
cbexzbcf(1:iefbc,1:jebc,kbbc+ke)=cb1;
caeyzbcf(2:iefbc,1:jebc,kbbc+ke)=ca1;
cbeyzbcf(2:iefbc,1:jebc,kbbc+ke)=cb1;
caexzbcb(1:iefbc,1:jebc,kbbc+ke)=ca1;
cbexzbcb(1:iefbc,1:jebc,kbbc+ke)=cb1;
caeyzbcb(2:iefbc,1:jebc,kbbc+ke)=ca1;
cbeyzbcb(2:iefbc,1:jebc,kbbc+ke)=cb1;
caexzbcl(1:iebc,1:jb,kbbc+ke)=ca1;
cbexzbcl(1:iebc,1:jb,kbbc+ke)=cb1;
caeyzbcl(2:iebc,1:je,kbbc+ke)=ca1;
cbeyzbcl(2:iebc,1:je,kbbc+ke)=cb1;
caexzbcr(1:iebc,1:jb,kbbc+ke)=ca1;
cbexzbcr(1:iebc,1:jb,kbbc+ke)=cb1;
caeyzbcr(1:iebc,1:je,kbbc+ke)=ca1;
cbeyzbcr(1:iebc,1:je,kbbc+ke)=cb1;



for k=1:kebc                               % 与sigmazs有关的量
  z1=k*ds;
  z2=(k-1)*ds;
  sigmaz=bcfactor*(z1^(orderbc+1)-z2^(orderbc+1));
  sigmazs=sigmaz*(muz/(epsz*eps(1)));
  da1=exp(-sigmazs*dt/muz);
  db1=(1.0-da1)/(sigmazs*ds);  
  dahxzbct(1:ib,1:je,k)=da1;
  dbhxzbct(1:ib,1:je,k)=db1;
  dahyzbct(1:ie,1:jb,k)=da1;
  dbhyzbct(1:ie,1:jb,k)=db1;   
  dahxzbcf(2:iefbc,1:jebc,kebc+ke+k)=da1;
  dbhxzbcf(2:iefbc,1:jebc,kebc+ke+k)=db1;   
  dahyzbcf(1:iefbc,1:jebc,kebc+ke+k)=da1;
  dbhyzbcf(1:iefbc,1:jebc,kebc+ke+k)=db1;   
  dahxzbcb(2:iefbc,1:jebc,kebc+ke+k)=da1;
  dbhxzbcb(2:iefbc,1:jebc,kebc+ke+k)=db1;   
  dahyzbcb(1:iefbc,1:jebc,kebc+ke+k)=da1;
  dbhyzbcb(1:iefbc,1:jebc,kebc+ke+k)=db1;   
  dahxzbcl(2:iebc,1:je,kebc+ke+k)=da1;
  dbhxzbcl(2:iebc,1:je,kebc+ke+k)=db1;   
  dahyzbcl(1:iebc,1:jb,kebc+ke+k)=da1;
  dbhyzbcl(1:iebc,1:jb,kebc+ke+k)=db1;   
  dahxzbcr(1:iebc,1:je,kebc+ke+k)=da1;
  dbhxzbcr(1:iebc,1:je,kebc+ke+k)=db1;   
  dahyzbcr(1:iebc,1:jb,kebc+ke+k)=da1;
  dbhyzbcr(1:iebc,1:jb,kebc+ke+k)=db1;   

    %与sigmaz、sigmazs无关的量
  caexybct(1:ie,2:je,k)=ca(1);                   
  cbexybct(1:ie,2:je,k)=cb(1);  
  caeyxbct(2:ie,1:je,k)=ca(1);                   
  cbeyxbct(2:ie,1:je,k)=cb(1);   
  caezxbct(2:ie,1:jb,k)=ca(1);                   
  cbezxbct(2:ie,1:jb,k)=cb(1);  
  caezybct(1:ib,2:je,k)=ca(1);                   
  cbezybct(1:ib,2:je,k)=cb(1);      
  dahxybct(1:ib,1:je,k)=da(1);  
  dbhxybct(1:ib,1:je,k)=db(1);  
  dahyxbct(1:ie,1:jb,k)=da(1);
  dbhyxbct(1:ie,1:jb,k)=db(1);   
  dahzxbct(1:ie,1:je,k)=da(1);
  dbhzxbct(1:ie,1:je,k)=db(1);  
  dahzybct(1:ie,1:je,k)=da(1);
  dbhzybct(1:ie,1:je,k)=db(1); 
end


% ***********************************************************************
%     Movie initialization
% ***********************************************************************

tview(:,:)=real(hz(:,:,ks));
sview(:,:)=real(hz(is,:,:));

subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview');
shading flat;
caxis([-1.0 1.0]);
colorbar;
axis image;
title(['real Hz(i,j,k=',ks,') ', 'time step = 0']);
xlabel('i coordinate');
ylabel('j coordinate');

subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview');
shading flat;
caxis([-1.0 1.0]);
colorbar;
axis image;