Matlab一维FDTD卷积边界条件CPML

处理倏逝波，常规的PML会产生大的误差，这就需要卷积边界条件（CPML）

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% One-Dimensional_FDTD CPML ABC
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
clear;
clear all;
clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Total number of time steps
MaxTime  = 130;
%--------------------------------    
epsR = 1.0;
sigM1 = 0.0 ; 
%--------------------------------    
c0 = 2.99792458E8;
mu0 = 4.0 * pi * 1.0E-7; 
eps0 = 1.0/（c0*c0*mu0）;  
%--------------------------------  
freq=1.9341e+14;
lambda=c0/freq;         
%--------------------------------  
dxyz = .25E-3;
%dxyz = lambda/8;
dt =   dxyz/（c0）;
%--------------------------------  
% Specify the Impulsive Source
 %  tw = 53.0E-12; tO = 4.0*tw;
   tw=6; t0 = 20;
% PML thickness in each direction
kBC = 40; 
%--------------------------------  
% Size of main grid 
ke=120; 
kh=ke+1;   
% Size of total computational domain      
keT=ke+2*kBC;        
khT=keT+1;
%--------------------------------  
Ex =zeros（khT1）;  
Hy =zeros（khT1）; 
CA =zeros（khT1）;
CB =zeros（khT1）;
sig =zeros（khT1）;
epsi=zeros（khT1）;
%--------------------------------  
ksrc=round（keT/2）; 
%--------------------------------  
% Specify the CPML Order and Other Parameters 
m = 3; ma = 1 ;
sigZmax =  （0.8*（m+1）/（dxyz*（mu0/eps0*epsR）^0.5））;
aZmax = 0.05;
kZmax = 1.0;
%--------------------------------      
record_grid = MaxTime;
NoFig=99;
Nplt_jmp=5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
psi_Exz1=zeros（khT1）;psi_Exz2=zeros（khT1）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
psi_Hyz1=zeros（khT1）;psi_Hyz2=zeros（khT1）;       
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
bEz=zeros（kBC1）;cEz=zeros（kBC1） ;
A_EzBC=zeros（kBC1） ;
sigEzBC=zeros（kBC1） ;
K_EzBC=zeros（kBC1） ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bHz=zeros（kBC-11）;cHz=zeros（kBC-11） ;
A_HzBC=zeros（kBC-11）;
sigHz_BC=zeros（kBC-1）;
K_HzBC=zeros（kBC-11）;
% Denominators for the update equations 
F_ez=zeros（keT1）;
F_hz=zeros（keT1）;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  INITIALIZE VARIABLES
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
sig（:） = sigM1;
epsi（:） = epsR*eps0;  
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  SET CPML PARAMETERS IN EACH DIRECTION
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%  z-dir
%
   for k = 1:kBC
      sigEzBC（k） = sigZmax * （ （kBC - k ） / （kBC - 1.0） ）^m;
      K_EzBC（k） = 1.0+（kZmax-1.0）*（（kBC - k） / （kBC - 1.0））^m;
      A_EzBC（k） = aZmax*（（k-1）/（kBC-1.0））^ma;
      bEz（k） = exp（-（sigEzBC（k） / K_EzBC（k） + A_EzBC（k））*dt/eps0）;
      if （（sigEzBC（k） == 0.0） &&    ...
         （A_EzBC（k） == 0.0） && （k == kBC）） 
         cEz（k） = 0.0;
      else
         cEz（k） = sigEzBC（k）*（bEz（k）-1.0）/ （sigE

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件      43056  2008-12-30 00:56  CPML_1D\CPML1D_055.jpg

     文件      18897  2008-12-30 00:56  CPML_1D\CPML1D_105.jpg

     文件      41579  2008-12-30 00:55  CPML_1D\CPML1D_130.jpg

     文件       8633  2008-12-30 00:53  CPML_1D\CPML_1D.asv

     文件       8628  2008-12-30 01:00  CPML_1D\CPML_1D.m

     目录          0  2009-01-02 23:51  CPML_1D

----------- ---------  ---------- -----  ----

               121011                    7