FDTD_2D_CPML_S2.zip_2D FDTD GPR_GPR_GPR-FDTD-2D_cpml_gpr fdtd

% TM simulation for the FDTD method (CPML) % 此版本代码alpha在外边界处衰减至0 %----------------------------------------------------------------------% tic; close all; clear; clc; % MOVIE = VideoWriter('CPML Movie.avi'); %----------------------------------------------------------------------% c0 = 2.99792458e8; % 光速 eps0 = 8.854187818e-12; % 真空中的介电常数 mu0 = 4*pi*1e-7; % 真空中的磁导率 %----------------------------------------------------------------------% order = 4; % PML层中的指数系数 thickness = 10; % 10比8效果更佳 ij=100; % 模拟区域边长 xdim = ij+2*thickness; % x方向（纵向）总长 ydim = ij+2*thickness; % y方向（横向）总长 %----------------------------------------------------------------------% freq = 4e8; % 主频 xexcitation = xdim/2; % 源加载的位置 yexcitation =ydim/2; % 源加载的位置 %----------------------------------------------------------------------% % 模拟区域物性参数 eps1=3.5*ones(xdim,ydim); % 相对介电常数 % eps1(50:60,100:120)=8; % eps1(140:160,1:end)=32; sigma1=0.0025*ones(xdim,ydim); % 电导率 %0.05在减少噪声方面比0.005更佳 mu1= ones(xdim,ydim); % 相对磁导率 %----------------------------------------------------------------------% velocityfield=c0./sqrt(eps1(1:xdim,1:ydim).*mu1(1:xdim,1:ydim)); c_min=min(min(velocityfield)); %介质中最小波速 c_max=max(max(velocityfield)); %介质中最大波速 %----------------------------------------------------------------------% n_tstep = 2000; % 时间总步长 dt = 2.5e-11; % dt %----------------------------------------------------------------------% source66=2; %----------------------------------------------------------------------% if source66==1 %-------------------------------雷克子波-------------------------------% omega=2*pi*freq;%omega alfa=0.93*omega;%alfa x=0:1:n_tstep-1; t=x*dt; source=(t).*(t).*exp(-alfa*(t)).*sin(omega*(t)); source=source/max(abs(source)); %----------------------------------------------------------------------% nn = 2^nextpow2(n_tstep); W = abs(fftshift(fft(source,nn))); W = W./max(W); W = W(nn/2+1:end); fn = 0.5/dt; df = 2*fn/nn; f = -fn:df:fn-df; f = f(nn/2+1:end); threshold=0.02; fmax = f(find(W>=threshold, 1, 'last' )); wlmin = c_min/(fmax); ds1=wlmin/12; %ds<ds1 %----------------------------------------------------------------------% end if source66==2 %------------------------------哈里斯脉冲------------------------------% % blackharrispulse % 哈里斯窗的导数 % Note that their formulation has been changed here to T = 1.14/fr, % such that fr represents the approximate dominant frequency of the resulting pulse. % compute the Blackman-Harris window as specified in Chen et al. (1997) x=0:1:n_tstep-1; dss=c_min/freq/20; t=x*dt; a = [0.35322222 -0.488 0.145 -0.010222222]; % T = 1.14/freq; T = 1.55/freq; window = zeros(1,n_tstep); for n=0:3 window = window + a(n+1)*cos(2*n*pi*t./T); end window(t>=T) = 0; % for the pulse, approximate the window's derivative（导数） and normalize（归一化） blackharrispulse = [window(2:end) 0] - window(1:end); % 导数|差分 blackharrispulse = blackharrispulse./max(abs(blackharrispulse)); %归一化 % %----------------------------------------------------------------------% nn = 2^nextpow2(n_tstep); W = abs(fftshift(fft(blackharrispulse,nn))); W = W./max(W); W = W(nn/2+1:end); fn = 0.5/dt; df = 2*fn/nn; f = -fn:df:fn-df; f = f(nn/2+1:end); threshold=0.02; % 阈值 fmax = f(find(W>=threshold, 1, 'last' )); wavelengthmin = c_min/(fmax); ds1=wavelengthmin/12; % ds<ds1 %----------------------------------------------------------------------% end %----------------------------------------------------------------------% ds = 0.99*ds1; % ds<ds1 dt1= ds/(sqrt(2)*c_max); % dt<dt1 % dt=sourcedt; if dt>dt1 msgbox( 'dt应小于dt1 ','出错!'); end %----------------------------------------------------------------------% xekappa = ones(1,xdim); % grading of kappa on the electric grid in x-direction xhkappa = ones(1,xdim); % grading of kappa on the magnetic grid in x-direction ixekappa = ones(1,xdim); % grading of kappa on the electric grid in x-direction ixhkappa = ones(1,xdim); % grading of kappa on the magnetic grid in x-direction xesigma = zeros(1,xdim); % grading of sigma on the electric grid in x-direction xhsigma = zeros(1,xdim); % grading of sigma of the magnetic grid in x-direction xealpha = zeros(1,xdim); % alpha就是课本的a_w xhalpha = zeros(1,xdim); % grading of alpha on the magnetic grid in x-direction yekappa = ones(1,ydim); % grading of kappa on the electric grid in y-direction yhkappa = ones(1,ydim); % grading of kappa on the magnetic grid in y-direction iyekappa = ones(1,ydim); % grading of kappa on the electric grid in y-direction iyhkappa = ones(1,ydim); % grading of kappa on the magnetic grid in y-direction yesigma = zeros(1,ydim); % grading of sigma on the electric grid in y-direction yhsigma = zeros(1,ydim); % grading of sigma on the magnetic grid in y-direction yealpha = zeros(1,ydim); % grading of alpha on the electric grid in y-direction yhalpha = zeros(1,ydim); % grading of alpha on the magnetic grid in y-direction %-------------------------------------------------------------------------------------------------% sigma_max = (order+1)/(150*sqrt(eps1(thickness+1,thickness+1))*pi*ds); % PML层的最大电阻率(在有些参考书上前面有乘以0.7-1.5的系数) kappa_max = 5; % kappa最大值 ,kappa值的引入是为了改善PML对表面波的吸收特性 alpha_max = 0.005; % 0-0.05， alpha值是为了改善PML对低频分量的吸收特性。 %-------------------------------------------------------------------------------------------------% % 1:10 for i=1:thickness xesigma(1,i)=sigma_max*(((thickness+1)-i)/thickness)^order; % 代入后其实就是(11-i)*ds xekappa(1,i)=1+(kappa_max-1)*(((thickness+1)-i)/thickness)^order; xealpha(1,i)=alpha_max*((i-1)/(thickness-1)); %简简单单的等比数列 end % 51:60 for i=xdim-thickness+1:xdim % (xdim-thickness-1)*ds=(xdim-thickness-1)*ds ; 故i=xdim-thickness也相当于从xdim-thickness+1 ; xesigma(1,i)=sigma_max*((i-(xdim-thickness))/thickness)^order; % 代入后其实就是(i-50)*ds xekappa(1,i)=1+(kappa_max-1)*((i-(xdim-thickness))/thickness)^order; xealpha(1,i)=alpha_max*((xdim-i)/(thickness-1)); end % 1:10 for i=1:thickness % xhsigma(1,i)=sigma_max*((i-(thickness+0.5))/thickness)^order; % (10.5-i)*ds xhkappa(1,i)=1+(kappa_max-1)*((i-(thickness+0.5))/thickness)^order; xhalpha(1,i)=alpha_max*((i-0.5)/(thickness-1)); end % 50:59 for i=xdim-thickness:xdim-1 % (xdim-thickness-1)*ds = (xdim-thickness-1)*ds; xhsigma(1,i)=sigma_max*((i-(xdim-thickness-0.5))/thickness)^order; %(49.5-i)*ds xhkappa(1,i)=1+(kappa_max - 1)*((i-(xdim-thickness-0.5))/thickness)^order; xhalpha(1,i)=alpha_max*((xdim-0.5-i)/(thickness-1)); end for i=1:xdim ixekappa(1,i)=1/xekappa(1,i); ixhkappa(1,i)=1/xhkappa(1,i); end %-------------------------------------------------------------------------------------------------% %1:10 for j=1:thickness yesigma(1,j)=sigma_max*(((thickness+1)-j)/thickness)^order; yekappa(1,j)=1+(kappa_max-1)*(((thickness+1)-j)/thickness)^order; yealpha(1,j)=alpha_max*((j-1)/(thickness-1)); end %51:60 for j=ydim-thickness+1:ydim yesigma(1,j)=sigma_max*((j-(ydim-thickness))/thickness)^order; yekappa(1,j)=1+(kappa_max-1)*((j-(ydim-thickness))/thickness)^order; yealpha(1,j)=alpha_max*((ydim-j)/(thickness-1)); end %1:10 for j=1:thickness yhsigma(1,j)=sigma_max*(((thickness+0.5)-j)/thickness)^order; yhkappa(1,j)=1+(kappa_max-1)*(((thickness+0.5)-j)/thickness)^order; yhalpha(1,j)=alpha_max*((j-0.5)/(thickness-1)); end %50:59 for j=ydim-thickness:ydim-1 yhsigma(1,j)=sigma_ma