二维TE波紧差分格式的ADI-FDTD方法，采用高斯源，引入PML吸收边界条件matlab实现

二维TE波紧差分格式的ADI-FDTD方法，采用高斯源，引入PML吸收边界条件的matlab实现如下：

clc
clear all
close all

%% 参数设置
c = 3e8;         % 光速
f = 10e9;        % 频率
lambda = c/f;    % 波长
dx = lambda/20;  % 空间步长
dy = lambda/20;  % 空间步长
dt = dx/(2*c);   % 时间步长
T = 500*dt;      % 总时间
m = round(T/dt);  % 总时间步数
N = 80;          % 空间离散点数
M = 80;          % 空间离散点数

%% 材料参数设置
eps0 = 8.85e-12;  % 真空介电常数
epsr = 4;         % 相对介电常数
eps = eps0 * epsr;
sigma = 0.01;     % 电导率
mu0 = 1.26e-6*pi; % 真空磁导率
mur = 1;          % 相对磁导率
mu = mu0 * mur;

%% 初始化电场和磁场
Ex = zeros(N, M);
Ey = zeros(N, M);
Hz = zeros(N, M);

%% 计算PML参数
npml = 20;       % PML层数
Rmax = 1e-8;     % 最大反射系数
m = 4;           % PML阶数
sigmax = -(m+1)*eps0*c*log(Rmax)/(2*dx*npml);  % PML电导率
kappamax = 1;    % PML介质导电率
kappa = 1 + (kappamax-1)*((1:npml)/npml).^m;    % PML介质导电率
sigma_x = zeros(N, M);
sigma_y = zeros(N, M);
for i = 1:npml
    sigma_x(:,i) = sigmax*(npml-i+0.5)/npml^m;
    sigma_x(:,M-i+1) = sigmax*(npml-i+0.5)/npml^m;
    sigma_y(i,:) = sigmax*(npml-i+0.5)/npml^m;
    sigma_y(N-i+1,:) = sigmax*(npml-i+0.5)/npml^m;
end
kappa_x = ones(N,M);
kappa_y = ones(N,M);
for i = 1:npml
    kappa_x(:,i) = 1 + (kappa(i)-1)*((npml-i+0.5)/npml).^m;
    kappa_x(:,M-i+1) = 1 + (kappa(i)-1)*((npml-i+0.5)/npml).^m;
    kappa_y(i,:) = 1 + (kappa(i)-1)*((npml-i+0.5)/npml).^m;
    kappa_y(N-i+1,:) = 1 + (kappa(i)-1)*((npml-i+0.5)/npml).^m;
end

%% 计算系数
c1 = (1 - sigma_x*dt./(2*eps*kappa_x))./(1 + sigma_x*dt./(2*eps*kappa_x));
c2 = dt./(eps*kappa_x*dx);
c3 = (1 - sigma_y*dt./(2*eps*kappa_y))./(1 + sigma_y*dt./(2*eps*kappa_y));
c4 = dt./(eps*kappa_y*dy);
c5 = (1 - sigma*dt./(2*mu))./(1 + sigma*dt./(2*mu));
c6 = dt./(mu*dx);
c7 = dt./(mu*dy);

%% 高斯源参数
s0 = 1;
x0 = N/2*dx;
y0 = M/2*dy;
sigmax = lambda/4;
sigmay = lambda/4;

%% 进行时间步进
for n = 1:m
    % 更新Hz场
    Hz(:,2:M) = c5.*Hz(:,2:M) + c6.*(Ex(:,2:M)-Ex(:,1:M-1)) - c7.*(Ey(2:N,:)-Ey(1:N-1,:));
    
    % 更新Ex场
    Ex(2:N,:) = c1(2:N,:).*Ex(2:N,:) + c2(2:N,:).*(Hz(2:N,:)-Hz(1:N-1,:)) - c3(2:N,:).*sigma_y(2:N,:).*Ey(2:N,:);
    Ex(:,1) = 0;
    Ex(:,M) = 0;
    
    % 更新Ey场
    Ey(:,2:M) = c3(:,2:M).*Ey(:,2:M) + c4(:,2:M).*(Hz(:,2:M)-Hz(:,1:M-1)) - c1(:,2:M).*sigma_x(:,2:M).*Ex(:,2:M);
    Ey(1,:) = 0;
    Ey(N,:) = 0;
    
    % 添加高斯源
    t = n*dt;
    f = exp(-((t-5e-9)/sigmax)^2)*s0*dy;
    g = exp(-((t-5e-9)/sigmay)^2)*s0*dx;
    Ey(round(x0/dx),round(y0/dy)) = Ey(round(x0/dx),round(y0/dy)) + f;
    Ex(round(x0/dx),round(y0/dy)) = Ex(round(x0/dx),round(y0/dy)) + g;
    
    % 绘图
    imagesc(Hz');
    colormap('jet');
    colorbar;
    title(sprintf('Time = %.2e s', t));
    xlabel('x');
    ylabel('y');
    pause(0.01);
end

其中，PML吸收边界条件的处理方法是在边界处引入一个带有吸收特性的介质，通过调整其电导率和介质导电率来达到吸收边界的效果。在此代码中，我们采用了一阶的PML吸收边界条件，并使用了梯形规则进行时间离散化。同时，为了模拟高斯源的情况，我们在代码中添加了一个高斯源，并在每个时间步内进行更新。最终，我们使用`imagesc`函数将模拟结果进行可视化展示。