threedimensional_fdtdcpml.m

FDTD three-dimensional CPML

         den_ey(j) = 1.0/dy;
      end
   end
   kk = nzPML_2;
   for k = 1:Kmax-1
      if (k <= nzPML_1) 
         den_ez(k) = 1.0/(kappae_z_PML_1(k)*dz);
      elseif (k >= Kmax-nzPML_2)
         den_ez(k) = 1.0/(kappae_z_PML_2(kk)*dz);
         kk = kk - 1;
      else
         den_ez(k) = 1.0/dz;
      end
   end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  BEGIN TIME STEP
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
for n = 1:nmax;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  UPDATE Ex
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   for k = 1:Kmax-1
      for i = 1:Imax-1
    for j = 2:Jmax-1
              if (i >= istart-1 && i <= iend && j >= jstart && j <= jend && k >= kstart && k <= kend)
                 Ex(i,j,k) = 0.0;
              else
          Ex(i,j,k) = CA(i,j,k) * Ex(i,j,k) + CB(i,j,k) * ( (Hz(i,j,k) - Hz(i,j-1,k))*den_ey(j)  +  (Hy(i,j,k) - Hy(i,j,k+1))*den_ez(k) );
              end
    end
      end
      for i = 1:Imax-1
%.....................................................................
%  PML for bottom Ex, j-direction
%.....................................................................
         for j = 2:nyPML_1
         psi_Exy_1(i,j,k) = be_y_1(j)*psi_Exy_1(i,j,k)+ ce_y_1(j) *(Hz(i,j,k) - Hz(i,j-1,k))/dy;
         Ex(i,j,k) = Ex(i,j,k) + CB(i,j,k)*psi_Exy_1(i,j,k);
         end
%.....................................................................
%  PML for top Ex, j-direction
%.....................................................................
         jj = nyPML_2;
         for j = Jmax+1-nyPML_2:Jmax-1
         psi_Exy_2(i,jj,k) = be_y_2(jj)*psi_Exy_2(i,jj,k)+ ce_y_2(jj) *(Hz(i,j,k) - Hz(i,(j-1),k))/dy;
         Ex(i,j,k) = Ex(i,j,k) + CB(i,j,k)*psi_Exy_2(i,jj,k);
            jj = jj-1;
         end
      end
   end
   for i = 1:Imax-1
      for j = 2:Jmax-1
%.....................................................................
%  PML for bottom Ex, k-direction
%.....................................................................
         for k = 1:nzPML_1
         psi_Exz_1(i,j,k) = be_z_1(k)*psi_Exz_1(i,j,k)+ ce_z_1(k) *(Hy(i,j,k) - Hy(i,j,k+1))/dz;
         Ex(i,j,k) = Ex(i,j,k) + CB(i,j,k)*psi_Exz_1(i,j,k);
         end
%.....................................................................
%  PML for top Ex, k-direction
%.....................................................................
         kk = nzPML_2;
         for k = Kmax-nzPML_2:Kmax-1
         psi_Exz_2(i,j,kk) = be_z_2(kk)*psi_Exz_2(i,j,kk)+ ce_z_2(kk) *(Hy(i,j,k) - Hy(i,j,k+1))/dz;
         Ex(i,j,k) = Ex(i,j,k) + CB(i,j,k)*psi_Exz_2(i,j,kk);
            kk = kk-1;
         end
      end
   end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  UPDATE Ey
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   for k = 1:Kmax-1
      for i = 2:Imax-1
    for j = 1:Jmax-1
              if (i >= istart && i <= iend && j >= jstart-1 && j <= jend && k >= kstart && k <= kend)
                 Ey(i,j,k) = 0.0;
              else
                 Ey(i,j,k) = CA(i,j,k) * Ey(i,j,k) + CB(i,j,k) *( (Hz(i-1,j,k) - Hz(i,j,k))*den_ex(i) +  (Hx(i,j,k+1) - Hx(i,j,k))*den_ez(k) );
              end
    end
      end
      for j = 1:Jmax-1
%.....................................................................
%  PML for bottom Ey, i-direction
%.....................................................................
         for i = 2:nxPML_1
       psi_Eyx_1(i,j,k) = be_x_1(i)*psi_Eyx_1(i,j,k)+ ce_x_1(i)*(Hz(i-1,j,k) - Hz(i,j,k))/dx;
       Ey(i,j,k) = Ey(i,j,k) + CB(i,j,k)*psi_Eyx_1(i,j,k);
         end
%.....................................................................
%  PML for top Ey, i-direction
%.....................................................................
         ii = nxPML_2;
         for i = Imax+1-nxPML_2:Imax-1
       psi_Eyx_2(ii,j,k) = be_x_2(ii)*psi_Eyx_2(ii,j,k)+ ce_x_2(ii)*(Hz(i-1,j,k) - Hz(i,j,k))/dx;
       Ey(i,j,k) = Ey(i,j,k) + CB(i,j,k)*psi_Eyx_2(ii,j,k);
            ii = ii-1;
         end
      end
   end
   for i = 2:Imax-1
      for j = 1:Jmax-1
%.....................................................................
%  PML for bottom Ey, k-direction
%.....................................................................
         for k = 1:nzPML_1
       psi_Eyz_1(i,j,k) = be_z_1(k)*psi_Eyz_1(i,j,k)+ ce_z_1(k)*(Hx(i,j,k+1) - Hx(i,j,k))/dz;
       Ey(i,j,k) = Ey(i,j,k) + CB(i,j,k)*psi_Eyz_1(i,j,k);
         end
%.....................................................................
%  PML for top Ey, k-direction
%.....................................................................
         kk = nzPML_2;
         for k = Kmax-nzPML_2:Kmax-1;
     psi_Eyz_2(i,j,kk) = be_z_2(kk)*psi_Eyz_2(i,j,kk)+ ce_z_2(kk)*(Hx(i,j,k+1) - Hx(i,j,k))/dz;
     Ey(i,j,k) = Ey(i,j,k) + CB(i,j,k)*psi_Eyz_2(i,j,kk);
            kk = kk-1;
         end
      end
   end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  UPDATE Ez
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   for k = 2:Kmax-1
      for i = 2:Imax-1
         for j = 2:Jmax-1
            Ez(i,j,k) = CA(i,j,k) * Ez(i,j,k) + CB(i,j,k)* ((Hy(i,j,k) - Hy(i-1,j,k))*den_ex(i) + (Hx(i,j-1,k) - Hx(i,j,k))*den_ey(j));
         end
      end
      for j = 2:Jmax-1
%.....................................................................
%  PML for bottom Ez, x-direction
%.....................................................................
         for i = 2:nxPML_1
          psi_Ezx_1(i,j,k) = be_x_1(i)*psi_Ezx_1(i,j,k)+ ce_x_1(i) *(Hy(i,j,k) - Hy(i-1,j,k))/dx;
          Ez(i,j,k) = Ez(i,j,k) + CB(i,j,k)*psi_Ezx_1(i,j,k);
         end
%.....................................................................
%  PML for top Ez, x-direction
%.....................................................................
         ii = nxPML_2;
         for i = Imax+1-nxPML_2:Imax-1
          psi_Ezx_2(ii,j,k) = be_x_2(ii)*psi_Ezx_2(ii,j,k)+ ce_x_2(ii) *(Hy(i,j,k) - Hy(i-1,j,k))/dx;
          Ez(i,j,k) = Ez(i,j,k) + CB(i,j,k)*psi_Ezx_2(ii,j,k);
            ii = ii-1;
         end
      end
      for i = 2:Imax-1
%.....................................................................
%  PML for bottom Ez, y-direction
%.....................................................................
         for j = 2:nyPML_1
            psi_Ezy_1(i,j,k) = be_y_1(j)*psi_Ezy_1(i,j,k)+ ce_y_1(j)*(Hx(i,j-1,k) - Hx(i,j,k))/dy;
            Ez(i,j,k) = Ez(i,j,k) + CB(i,j,k)*psi_Ezy_1(i,j,k);
         end
%.....................................................................
%  PML for top Ez, y-direction
%.....................................................................
         jj = nyPML_2;
         for j = Jmax+1-nyPML_2:Jmax-1
            psi_Ezy_2(i,jj,k) = be_y_2(jj)*psi_Ezy_2(i,jj,k)+ ce_y_2(jj)*(Hx(i,j-1,k) - Hx(i,j,k))/dy;
            Ez(i,j,k) = Ez(i,j,k) + CB(i,j,k)*psi_Ezy_2(i,jj,k);
            jj = jj-1;
         end
      end
   end
%-----------------------------------------------------------------------
%   SOURCE
%-----------------------------------------------------------------------
Ez(isource01,jsource01,1:Kmax)=Ez(isource01,jsource01,1:Kmax)+srcconst*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));
%Hz(isource02,jsource02,2:Kmax-1)=Hz(isource02,jsource02,2:Kmax-1)+srcconst*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  UPDATE Hx
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   for k = 2:Kmax-1
      for i = 1:Imax-1
    for j = 1:Jmax-1
       Hx(i,j,k) = DA * Hx(i,j,k) + DB * ( (Ez(i,j,k) - Ez(i,j+1,k))*den_hy(j)  +  (Ey(i,j,k) - Ey(i,j,k-1))*den_hz(k) );
    end
      end
      
      for i = 1:Imax-1
%.....................................................................
%  PML for bottom Hx, j-direction
%.....................................................................
         for j = 1:nyPML_1-1
         psi_Hxy_1(i,j,k) = bh_y_1(j)*psi_Hxy_1(i,j,k)+ ch_y_1(j) *(Ez(i,j,k) - Ez(i,j+1,k))/dy;
         Hx(i,j,k) = Hx(i,j,k) + DB*psi_Hxy_1(i,j,k);
         end
%.....................................................................
%  PML for top Hx, j-direction
%.....................................................................
         jj = nyPML_2-1;
         for j = Jmax+1-nyPML_2:Jmax-1
         psi_Hxy_2(i,jj,k) = bh_y_2(jj)*psi_Hxy_2(i,jj,k) + ch_y_2(jj) *(Ez(i,j,k) - Ez(i,(j+1),k))/dy;
         Hx(i,j,k) = Hx(i,j,k) + DB*psi_Hxy_2(i,jj,k);
            jj = jj-1;
         end
      end
   end
   for i = 1:Imax-1
      for j = 1:Jmax-1
%.....................................................................
%  PML for bottom Hx, k-direction
%.....................................................................
         for k = 2:nzPML_1
         psi_Hxz_1(i,j,k-1) = bh_z_1(k-1)*psi_Hxz_1(i,j,k-1)+ ch_z_1(k-1) *(Ey(i,j,k) - Ey(i,j,k-1))/dz;
         Hx(i,j,k) = Hx(i,j,k) + DB*psi_Hxz_1(i,j,k-1);
         end
%.....................................................................
%  PML for top Hx, k-direction
%.....................................................................
         kk = nzPML_2-1;
         for k = Kmax+1-nzPML_2:Kmax-1
         psi_Hxz_2(i,j,kk) = bh_z_2(kk)*psi_Hxz_2(i,j,kk)+ ch_z_2(kk) *(Ey(i,j,k) - Ey(i,j,k-1))/dz;
         Hx(i,j,k) = Hx(i,j,k) + DB*psi_Hxz_2(i,j,kk);
           kk = kk-1;
         end
      end
   end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  UPDATE Hy
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   for k = 2:Kmax-1
      for i = 1:Imax-1
    for j = 1:Jmax-1
            Hy(i,j,k) = DA * Hy(i,j,k) + DB *( (Ez(i+1,j,k) - Ez(i,j,k))*den_hx(i) + (Ex(i,j,k-1) - Ex(i,j,k))*den_hz(k) );
    end
      end
      for j = 1:Jmax-1
%.....................................................................
%  PML for bottom Hy, i-direction
%.....................................................................
         for i = 1:nxPML_1-1
       psi_Hyx_1(i,j,k) = bh_x_1(i)*psi_Hyx_1(i,j,k)+ ch_x_1(i)*(Ez(i+1,j,k) - Ez(i,j,k))/dx;
       Hy(i,j,k) = Hy(i,j,k) + DB*psi_Hyx_1(i,j,k);
         end
%.....................................................................
%  PML for top Hy, i-direction
%.....................................................................
         ii = nxPML_2-1;
         for i = Imax+1-nxPML_2:Imax-1
       psi_Hyx_2(ii,j,k) = bh_x_2(ii)*psi_Hyx_2(ii,j,k)+ ch_x_2(ii)*(Ez(i+1,j,k) - Ez(i,j,k))/dx;
       Hy(i,j,k) = Hy(i,j,k) + DB*psi_Hyx_2(ii,j,k);
            ii = ii-1;
         end
      end
   end
   for i = 1:Imax-1
      for j = 1:Jmax-1
%.....................................................................
%  PML for bottom Hy, k-direction
%.....................................................................
         for k = 2:nzPML_1
       psi_Hyz_1(i,j,k-1) = bh_z_1(k-1)*psi_Hyz_1(i,j,k-1)+ ch_z_1(k-1)*(Ex(i,j,k-1) - Ex(i,j,k))/dz;
       Hy(i,j,k) = Hy(i,j,k) + DB*psi_Hyz_1(i,j,k-1);
         end
%.....................................................................
%  PML for top Hy, k-direction
%.....................................................................
         kk = nzPML_2-1;
         for k = Kmax+1-nzPML_2:Kmax-1
     psi_Hyz_2(i,j,kk) = bh_z_2(kk)*psi_Hyz_2(i,j,kk)+ ch_z_2(kk)*(Ex(i,j,k-1) - Ex(i,j,k))/dz;
     Hy(i,j,k) = Hy(i,j,k) + DB*psi_Hyz_2(i,j,kk);
            kk = kk-1;
         end
      end
   end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  UPDATE Hz
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   for k = 1:Kmax-1
      for i = 1:Imax-1
        for j = 1:Jmax-1
            Hz(i,j,k) = DA * Hz(i,j,k) + DB* ((Ey(i,j,k) - Ey(i,j,k))*den_hx(i) +  (Ex(i,j+1,k) - Ex(i,j,k))*den_hy(j));
        end
      end
      for j = 1:Jmax-1
%.....................................................................
%  PML for bottom Hz, x-direction
%.....................................................................
         for i = 1:nxPML_1-1
          psi_Hzx_1(i,j,k) = bh_x_1(i)*psi_Hzx_1(i,j,k)+ ch_x_1(i) *(Ey(i,j,k) - Ey(i+1,j,k))/dx;
          Hz(i,j,k) = Hz(i,j,k) + DB*psi_Hzx_1(i,j,k);
         end
%.....................................................................
%  PML for top Hz, x-direction
%.....................................................................
         ii = nxPML_2-1;
         for i = Imax+1-nxPML_2:Imax-1
          psi_Hzx_2(ii,j,k) = bh_x_2(ii)*psi_Hzx_2(ii,j,k)+ ch_x_2(ii) *(Ey(i,j,k) - Ey(i,j,k))/dx;
          Hz(i,j,k) = Hz(i,j,k) + DB*psi_Hzx_2(ii,j,k);
            ii = ii-1;
         end
      end
      for i = 1:Imax-1
%.....................................................................
%  PML for bottom Hz, y-direction
%.....................................................................
         for j = 1:nyPML_1-1
            psi_Hzy_1(i,j,k) = bh_y_1(j)*psi_Hzy_1(i,j,k)+ ch_y_1(j)*(Ex(i,j+1,k) - Ex(i,j,k))/dy;
            Hz(i,j,k) = Hz(i,j,k) + DB*psi_Hzy_1(i,j,k);
         end
%.....................................................................
%  PML for top Hz, y-direction
%.....................................................................
         jj = nyPML_2-1;
         for j = Jmax+1-nyPML_2:Jmax-1
            psi_Hzy_2(i,jj,k) = bh_y_2(jj)*psi_Hzy_2(i,jj,k)+ ch_y_2(jj)*(Ex(i,j+1,k) - Ex(i,j,k))/dy;
            Hz(i,j,k) = Hz(i,j,k) + DB*psi_Hzy_2(i,jj,k);
            jj = jj-1;
         end
      end
   end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%   VISUALIZATION
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if mod(n,Nplt_jmp)==0;
NoFig=NoFig+1;
timestep=int2str(n);

tviewez(:,:)=Ez(:,:,ksource01);
tviewhy(:,:)=Hy(:,:,ksource01);
tviewhx(:,:)=Hx(:,:,ksource01);
tviewhz(:,:)=Hz(:,:,ksource02);
sview(:,:)=Ez(:,jsource01,:);

subplot('position',[0 0.55 0.4 0.4]),pcolor(tviewez');
shading flat;
colorbar;
axis('square')
title(['TM-Ez, time step =',timestep]);
xlabel('i coordinate');
ylabel('j coordinate');

subplot('position',[0 0.05 0.4 0.4]),pcolor(tviewhy');
shading flat;
colorbar;
axis('square')
title(['TM-Hy']);
xlabel('i coordinate');
ylabel('j coordinate');

subplot('position',[0.5 0.55 0.4 0.4]),pcolor(tviewhx');
shading flat;
colorbar;
axis('square')
title(['TM-Hx']);
xlabel('i coordinate');
ylabel('j coordinate');
%------------------------------------------------------
%
subplot('position',[0.55 0.05 0.4 0.4]),pcolor(sview');
shading flat;
colorbar;
rect(1:2)=[0 0];
title(['Ez - k coordinate']);
xlabel('i coordinate');
ylabel('k coordinate');
%}

%{
load view_grid.dat
for j=1:Jmax
    image(1:Imax,j)=view_grid(1+(j-1)*Imax:J*Imax);
end
pcolor(image');shading flat;
%}
%--------------------------------------------------------------------------
%{
x=ieT:ihT;
y=jeT:jhT;
timestep=int2str(n);
[X Y]=meshgrid(x,y);
subplot(1,2,1),pcolor(Ez'); 
mesh(Ez)
axis([-inf inf -inf inf -0.05 0.05]);
xlabel('i-axis')
ylabel('j-axis')
zlabel('Ez')
title(['