⭐ 欢迎来到虫虫下载站! | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站

📄 threedimensional_fdtdcpml.m

📁 FDTD three-dimensional CPML
💻 M
📖 第 1 页 / 共 2 页
字号:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FDTD_3D_CPML
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
clear all;
clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nmax =220;
  
epsR = 1.0;
sigM1 = 0.0; 

pi = 3.14159265358979;
c0 = 2.99792458E8;
muO = 4.0 * pi * 1.0E-7; 
epsO = 1.0/(c0*c0*muO);
  
freq=1.9341e+14;
lambda=c0/freq;         
omega=2.0*pi*freq;

dx = 1.0E-3; dy = 1.0E-3; dz = 1.0E-3;
dt = 0.99 / (c0*(1.0/dx^2+1.0/dy^2+1.0/dz^2)^0.5);

tw = 53.0E-12;      tO = 4.0*tw;

Imax = 50; 
Jmax = 50;
Kmax = 50;
%
istart = (Imax-1)/2-11; 
iend = istart+24; 

jstart = istart;
jend = jstart + 99; 

kstart = Kmax/2; 
kend = kstart;

isource = istart; 
jsource = jstart; 
ksource = kend+1;

irecv1 = iend; 
jrecv1 = jend+1; 
krecv1 = kend;
%}
%
isource01 = Imax/2; 
jsource01 = Jmax/2; 
ksource01 = Kmax/2;
isource02 = (Imax)/2; 
jsource02 = (Jmax)/2; 
ksource02 = (Kmax)/2;
%}
% PML thickness in each direction
nxPML_1 = 20; nxPML_2 = nxPML_1; nyPML_1 = nxPML_1;
nyPML_2 = nxPML_1; nzPML_1 = nxPML_1; nzPML_2 = nxPML_2;
%-------------------------
m = 3; ma = 1;

sig_x_max = 0.75 * (0.8*(m+1)/(dx*(muO/epsO*epsR)^0.5));
sig_y_max = 0.75 * (0.8*(m+1)/(dy*(muO/epsO*epsR)^0.5));
sig_z_max = 0.75 * (0.8*(m+1)/(dz*(muO/epsO*epsR)^0.5));
alpha_x_max = 0.24;
alpha_y_max = alpha_x_max; 
alpha_z_max = alpha_x_max;
kappa_x_max = 15.0;
kappa_y_max = kappa_x_max; 
kappa_z_max = kappa_x_max;

Ex =zeros(Imax-1, Jmax, Kmax-1);Ey =zeros(Imax-1, Jmax, Kmax-1);Ez =zeros(Imax, Jmax, Kmax);
Hx =zeros(Imax,Jmax-1, Kmax);Hy =zeros(Imax-1, Jmax, Kmax);Hz =zeros(Imax-1, Jmax-1, Kmax-1); 
      
CA  =zeros(Imax, Jmax, Kmax);
CB  =zeros(Imax, Jmax, Kmax);
sig =zeros(Imax, Jmax, Kmax);
eps=zeros(Imax, Jmax, Kmax); 

NoFig=29;
Nplt_jmp=5;
%***********************************************************************
%     Differentiated Gaussian pulse excitation
%***********************************************************************
rtau=50.0e-12;
tau=rtau/dt;
ndelay=3*tau;
srcconst=-dt*3.0e+11;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  CPML 
psi_Exy_1=zeros(Imax-1,nyPML_1,Kmax-1);psi_Exy_2=zeros(Imax-1,nyPML_2,Kmax-1);

psi_Exz_1=zeros(Imax-1,Jmax,nzPML_1);psi_Exz_2=zeros(Imax-1,Jmax,nzPML_2);

psi_Eyz_1=zeros(Imax,Jmax-1,nzPML_1);psi_Eyz_2=zeros(Imax,Jmax-1,nzPML_2);

psi_Eyx_1=zeros(nxPML_1,Jmax-1,Kmax-1);psi_Eyx_2=zeros(nxPML_2,Jmax-1,Kmax-1);

psi_Ezx_1=zeros(nxPML_1,Jmax,Kmax);psi_Ezx_2=zeros(nxPML_2,Jmax,Kmax);

psi_Ezy_1=zeros(Imax,nyPML_1,Kmax);psi_Ezy_2=zeros(Imax,nyPML_1,Kmax);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
psi_Hxy_1=zeros(Imax,nyPML_1-1,Kmax);psi_Hxy_2=zeros(Imax,nyPML_2-1,Kmax);

psi_Hxz_1=zeros(Imax,Jmax-1,nzPML_1-1);psi_Hxz_2=zeros(Imax,Jmax-1,nzPML_2-1);

psi_Hyz_1=zeros(Imax-1,Jmax,nzPML_1-1);psi_Hyz_2=zeros(Imax-1,Jmax,nzPML_2-1);

psi_Hyx_1=zeros(nxPML_1-1,Jmax,Kmax);psi_Hyx_2=zeros(nxPML_2-1,Jmax,Kmax); 

psi_Hzx_1=zeros(nxPML_1-1,Jmax-1,Kmax-1);psi_Hzx_2=zeros(nxPML_2-1,Jmax-1,Kmax-1);

psi_Hzy_1=zeros(Imax-1,nyPML_1-1,Kmax-1);psi_Hzy_2=zeros(Imax-1,nyPML_2-1,Kmax-1);        
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
% E-Field
% x-dir
be_x_1=zeros(nxPML_1,1);
ce_x_1=zeros(nxPML_1,1);
alphae_x_PML_1=zeros(nxPML_1,1);
sige_x_PML_1=zeros(nxPML_1,1);
kappae_x_PML_1=zeros(nxPML_1,1);
be_x_2=zeros(nxPML_2,1);
ce_x_2=zeros(nxPML_2,1);
alphae_x_PML_2=zeros(nxPML_2,1);
sige_x_PML_2=zeros(nxPML_2,1);
kappae_x_PML_2=zeros(nxPML_2,1);
% y-dir
be_y_1=zeros(nyPML_1,1);
ce_y_1=zeros(nyPML_1,1);
alphae_y_PML_1=zeros(nyPML_1,1);
sige_y_PML_1=zeros(nyPML_1,1);
kappae_y_PML_1=zeros(nyPML_1,1);
be_y_2=zeros(nyPML_2,1);
ce_y_2=zeros(nyPML_2,1);
alphae_y_PML_2=zeros(nyPML_2,1);
sige_y_PML_2=zeros(nyPML_2,1);
kappae_y_PML_2=zeros(nyPML_2,1);
% z-dir
be_z_1=zeros(nzPML_1,1);
ce_z_1=zeros(nzPML_1,1);
alphae_z_PML_1=zeros(nzPML_1,1);
sige_z_PML_1=zeros(nzPML_1,1);
kappae_z_PML_1=zeros(nzPML_1,1);
be_z_2=zeros(nzPML_2,1);
ce_z_2=zeros(nzPML_2,1);
alphae_z_PML_2=zeros(nzPML_2,1);
sige_z_PML_2=zeros(nzPML_2,1);
kappae_z_PML_2=zeros(nzPML_2,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% H-Field
% x-dir
bh_x_1=zeros(nxPML_1-1,1);
ch_x_1=zeros(nxPML_1-1,1);
alphah_x_PML_1=zeros(nxPML_1-1,1);
sigh_x_PML_1=zeros(nxPML_1-1,1);
kappah_x_PML_1=zeros(nxPML_1-1,1);
bh_x_2=zeros(nxPML_2-1,1);
ch_x_2=zeros(nxPML_2-1,1);
alphah_x_PML_2=zeros(nxPML_2-1,1);
sigh_x_PML_2=zeros(nxPML_2-1,1);
kappah_x_PML_2=zeros(nxPML_2-1,1);
% y-dir
bh_y_1=zeros(nyPML_1-1,1);
ch_y_1=zeros(nyPML_1-1,1);
alphah_y_PML_1=zeros(nyPML_1-1,1);
sigh_y_PML_1=zeros(nyPML_1-1,1);
kappah_y_PML_1=zeros(nyPML_1-1,1);
bh_y_2=zeros(nyPML_2-1,1);
ch_y_2=zeros(nyPML_2-1,1);
alphah_y_PML_2=zeros(nyPML_2-1,1);
sigh_y_PML_2=zeros(nyPML_2-1,1);
kappah_y_PML_2=zeros(nyPML_2-1,1);
% z-dir
bh_z_1=zeros(nzPML_1-1,1);
ch_z_1=zeros(nzPML_1-1,1);
alphah_z_PML_1=zeros(nzPML_1-1,1);
sigh_z_PML_1=zeros(nzPML_1-1,1);
kappah_z_PML_1=zeros(nzPML_1-1,1);
bh_z_2=zeros(nzPML_2-1,1);
ch_z_2=zeros(nzPML_2-1,1);
alphah_z_PML_2=zeros(nzPML_2-1,1);
sigh_z_PML_2=zeros(nzPML_2-1,1);
kappah_z_PML_2=zeros(nzPML_2-1,1);

% Denominators for the update equations 
den_ex=zeros(Imax-1,1);
den_hx=zeros(Imax-1,1);
 
den_ey=zeros(Jmax-1,1);
den_hy=zeros(Jmax-1,1);
 
den_ez=zeros(Kmax-1,1);
den_hz=zeros(Kmax-1,1);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  INITIALIZE VARIABLES
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
sig(:,:,:) = sigM1;
eps(:,:,:) = epsR*epsO;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  SET CPML PARAMETERS IN EACH DIRECTION
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  x-dir
%
   for i = 1:nxPML_1
      sige_x_PML_1(i) = sig_x_max * ( (nxPML_1 - i) / (nxPML_1 - 1.0) )^m;
      alphae_x_PML_1(i) = alpha_x_max*((i-1.0)/(nxPML_1-1.0))^ma;
      kappae_x_PML_1(i) = 1.0+(kappa_x_max-1.0)*((nxPML_1 - i) / (nxPML_1 - 1.0))^m;
      be_x_1(i) = exp(-(sige_x_PML_1(i) / kappae_x_PML_1(i) + alphae_x_PML_1(i))*dt/epsO);
      if ((sige_x_PML_1(i) == 0.0) && (alphae_x_PML_1(i) == 0.0) && (i == nxPML_1)) 
         ce_x_1(i) = 0.0;
      else
         ce_x_1(i) = sige_x_PML_1(i)*(be_x_1(i)-1.0)/(sige_x_PML_1(i)+kappae_x_PML_1(i)*alphae_x_PML_1(i))/ kappae_x_PML_1(i);
      end
   end
   for i = 1:nxPML_1-1
      sigh_x_PML_1(i) = sig_x_max * ( (nxPML_1 - i - 0.5)/(nxPML_1-1.0))^m;
      alphah_x_PML_1(i) = alpha_x_max*((i-0.5)/(nxPML_1-1.0))^ma;
      kappah_x_PML_1(i) = 1.0+(kappa_x_max-1.0)*((nxPML_1 - i - 0.5) / (nxPML_1 - 1.0))^m;
      bh_x_1(i) = exp(-(sigh_x_PML_1(i) / kappah_x_PML_1(i) + alphah_x_PML_1(i))*dt/epsO);
      ch_x_1(i) = sigh_x_PML_1(i)*(bh_x_1(i)-1.0)/(sigh_x_PML_1(i)+kappah_x_PML_1(i)*alphah_x_PML_1(i))/ kappah_x_PML_1(i);
   end
%-------------------
   for i = 1:nxPML_2
      sige_x_PML_2(i) = sig_x_max * ( (nxPML_2 - i) / (nxPML_2 - 1.0) )^m;
      alphae_x_PML_2(i) = alpha_x_max*((i-1.0)/(nxPML_2-1.0))^ma;
      kappae_x_PML_2(i) = 1.0+(kappa_x_max-1.0)*((nxPML_2 - i) / (nxPML_2 - 1.0))^m;
      be_x_2(i) = exp(-(sige_x_PML_2(i) / kappae_x_PML_2(i) + alphae_x_PML_2(i))*dt/epsO);
      if ((sige_x_PML_2(i) == 0.0) && (alphae_x_PML_2(i) == 0.0) && (i == nxPML_2)) 
         ce_x_2(i) = 0.0;
      else
         ce_x_2(i) = sige_x_PML_2(i)*(be_x_2(i)-1.0)/(sige_x_PML_2(i)+kappae_x_PML_2(i)*alphae_x_PML_2(i))/ kappae_x_PML_2(i);
      end
   end
   for i = 1:nxPML_2-1
      sigh_x_PML_2(i) = sig_x_max * ( (nxPML_2 - i - 0.5)/(nxPML_2-1.0))^m;
      alphah_x_PML_2(i) = alpha_x_max*((i-0.5)/(nxPML_2-1.0))^ma;
      kappah_x_PML_2(i) = 1.0+(kappa_x_max-1.0)*((nxPML_2 - i - 0.5) / (nxPML_2 - 1.0))^m;
      bh_x_2(i) = exp(-(sigh_x_PML_2(i) / kappah_x_PML_2(i) + alphah_x_PML_2(i))*dt/epsO);
      ch_x_2(i) = sigh_x_PML_2(i)*(bh_x_2(i)-1.0)/(sigh_x_PML_2(i)+kappah_x_PML_2(i)*alphah_x_PML_2(i))/ kappah_x_PML_2(i);
   end
%  y-dir
%
   for j = 1:nyPML_1
      sige_y_PML_1(j) = sig_y_max * ( (nyPML_1 - j ) / (nyPML_1 - 1.0) )^m;
      alphae_y_PML_1(j) = alpha_y_max*((j-1)/(nyPML_1-1.0))^ma;
      kappae_y_PML_1(j) = 1.0+(kappa_y_max-1.0)*((nyPML_1 - j) / (nyPML_1 - 1.0))^m;
      be_y_1(j) = exp(-(sige_y_PML_1(j) / kappae_y_PML_1(j) + alphae_y_PML_1(j))*dt/epsO);
      if ((sige_y_PML_1(j) == 0.0) && (alphae_y_PML_1(j) == 0.0) && (j == nyPML_1))
         ce_y_1(j) = 0.0;
      else
         ce_y_1(j) = sige_y_PML_1(j)*(be_y_1(j)-1.0)/(sige_y_PML_1(j)+kappae_y_PML_1(j)*alphae_y_PML_1(j))/ kappae_y_PML_1(j);
      end
   end
   for j = 1:nyPML_1-1
      sigh_y_PML_1(j) = sig_y_max * ( (nyPML_1 - j - 0.5)/(nyPML_1-1.0))^m;
      alphah_y_PML_1(j) = alpha_y_max*((j-0.5)/(nyPML_1-1.0))^ma;
      kappah_y_PML_1(j) = 1.0+(kappa_y_max-1.0)*((nyPML_1 - j - 0.5) / (nyPML_1 - 1.0))^m;
      bh_y_1(j) = exp(-(sigh_y_PML_1(j) / kappah_y_PML_1(j) + alphah_y_PML_1(j))*dt/epsO);
      ch_y_1(j) = sigh_y_PML_1(j)*(bh_y_1(j)-1.0)/(sigh_y_PML_1(j)+kappah_y_PML_1(j)*alphah_y_PML_1(j))/ kappah_y_PML_1(j);
   end
%-------------------
   for j = 1:nyPML_2
      sige_y_PML_2(j) = sig_y_max * ( (nyPML_2 - j ) / (nyPML_2 - 1.0) )^m;
      alphae_y_PML_2(j) = alpha_y_max*((j-1)/(nyPML_2-1.0))^ma;
      kappae_y_PML_2(j) = 1.0+(kappa_y_max-1.0)*((nyPML_2 - j) / (nyPML_2 - 1.0))^m;
      be_y_2(j) = exp(-(sige_y_PML_2(j) / kappae_y_PML_2(j) + alphae_y_PML_2(j))*dt/epsO);
      if ((sige_y_PML_2(j) == 0.0) && (alphae_y_PML_2(j) == 0.0) && (j == nyPML_2))
         ce_y_2(j) = 0.0;
      else
         ce_y_2(j) = sige_y_PML_2(j)*(be_y_2(j)-1.0)/(sige_y_PML_2(j)+kappae_y_PML_2(j)*alphae_y_PML_2(j))/ kappae_y_PML_2(j);
      end
   end
   for j = 1:nyPML_2-1
      sigh_y_PML_2(j) = sig_y_max * ( (nyPML_2 - j - 0.5)/(nyPML_2-1.0))^m;
      alphah_y_PML_2(j) = alpha_y_max*((j-0.5)/(nyPML_2-1.0))^ma;
      kappah_y_PML_2(j) = 1.0+(kappa_y_max-1.0)*((nyPML_2 - j - 0.5) / (nyPML_2 - 1.0))^m;
      bh_y_2(j) = exp(-(sigh_y_PML_2(j) / kappah_y_PML_2(j) + alphah_y_PML_2(j))*dt/epsO);
      ch_y_2(j) = sigh_y_PML_2(j)*(bh_y_2(j)-1.0)/(sigh_y_PML_2(j)+kappah_y_PML_2(j)*alphah_y_PML_2(j))/ kappah_y_PML_2(j);
   end
%  k-dir
% 
   for k = 1:nzPML_1
      sige_z_PML_1(k) = sig_z_max * ( (nzPML_1 - k ) / (nzPML_1 - 1.0) )^m;
      alphae_z_PML_1(k) = alpha_z_max*((k-1)/(nzPML_1-1.0))^ma;
      kappae_z_PML_1(k) = 1.0+(kappa_z_max-1.0)*((nzPML_1 - k) / (nzPML_1 - 1.0))^m;
      be_z_1(k) = exp(-(sige_z_PML_1(k) / kappae_z_PML_1(k) + alphae_z_PML_1(k))*dt/epsO);
      if ((sige_z_PML_1(k) == 0.0) && (alphae_z_PML_1(k) == 0.0) && (k == nzPML_1))
         ce_z_1(k) = 0.0;
      else
         ce_z_1(k) = sige_z_PML_1(k)*(be_z_1(k)-1.0)/(sige_z_PML_1(k)+kappae_z_PML_1(k)*alphae_z_PML_1(k))/ kappae_z_PML_1(k);
      end
   end
   for k = 1:nzPML_1-1
      sigh_z_PML_1(k) = sig_z_max * ( (nzPML_1 - k - 0.5)/(nzPML_1-1.0))^m;
      alphah_z_PML_1(k) = alpha_z_max*((k-0.5)/(nzPML_1-1.0))^ma;
      kappah_z_PML_1(k) = 1.0+(kappa_z_max-1.0)*((nzPML_1 - k - 0.5) / (nzPML_1 - 1.0))^m;
      bh_z_1(k) = exp(-(sigh_z_PML_1(k) / kappah_z_PML_1(k) + alphah_z_PML_1(k))*dt/epsO);
      ch_z_1(k) = sigh_z_PML_1(k)*(bh_z_1(k)-1.0)/(sigh_z_PML_1(k)+kappah_z_PML_1(k)*alphah_z_PML_1(k))/ kappah_z_PML_1(k);
   end   
 %-------------------  
   for k = 1:nzPML_2
      sige_z_PML_2(k) = sig_z_max * ( (nzPML_2 - k ) / (nzPML_2 - 1.0) )^m;
      alphae_z_PML_2(k) = alpha_z_max*((k-1)/(nzPML_2-1.0))^ma;
      kappae_z_PML_2(k) = 1.0+(kappa_z_max-1.0)*((nzPML_2 - k) / (nzPML_2 - 1.0))^m;
      be_z_2(k) = exp(-(sige_z_PML_2(k) / kappae_z_PML_2(k) + alphae_z_PML_2(k))*dt/epsO);
      if ((sige_z_PML_2(k) == 0.0) && (alphae_z_PML_2(k) == 0.0) && (k == nzPML_2))
         ce_z_2(k) = 0.0;
      else
         ce_z_2(k) = sige_z_PML_2(k)*(be_z_2(k)-1.0)/(sige_z_PML_2(k)+kappae_z_PML_2(k)*alphae_z_PML_2(k))/ kappae_z_PML_2(k);
      end
   end
   for k = 1:nzPML_2-1
      sigh_z_PML_2(k) = sig_z_max * ( (nzPML_2 - k - 0.5)/(nzPML_2-1.0))^m;
      alphah_z_PML_2(k) = alpha_z_max*((k-0.5)/(nzPML_2-1.0))^ma;
      kappah_z_PML_2(k) = 1.0+(kappa_z_max-1.0)*((nzPML_2 - k - 0.5) / (nzPML_2 - 1.0))^m;
      bh_z_2(k) = exp(-(sigh_z_PML_2(k) / kappah_z_PML_2(k) + alphah_z_PML_2(k))*dt/epsO);
      ch_z_2(k) = sigh_z_PML_2(k)*(bh_z_2(k)-1.0)/(sigh_z_PML_2(k)+kappah_z_PML_2(k)*alphah_z_PML_2(k))/ kappah_z_PML_2(k);
   end   
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  FILL IN UPDATING COEFFICIENTS
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   DA = 1.0;
   DB = (dt/(muO));
   for i = 1:Imax
      for j = 1:Jmax
         for k = 1:Kmax
            CA(i,j,k) = (1.0 - sig(i,j,k)*dt / (2.0*eps(i,j,k))) / (1.0 + sig(i,j,k) * dt / (2.0*eps(i,j,k)));
            CB(i,j,k) = (dt/(eps(i,j,k))) / (1.0 + sig(i,j,k)*dt / (2.0*eps(i,j,k)));
         end
      end
   end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  FILL IN DENOMINATORS FOR FIELD UPDATES
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   ii = nxPML_2-1;
   for i = 1:Imax-1
      if (i <= nxPML_1-1)
         den_hx(i) = 1.0/(kappah_x_PML_1(i)*dx);
      elseif (i >= Imax+1-nxPML_2)
         den_hx(i) = 1.0/(kappah_x_PML_2(ii)*dx);
         ii = ii-1;
      else
         den_hx(i) = 1.0/dx;
      end
   end
   jj = nyPML_2-1;
   for j = 1:Jmax-1
      if (j <= nyPML_1-1) 
         den_hy(j) = 1.0/(kappah_y_PML_1(j)*dy);
      elseif (j >= Jmax+1-nyPML_2) 
         den_hy(j) = 1.0/(kappah_y_PML_2(jj)*dy);
         jj = jj-1;
      else
         den_hy(j) = 1.0/dy;
      end
   end
   kk = nzPML_2-1;
   for k = 2:Kmax-1
      if (k <= nzPML_1) 
         den_hz(k) = 1.0/(kappah_z_PML_1(k-1)*dz);
      elseif (k >= Kmax+1-nzPML_2)
         den_hz(k) = 1.0/(kappah_z_PML_2(kk)*dz);
         kk = kk - 1;
      else
         den_hz(k) = 1.0/dz;
      end
   end
 %-------------------  
   ii = nxPML_2;
   for i = 1:Imax-1
      if (i <= nxPML_1)
         den_ex(i) = 1.0/(kappae_x_PML_1(i)*dx);
      elseif (i >= Imax+1-nxPML_2)
         den_ex(i) = 1.0/(kappae_x_PML_2(ii)*dx);
         ii = ii-1;
      else
         den_ex(i) = 1.0/dx;
      end
   end
   jj = nyPML_2;
   for j = 1:Jmax-1
      if (j <= nyPML_1)
         den_ey(j) = 1.0/(kappae_y_PML_1(j)*dy);
      elseif (j >= Jmax+1-nyPML_2)
         den_ey(j) = 1.0/(kappae_y_PML_2(jj)*dy);
         jj = jj-1;
      else

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -