matlab写一个三维FDTD，并展示颜色图

以下是一个简单的三维FDTD的MATLAB代码示例，并展示了一个颜色图：

% 三维FDTD示例代码

% 定义常量
c0 = 3e8; % 光速
eps0 = 8.854e-12; % 真空介电常数
mu0 = 4*pi*1e-7; % 真空磁导率

% 定义空间参数
Lx = 1; % x方向长度
Ly = 2; % y方向长度
Lz = 3; % z方向长度
dx = 0.02; % x方向网格间距
dy = 0.02; % y方向网格间距
dz = 0.02; % z方向网格间距
Nx = Lx/dx; % x方向网格数
Ny = Ly/dy; % y方向网格数
Nz = Lz/dz; % z方向网格数
dt = 1/(2*c0/sqrt(1/dx^2+1/dy^2+1/dz^2)); % 时间步长

% 定义电磁场数组
Ex = zeros(Nx+1, Ny, Nz);
Ey = zeros(Nx, Ny+1, Nz);
Ez = zeros(Nx, Ny, Nz+1);
Hx = zeros(Nx, Ny+1, Nz+1);
Hy = zeros(Nx+1, Ny, Nz+1);
Hz = zeros(Nx+1, Ny+1, Nz);

% 定义介质参数
eps_r = ones(Nx, Ny, Nz); % 介电常数
mu_r = ones(Nx, Ny, Nz); % 磁导率

% 定义激励源
f0 = 100e6; % 激励频率
w0 = 2*pi*f0; % 激励角频率
t0 = 3/w0; % 激励峰值时间
sigma = 0.5/w0; % 激励高斯宽度
source = zeros(Nx, Ny, Nz); % 激励源数组
for i=1:Nx
    for j=1:Ny
        for k=1:Nz
            source(i,j,k) = exp(-((i*dx-Lx/2)^2+(j*dy-Ly/2)^2+(k*dz-Lz/2)^2)/sigma^2)*sin(w0*(t0-dt)); % 高斯脉冲
        end
    end
end

% 定义边界条件
boundary = 'pml'; % 边界条件类型
if strcmp(boundary,'pml')
    pml_thickness = 10*dx; % PML厚度
    pml_strength = 1; % PML吸收强度（一般取1）
    pml_m = 4; % PML阶数
    pml_kappa = 1; % PML参数（一般取1）
    pml_sigma_max = -(pml_m+1)*log(pml_strength)/(2*impedance0*pml_thickness); % PML最大电导率
    pml_kappa_max = (pml_m+1)*pml_sigma_max/pml_kappa; % PML最大介电率
    pml_sigma_x = zeros(Nx+1,Ny,Nz); % x方向电导率
    pml_sigma_y = zeros(Nx,Ny+1,Nz); % y方向电导率
    pml_sigma_z = zeros(Nx,Ny,Nz+1); % z方向电导率
    pml_kappa_x = ones(Nx+1,Ny,Nz); % x方向介电率
    pml_kappa_y = ones(Nx,Ny+1,Nz); % y方向介电率
    pml_kappa_z = ones(Nx,Ny,Nz+1); % z方向介电率
    for i=1:pml_thickness/dx
        pml_sigma_x(i,:,:) = pml_sigma_max*(i*dx/(pml_thickness))^pml_m;
        pml_sigma_x(Nx+1-i,:,:) = pml_sigma_max*(i*dx/(pml_thickness))^pml_m;
        pml_kappa_x(i,:,:) = 1 + (pml_kappa_max-1)*((i*dx/(pml_thickness))^pml_m);
        pml_kappa_x(Nx+1-i,:,:) = 1 + (pml_kappa_max-1)*((i*dx/(pml_thickness))^pml_m);
    end
    for j=1:pml_thickness/dy
        pml_sigma_y(:,j,:) = pml_sigma_max*(j*dy/(pml_thickness))^pml_m;
        pml_sigma_y(:,Ny+1-j,:) = pml_sigma_max*(j*dy/(pml_thickness))^pml_m;
        pml_kappa_y(:,j,:) = 1 + (pml_kappa_max-1)*((j*dy/(pml_thickness))^pml_m);
        pml_kappa_y(:,Ny+1-j,:) = 1 + (pml_kappa_max-1)*((j*dy/(pml_thickness))^pml_m);
    end
    for k=1:pml_thickness/dz
        pml_sigma_z(:,:,k) = pml_sigma_max*(k*dz/(pml_thickness))^pml_m;
        pml_sigma_z(:,:,Nz+1-k) = pml_sigma_max*(k*dz/(pml_thickness))^pml_m;
        pml_kappa_z(:,:,k) = 1 + (pml_kappa_max-1)*((k*dz/(pml_thickness))^pml_m);
        pml_kappa_z(:,:,Nz+1-k) = 1 + (pml_kappa_max-1)*((k*dz/(pml_thickness))^pml_m);
    end
elseif strcmp(boundary,'periodic')
    % 周期性边界条件
    Ex(1,:,:) = Ex(Nx,:,:);
    Ex(Nx+1,:,:) = Ex(2,:,:);
    Ey(:,1,:) = Ey(:,Ny,:);
    Ey(:,Ny+1,:) = Ey(:,2,:);
    Ez(:,:,1) = Ez(:,:,Nz);
    Ez(:,:,Nz+1) = Ez(:,:,2);
    Hx(:,1,:) = Hx(:,Ny,:);
    Hx(:,Ny+1,:) = Hx(:,2,:);
    Hy(1,:,:) = Hy(Nx,:,:);
    Hy(Nx+1,:,:) = Hy(2,:,:);
    Hz(:,:,1) = Hz(:,:,Nz);
    Hz(:,:,Nz+1) = Hz(:,:,2);
end

% 开始时间迭代
for n=1:1000
    % 更新H场
    Hx = Hx + (dt/(mu_r*mu0*dx)) .* (Ez(:,:,2:Nz+1)-Ez(:,:,1:Nz));
    Hy = Hy + (dt/(mu_r*mu0*dy)) .* (Ex(:,1:Ny+1,2:Nz+1)-Ex(:,1:Ny+1,1:Nz));
    Hz = Hz + (dt/(mu_r*mu0*dz)) .* (Ey(2:Nx+1,:,1:Nz+1)-Ey(1:Nx,:,1:Nz+1));
    % 处理PML边界条件
    if strcmp(boundary,'pml')
        Hx(:,1:pml_thickness/dy,:) = Hx(:,1:pml_thickness/dy,:) ...
            - (dt/(mu_r(:,1:pml_thickness/dy,:)*mu0*dy)) .* (Ez(:,1:pml_thickness/dy,2:Nz+1)-Ez(:,1:pml_thickness/dy,1:Nz)) ...
            - pml_sigma_x(:,1:pml_thickness/dy,:) .* Hx(:,1:pml_thickness/dy,:);
        Hx(:,Ny+1-pml_thickness/dy+1:Ny,:) = Hx(:,Ny+1-pml_thickness/dy+1:Ny,:) ...
            - (dt/(mu_r(:,Ny+1-pml_thickness/dy+1:Ny,:)*mu0*dy)) .* (Ez(:,Ny-pml_thickness/dy+1:Ny,2:Nz+1)-Ez(:,Ny-pml_thickness/dy+1:Ny,1:Nz)) ...
            - pml_sigma_x(:,Ny-pml_thickness/dy+1:Ny,:) .* Hx(:,Ny-pml_thickness/dy+1:Ny,:);
        Hx(:,:,1:pml_thickness/dz) = Hx(:,:,1:pml_thickness/dz) ...
            - (dt/(mu_r(:,:,1:pml_thickness/dz)*mu0*dz)) .* (Ey(2:Nx+1,:,1:pml_thickness/dz)-Ey(1:Nx,:,1:pml_thickness/dz)) ...
            - pml_sigma_x(:,:,1:pml_thickness/dz) .* Hx(:,:,1:pml_thickness/dz);
        Hx(:,:,Nz+1-pml_thickness/dz+1:Nz) = Hx(:,:,Nz+1-pml_thickness/dz+1:Nz) ...
            - (dt/(mu_r(:,:,Nz+1-pml_thickness/dz+1:Nz)*mu0*dz)) .* (Ey(2:Nx+1,:,Nz-pml_thickness/dz+1:Nz)-Ey(1:Nx,:,Nz-pml_thickness/dz+1:Nz)) ...
            - pml_sigma_x(:,:,Nz-pml_thickness/dz+1:Nz) .* Hx(:,:,Nz-pml_thickness/dz+1:Nz);
        Hy(1:pml_thickness/dx,:,:) = Hy(1:pml_thickness/dx,:,:) ...
            - (dt/(mu_r(1:pml_thickness/dx,:,:)*mu0*dx)) .* (Ez(2:Nx+1,:,1:Nz+1)-Ez(1:Nx,:,1:Nz+1)) ...
            + pml_sigma_y(1:pml_thickness/dx,:,:) .* Hy(1:pml_thickness/dx,:,:);
        Hy(Nx+1-pml_thickness/dx+1:Nx,:,:) = Hy(Nx+1-pml_thickness/dx+1:Nx,:,:) ...
            - (dt/(mu_r(Nx+1-pml_thickness/dx+1:Nx,:,:)*mu0*dx)) .* (Ez(Nx-pml_thickness/dx+1:Nx,:,1:Nz+1)-Ez(Nx-pml_thickness/dx+1:Nx,:,1:Nz+1)) ...
            + pml_sigma_y(Nx-pml_thickness/dx+1:Nx,:,:) .* Hy(Nx-pml_thickness/dx+1:Nx,:,:);
        Hy(:,:,1:pml_thickness/dz) = Hy(:,:,1:pml_thickness/dz) ...
            - (dt/(mu_r(:,:,1:pml_thickness/dz)*mu0*dz)) .* (Ex(:,2:Ny+1,1:pml_thickness/dz)-Ex(:,1:Ny,1:pml_thickness/dz)) ...
            + pml_sigma_y(:,:,1:pml_thickness/dz) .* Hy(:,:,1:pml_thickness/dz);
        Hy(:,:,Nz+1-pml_thickness/dz+1:Nz) = Hy(:,:,Nz+1-pml_thickness/dz+1:Nz) ...
            - (dt/(mu_r(:,:,Nz+1-pml_thickness/dz+1:Nz)*mu0*dz)) .* (Ex(:,2:Ny+1,Nz-pml_thickness/dz+1:Nz)-Ex(:,1:Ny,Nz-pml_thickness/dz+1:Nz)) ...
            + pml_sigma_y(:,:,Nz-pml_thickness/dz+1:Nz) .* Hy(:,:,Nz-pml_thickness/dz+1:Nz);
        Hz(1:pml_thickness/dx,:,:) = Hz(1:pml_thickness/dx,:,:) ...
            - (dt/(mu_r(1:pml_thickness/dx,:,:)*mu0*dx)) .* (Ey(1:Nx+1,2:Ny+1,1:Nz+1)-Ey(1:Nx,2:Ny+1,1:Nz+1)) ...
            - pml_sigma_z(1:pml_thickness/dx,:,:) .* Hz(1:pml_thickness/dx,:,:);
        Hz(Nx+1-pml_thickness/dx+1:Nx,:,:) = Hz(Nx+1-pml_thickness/dx+1:Nx,:,:) ...
            - (dt/(mu_r(Nx+1-pml_thickness/dx+1:Nx,:,:)*mu0*dx)) .* (Ey(Nx-pml_thickness/dx+1:Nx,2:Ny+1,1:Nz+1)-Ey(Nx-pml_thickness/dx+1:Nx,2:Ny+1,1:Nz+1)) ...
            - pml_sigma_z(Nx-pml_thickness/dx+1:Nx,:,:) .* Hz(Nx-pml_thickness/dx+1:Nx,:,:);
        Hz(:,:,1:pml_thickness/dy) = Hz(:,:,1:pml_thickness/dy) ...
            - (dt/(mu_r(:,:,1:pml_thickness/dy)*mu0*dy)) .* (Ex(1:Nx+1,:,2:Nz+1)-Ex(1:Nx+1,:,1:Nz)) ...
            + pml_sigma_z(:,:,1:pml_thickness/dy) .* Hz(:,:,1:pml_thickness/dy);
        Hz(:,:,Nz+1-pml_thickness/dy+1:Nz) = Hz(:,:,Nz+1-pml_thickness/dy+1:Nz) ...
            - (dt/(mu_r(:,:,Nz+1-pml_thickness/dy+1:Nz)*mu0*dy)) .* (Ex(1:Nx+1,:,Nz-pml_thickness/dy+1:Nz)-Ex(1:Nx+1,:,Nz-pml_thickness/dy+1:Nz)) ...
            + pml_sigma_z(:,:,Nz-pml_thickness/dy+1:Nz) .* Hz(:,:,Nz-pml_thickness/dy+1:Nz);
    end
    % 更新E场
    Ex(:,2:Ny,:) = Ex(:,2:Ny,:) + (dt/(eps_r(:,2:Ny,:)*eps0*dy)) .* (Hz(:,2:Ny+1,:)-Hz(:,1:Ny-1,:));
    Ex(:,Ny,:) = Ex(:,Ny,:) + (dt/(eps_r(:,Ny,:)*eps0*dy)) .* (Hz(:,1,:)-Hz(:,Ny,:));
    Ey(2:Nx,:,:) = Ey(2:Nx,:,:) + (dt/(eps_r(2:Nx,:,:)*eps0*dx)) .* (Hz(2:Nx+1,:,:)-Hz(1:Nx-1,:,:));
    Ey(Nx,:,:) = Ey(Nx,:,:) + (dt/(eps_r(Nx,:,:)*eps0*dx)) .* (Hz(1,:,:)-Hz(Nx,:,:));
    Ez(:,:,2:Nz) = Ez(:,:,2:Nz) + (dt/(eps_r(:,:,2:Nz)*eps0*dz)) .* (Hy(:,:,2:Nz+1)-Hy(:,:,1:Nz-1));
    Ez(:,:,Nz) = Ez(:,:,Nz) + (dt/(eps_r(:,:,Nz)*eps0*dz)) .* (Hy(:,:,1)-Hy(:,:,Nz));
    % 添加激励源
    Ez(:,:,1) = Ez(:,:,1) + source(:,:,n);
    % 显示模拟结果
    if mod(n,10)==0
        slice_index = round(Nx/2);
        slice = reshape(sqrt(abs(Ex(slice_index,:,:)).^2+abs(Ey(slice_index,:,:)).^2+abs(Ez(slice_index,:,:)).^2),[Ny,Nz]);
        imagesc(slice');
        axis equal;
        axis off;
        colormap(jet);
        colorbar();
        drawnow;
    end
end