基于FDTD模拟一维磁场 电场

Question

% 定义初始常数
eps_0 = 8.854187817e-12; % 真空介电常数

mu_0 = 4*pi*1e-7; % 自由空间磁导率

c = 1/sqrt(mu_0*eps_0); % 光速
% 定义问题几何和参数
domain_size = 1; % 一维问题空间长度
dx = 1e-3; % 长度单位为米

dt = 3e-12; % 时间步长以秒为单位

number_of_time_steps = 2000; % 迭代次数
nx = round(domain_size/dx); % 一维问题空间中的单元数
source_position = 0.5; % 电流源Jz的位置
% 初始化
Ceze = zeros(nx+1,1);
Cezhy = zeros(nx+1,1);
Cezj = zeros(nx+1,1);
Ez = zeros(nx+1,1);
Jz = zeros(nx+1,1);
eps_r_z = ones (nx+1,1); % 自由空间
sigma_e_z = zeros(nx+1,1); % 自由空间
Chyh = zeros(nx,1);
Chyez = zeros(nx,1);
Chym = zeros(nx,1);
Hy = zeros(nx,1);
My = zeros(nx,1);
mu_r_y = ones (nx,1); % 自由空间
sigma_m_y = zeros(nx,1); % 自由空间
% 计算FDTD更新系数
Ceze = (2 * eps_r_z * eps_0 - dt * sigma_e_z) ...
./(2 * eps_r_z * eps_0 + dt * sigma_e_z);
Cezhy = (2 * dt / dx) ...
./(2 * eps_r_z * eps_0 + dt * sigma_e_z);
Cezj = (-2 * dt) ...
./(2 * eps_r_z * eps_0 + dt * sigma_e_z);
Chyh = (2 * mu_r_y * mu_0 - dt * sigma_m_y) ...
./(2 * mu_r_y * mu_0 + dt * sigma_m_y);
Chyez = (2 * dt / dx) ...
./(2 * mu_r_y * mu_0 + dt * sigma_m_y);
Chym = (-2 * dt) ...
./(2 * mu_r_y * mu_0 + dt * sigma_m_y);
% 定义高斯源波形
time = dt*[0:number_of_time_steps-1].';
Jz_waveform = exp(-((time-2e-10)/5e-11).^2);
source_position_index = round(nx*source_position/domain_size)+1;
% 初始化绘图的子例程
initialize_plotting_parameters;
for time_step = 1:number_of_time_steps
% 更新Jz
Jz(source_position_index) = Jz_waveform(time_step)
% 更新磁场
Hy(1:nx) = Chyh(1:nx) .* Hy(1:nx) ...
+ Chyez(1:nx) .* (Ez(2:nx+1) - Ez(1:nx)) ...
+ Chym(1:nx) .* My(1:nx);
% 更新电场
Ez(2:nx) = Ceze (2:nx) .* Ez(2:nx) ...
+ Cezhy(2:nx) .* (Hy(2:nx) - Hy(1:nx-1)) ...
+ Cezj(2:nx) .* Jz(2:nx);
Ez(1) = 0; % 在x=0m处应用PEC边界条件
Ez(nx+1) = 0; % 在x=1m处应用PEC边界条件
% 当前状态
plot_fields;
end
不知道该怎么改，求大佬指教