二维TE波CNDG-FDTD方法matlab实现

CNDG-FDTD方法是一种用于计算二维TE波的有限差分时域方法。下面是一份用MATLAB实现的二维TE波CNDG-FDTD方法的代码：

%参数设置
c = 3e8; %光速
dx = 1e-3; %网格尺寸
dt = dx/c; %时间步长
t = 1e-9; %总时间
f = 10e9; %频率
lambda = c/f; %波长
n = 1; %介质折射率
eps0 = 8.85e-12; %真空介电常数
eps = n^2 * eps0; %介质介电常数
mu0 = 4*pi*1e-7; %真空磁导率
mu = mu0; %介质磁导率

%计算网格数
L = 2*lambda; %计算区域边长
Nx = round(L/dx); %x方向网格数
Ny = round(L/dx); %y方向网格数

%初始化场分量
Ex = zeros(Nx,Ny);
Ey = zeros(Nx,Ny);
Hz = zeros(Nx,Ny);

%计算更新系数
m = 2;
a = (mu*dt)/(dx*m);
b = (eps*dt)/(dx*m);

%计算源
x0 = Nx/2;
y0 = Ny/2;
sigma = Nx/20;
omega = 2*pi*f;
t0 = 6*sigma*dt;
s = exp(-((1:Nx)-x0).^2/(2*sigma^2))' * exp(-((1:Ny)-y0).^2/(2*sigma^2));
Ez = s * sin(omega*0) * sin(pi*t/2/t0)^2;

%进行时间步进
for n = 1:round(t/dt)
    %计算Ex和Ey
    Ex(:,2:Ny) = Ex(:,2:Ny) + b/mu * (Hz(:,2:Ny) - Hz(:,1:Ny-1));
    Ey(2:Nx,:) = Ey(2:Nx,:) - b/mu * (Hz(2:Nx,:) - Hz(1:Nx-1,:));
    
    %计算Hz
    Hz(1:Nx-1,1:Ny-1) = Hz(1:Nx-1,1:Ny-1) + a/eps * ...
        (Ey(1:Nx-1,2:Ny) - Ey(1:Nx-1,1:Ny-1) - ...
         Ex(2:Nx,1:Ny-1) + Ex(1:Nx-1,1:Ny-1));
    
    %加入源
    Hz(:, :) = Hz(:, :) + Ez(:, :);
    
    %绘制场分量
    imagesc(Hz);
    title(sprintf('Time = %g s', n*dt));
    colorbar;
    axis equal tight;
    drawnow;
end

这份代码中，我们使用了一个高斯型源来产生TE波，并使用了CNDG-FDTD方法进行计算。在代码运行结束后，会绘制出Hz分量随时间的变化。