时域有限差分法仿真一维TE波在分裂场完全匹配层【源码】

%***********************************************************************
% 1-D FDTD code with PML
%***********************************************************************
%
% Program author: Mats Gustafsson
% Department of Electroscience
% Lund Institute of Technology
% Sweden
%
% Date of this version: October 2002
%
%
% Eref Ein Erans
% PEC PML objects PML PEC
% +-----+---+------+---------------------+-----+------+-->
% z=0 z=Lz
% r=1 Nref Nin Ntrans r=Nz+1
%
%***********************************************************************
% Physical constants
mu0 = 4e-7 * pi; % Permeability of vacuum磁导率
c0 = 299792458; % Speed of light in vacuum光速
eps0 = 1/c0^2/mu0; % Permittivity of vacuum电导率
eta0 = sqrt(mu0/eps0); % Wave impedance in vacuum波导率
GHz = 10^9;
mm = 10^(-3);

%***********************************************************************
% Parameter initiation
%***********************************************************************
Lz = 1; % Length in meters米
Tmax = 4*Lz/c0; % Time interval [0,Tmax]时间最大值
Dz = 1*mm; % spatial grid size空间网格尺寸

R = 1/sqrt(3); % stabillity number, grid cells per time step每一单位时间计算的网格单元数
freq = 5*GHz; % center frequency of source excitation激励源的中心频率
fmax = 9*GHz; % Max frequency to plot画图的最大频率
fmin = 1*GHz;

Nz = round(Lz/Dz); % Number of cells in the z-direction z轴上的单元格数量（round 四舍五入）
n_pict = round(Nz/20); % plot each n_pict step 每n_pict步画一次图
linew = 2; % linewidth in plot 画图的线宽

Dt = Dz*R/c0; % Time step 时间步长
Nt = round(Tmax/Dt); % Number of time steps 时间步长数量
Nfft = max(2^(round(log2(Nt)+1)),2*1024); % Number of points for the fft
z = linspace(0,Lz,Nz+1); % z-coordinates
t = Dt*[0:Nt-1]'; % time
lambda = c0/freq; % center wavelength of source excitation
ppw = lambda/Dz % sample points per wave length at the center frequency
omega = 2.0*pi*freq; % angular frequency


Nabc = 1*10; % Number of PML layers to the left and to the right
Nin = 4; % scattered field surface
Nref = 2; % near to far field surface
Ntrans = Nz-2;
Labc = Nabc*Dz; % PML length

%***********************************************************************
% Material parameters
%********************************

***************************************
Nmedia=3; % Number of material
media_par = [1.0 0.0; % vacuum: permitivity and conductivity
3.4 1*10^(-12); %
3.0 0];

objects = 0; % Number of objects
obj_z = Dz/2+[0.500 0.504; % object 1 between z1 and z2
0.624 0.628;
0.24 0.26];
obj_m = [2 2 1]; % Material of objects
obj_c = 'rrb'; % Color of object

epsr = ones(Nz-1,1);
sig = zeros(Nz-1,1);
obj_ind = [];
obj_ind(:,1) = ceil(obj_z(:,1)/Dz);
obj_ind(:,2) = floor(obj_z(:,2)/Dz);
obj_frac(:,1) = obj_ind(:,1)-obj_z(:,1)/Dz;
obj_frac(:,2) = -obj_ind(:,2)+obj_z(:,2)/Dz;
for n=1:objects
epsr(obj_ind(n,1):obj_ind(n,2)) = media_par(obj_m(n),1);
%epsr(obj_ind(n,1)) = obj_frac(n,1) * (media_par(obj_m(n),1)-1)+1;
%epsr(obj_ind(n,2)+1) = obj_frac(n,2) * (media_par(obj_m(n),1)-1)+1;
sig(obj_ind(n,1):obj_ind(n,2)) = media_par(obj_m(n),2);
%sig(obj_ind(n,1)) = obj_frac(n,1) * media_par(obj_m(n),2);
%sig(obj_ind(n,2)+1) = obj_frac(n,2) * media_par(obj_m(n),2);
end


%***********************************************************************
% 分配场矢量
%***********************************************************************
Ex = zeros(Nz+1 ,1);
Hy = zeros(Nz,1);

Eref = zeros(Nt,1);
Etrans = zeros(Nt,1);

% 吸收边界条件
E1abc = zeros(Nabc+1,1); % 左侧
E2abc = zeros(Nabc+1,1); % 右侧
H1abc = zeros(Nabc,1);
H2abc = zeros(Nabc,1);

zpmlH = linspace(0,Labc,Nabc)'/Labc;
zpmlE = linspace(Dz/2,Labc-Dz/2,Nabc-1)'/Labc;
abc_pow = 3;
sigma_max = -(abc_pow+1)*log(10^(-4))/(2*eta0*Labc);
sigmaH2 = sigma_max*zpmlH.^abc_pow*mu0/eps0;
sigmaE2 = sigma_max*zpmlE.^abc_pow;
sigmaE1 = sigmaE2((Nabc-1):-1:1);
sigmaH1 = sigmaH2(Nabc:-1:1);

figure(4)
plot(sigmaH2)

%***********************************************************************
% 波激励
%***********************************************************************
Hdelay = Dt/2-Dz/2/c0;

% 从左侧加入射脉冲
tau = 80.0e-12;
delay = 2.1*tau;
Ein = sin(omega*(t-delay)).*exp(-((t-delay).^2/tau^2));
Hin = 1/eta0*sin(omega*(t-delay+Hdelay)).*exp(-((t-delay+Hdelay).^2/tau^2));

% 谐波从左侧入射的时间
Einmax = max(abs(Ein));

%***********************************************************************
% 时间步进
%***********************************************************************
for n = 1:Nt;

% 从-Dt/2到Dt/2更新磁场
Hy = Hy - (Dt/mu0/Dz) * diff(Ex,1); % 主网格

Hy(Nin) = Hy(Nin) - (Dt/mu0/Dz) * Ein(n);

H1abc = (H1abc.*(mu0/Dt-sigmaH1/2) - diff(E1abc)/Dz) ./ (mu0/Dt+sigmaH1/2); % 从ABC向左
H2abc = (H2abc.*(mu0/Dt-sigmaH2/2) - diff(E2abc)/Dz) ./ (mu0/Dt+sigmaH2/2); % 从ABC向右

% 从0到Dt更新E
Ex(2:Nz) = (Ex(2:Nz).*(eps0*epsr/Dt-sig/2) - diff(Hy,1)/Dz)./(eps0*epsr/Dt+sig/2); % 除

去边界的主网格

Ex(Nin) = Ex(Nin) - Hin(n)/(Dz*eps0/Dt);

Ex(1) = Ex(1) - (Dt /eps0) * (Hy(1)-H1abc(Nabc))/Dz; %左边ABC和主网格之间的边界
E1abc(Nabc+1) = Ex(1); % 补充一点，以简化方案
E1abc(2:Nabc) = (E1abc(2:Nabc).*(eps0/Dt-sigmaE1/2) - diff(H1abc,1)/Dz) ./ (eps0/Dt+sigmaE1/2);

Ex(Nz+1) = Ex(Nz+1) - (Dt /eps0) * (H2abc(1)-Hy(Nz))/Dz;
E2abc(1) = Ex(Nz+1);
E2abc(2:Nabc) = (E2abc(2:Nabc).*(eps0/Dt-sigmaE2/2) - diff(H2abc,1)/Dz) ./ (eps0/Dt+sigmaE2/2);

% 在选定点采样电场
if mod(n,n_pict) == 0
figure(1);
subplot(2,1,1);
plot(z,Ex,'LineWidth',1);
title(['time is ',num2str(n*Dt*10^9),'ns'])
axis tight;
ax=axis;
axis([ax(1:2) -1.2*Einmax 1.2*Einmax]);
ax=axis;
ylabel('E-field');
for p=1:objects
line([obj_z(p,1) obj_z(p,1)], ax(3:4),'color',obj_c(p));
line([obj_z(p,2) obj_z(p,2)], ax(3:4),'color',obj_c(p));
line([obj_z(p,1) obj_z(p,2)], [ax(3) ax(3)],'color',obj_c(p),'LineWidth',4);
line([obj_z(p,1) obj_z(p,2)], [ax(4) ax(4)],'color',obj_c(p),'LineWidth',4);
end
subplot(2,1,2);
plot([E1abc(:)' Ex' E2abc(:)'],'LineWidth',1)
title(['time is ',num2str(n*Dt*10^9),'ns'])
ax = axis;
line([Nabc Nabc], ax(3:4),'color','r');
line(Nin+[Nabc Nabc], ax(3:4),'color','b');
line(Ntrans+2+[Nabc Nabc], ax(3:4),'color','b');
line([Nz+2+Nabc Nz+2+Nabc], ax(3:4),'color','r');
axis tight;
ylabel('E-field');
%plot(Eabc(2:Nabc,1:2));
pause(0.01);
end
Eref(n) = Ex(Nref);
Etrans(n) = Ex(Ntrans);
end

%***********************************************************************
% 后处理
%***********************************************************************
Ehin = fft(Ein,Nfft)/Nt*2;
Df = 1/(Dt*Nfft);
freqN = round(fmax/Df);
freqN1 = round(fmin/Df);
f = linspace(Df*freqN1,Df*freqN,freqN-freqN1+1)/10^9;
t = t*10^9; % 纳米秒（ns）

figure(2);
% 入射波
clf;
subplot(2,3,1);
plot(t,Ein,'r','LineWidth',linew)
axis tight;
xlabel('time');
ylabel('E-field');

subplot(2,3,4);
warning off;
A=plotyy(f,abs(Ehin(freqN1:freqN)),f,c0./f/Dz/10^9);
warning on;
set(A(2),'yLim',[0 40],'yTick',[0:5:40],'xlim',[fmin fmax]/10^9,'xGrid','on','yGrid','on','LineWidth',linew);
set(A(1),'xlim',[fmin fmax]/10^9,'LineWidth',linew);
set(get(A(2),'Ylabel'),'String','points/wavelength');
set(get(A(1),'Ylabel'),'String','input');
set(get(A(1),'Xlabel'),'String','frequency');

% 反射
subplot(2,3,2);
plot(t,Eref,'LineWidth',linew); axis tight;
xlabel('time [ns]');
ylabel('Reflected field');
subplot(2,3,5);
Ehref = fft(Eref,Nfft)/Nt*2;
plot(f,20*log10(abs(Ehref(freqN1:freqN)./Ehin(freqN1:freqN))),'LineWidth',linew); axis tight;
ax=axis;
axis([fmin/10^9 fmax/10^9 ax(3:4)]);
axis tight;
xlabel('frequency [GHz]');
ylabel('Refle

ction in dB');

% 传播
subplot(2,3,3);
plot(t,Etrans,'LineWidth',linew); axis tight;
xlabel('time [ns]');
ylabel('Transmitted field');
subplot(2,3,6);
Ehtrans = fft(Etrans,Nfft)/Nt*2;
plot(f,20*log10(abs(Ehtrans(freqN1:freqN)./Ehin(freqN1:freqN))),'LineWidth',linew); axis tight;
ax=axis;
axis([fmin/10^9 fmax/10^9 ax(3:4)]);
xlabel('frequency [GHz]');
ylabel('Transmission in dB');

tr=interp1(f,abs(Ehtrans(freqN1:freqN)./Ehin(freqN1:freqN)),5)
re=interp1(f,abs(Ehref(freqN1:freqN)./Ehin(freqN1:freqN)),5)
tr^2+re^2