Matlab一维FDTD卷积边界条件（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(khT,1); Hy =zeros(khT,1); CA =zeros(khT,1); CB =zeros(khT,1); sig =zeros(khT,1); epsi=zeros(khT,1); %-------------------------------- 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(khT,1);psi_Exz2=zeros(khT,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% psi_Hyz1=zeros(khT,1);psi_Hyz2=zeros(khT,1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% bEz=zeros(kBC,1);cEz=zeros(kBC,1) ; A_EzBC=zeros(kBC,1) ; sigEzBC=zeros(kBC,1) ; K_EzBC=zeros(kBC,1) ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% bHz=zeros(kBC-1,1);cHz=zeros(kBC-1,1) ; A_HzBC=zeros(kBC-1,1); sigHz_BC=zeros(kBC-1); K_HzBC=zeros(kBC-1,1); % Denominators for the update equations F_ez=zeros(keT,1); F_hz=zeros(keT,1); %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ % 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)/ (sigEzBC(k)+K_EzBC(k)*A_EzBC(k)) / K_EzBC(k); end end for k = 1:kBC-1 sigHz_BC(k) = sigZmax * ( (kBC - k - 0.5)/(kBC-1.0))^m; K_HzBC(k) = 1.0+(kZmax-1.0)*((kBC - k - 0.5) / (kBC - 1.0))^m; A_HzBC(k) = aZmax*((k-0.5)/(kBC-1.0))^ma; bHz(k) = exp(-(sigHz_BC(k) / K_HzBC(k) + A_HzBC(k))*dt/eps0); cHz(k) = sigHz_BC(k)*(bHz(k)-1.0)/ (sigHz_BC(k)+K_HzBC(k)*A_HzBC(k)) / K_HzBC(k); end %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ % FILL IN UPDATING COEFFICIENTS %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DA = 1.0; DB = dt/mu0; %------------------------- for k = 1:keT CA(k) = (1.0 - sig(k)*dt / (2.0*epsi(k))) / (1.0 + sig(k) * dt / (2.0*epsi(k))); CB(k) = (dt/epsi(k)) / (1.0 + sig(k)*dt / (2.0*epsi(k))); end %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ % FILL IN DENOMINATORS FOR FIELD UPDATES %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ % CPML_Psi kk = kBC-1; for k = 1:keT if k <= kBC-1 F_hz(k) = 1.0/(K_HzBC(k)*dxyz); elseif k >= khT+1-kBC F_hz(k) = 1.0/(K_HzBC(kk)*dxyz); kk = kk - 1; else F_hz(k) = 1.0/dxyz; end end kk = kBC; for k = 1:keT if k <= kBC F_ez(k) = 1.0/(K_EzBC(k)*dxyz); elseif k >= khT+1-kBC F_ez(k) = 1.0/(K_EzBC(kk)*dxyz); kk = kk - 1; else F_ez(k) = 1.0/dxyz; end end %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ % BEGIN TIME STEP %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ for nTimeStep = 1:MaxTime %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ % UPDATE Ex %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ for k = 2:keT Ex(k) = CA(k) * Ex(k) + CB(k) * (Hy(k-1) - Hy(k))*F_ez(k) ; end %..................................................................... % CPML for bottom Ex, k-direction %..................................................................... for k = 2:kBC psi_Exz1(k) = bEz(k)*psi_Exz1(k) + cEz(k) *(Hy(k-1) - Hy(k))/dxyz; Ex(k) = Ex(k) + CB(k)*psi_Exz1(k); end %..................................................................... % CPML for top Ex, k-direction %..................................................................... kk = kBC; for k = khT+1-kBC:keT psi_Exz2(kk) = bEz(kk)*psi_Exz2(kk) + cEz(kk) *(Hy(k-1) - Hy(k))/dxyz; Ex(k) = Ex(k) + CB(k)*psi_Exz2(kk); kk = kk-1; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %----------------------------------------------------------------------- % SOURCE %----------------------------------------------------------------------- source = exp(-.5*((t0-nTimeStep)/tw)^2.0) ; %source = sin(2*(nTimeStep*dt/t0^2.0)); Ex(ksrc) = Ex(ksrc)+source ; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ % UPDATE Hy %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ for k = 1:keT Hy(k) = DA * Hy(k) + DB * (Ex(k) - Ex(k+1))*F_hz(k) ; end %..................................................................... % CPML for bottom Hy, k-direction %..................................................................... for k = 1:kBC-1 psi_Hyz1(k) = bHz(k)*psi_Hyz1(k) + cHz(k)*(Ex(k) - Ex(k+1))/dxyz; Hy(k) = Hy(k) + DB*psi_Hyz1(k); end %..................................................................... % CPML for top Hy, k-direction %..................................................................... kk = kBC-1; for k = khT+1-kBC:keT psi_Hyz2(kk) = bHz(kk)*psi_Hyz2(kk) + cHz(kk)*(Ex(k) - Ex(k+1))/dxyz; Hy(k) = Hy(k) + DB*psi_Hyz2(kk); kk = kk-1; end %-------------------------------------------------------------------------- %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ % VISUALIZATION %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if mod(nTimeStep,Nplt_jmp)==0 NoFig=NoFig+1 ; timestep=int2str(nTimeStep); %---------------------------------------------------- %---------------------------------------------------- % subplot(3,1,1),plot(Ex); axis([1 keT -1 1]);