背景：
基於公式1.42（Ez分量）、1.43（Hy分量）的1D FDTD實現。
計算電場和磁場分量，該分量由z方向的電流片Jz產生，Jz位於兩個理想導體極板中間，兩個極板平行且向y和z方向無限延伸。
平行極板相距1m，差分網格Δx=1mm。
電流面密度導致分界面（電流薄層）磁場分量的不連續，在兩側產生Hy的波，每個強度為5×10^（-4）A/m。因為電磁波在自由空間中傳播，本征阻抗為η0。
顯示了從左右極板反射前后的傳播過程。本例中是PEC邊界，正切電場分量（Ez）在PEC表面消失。
觀察圖中step650到700的場的變換情況，經過PEC板的入射波和反射波的傳播特征。經過PEC反射后，Ez的極性變反，這是因為反射系數等於-1；
而磁場分量Hy沒有反向，反射系數等於1。
下面是書中的代碼（幾乎沒動）：
第1個是主程序：
%% ------------------------------------------------------------------------------
%% Output Info about this m-file
fprintf('\n****************************************************************\n');
fprintf('\n <FDTD 4 ElectroMagnetics with MATLAB Simulations> \n');
fprintf('\n Listing A.1 \n\n');
time_stamp = datestr(now, 31);
[wkd1, wkd2] = weekday(today, 'long');
fprintf(' Now is %20s, and it is %7s \n\n', time_stamp, wkd2);
%% ------------------------------------------------------------------------------
% This program demonstrates a one-dimensional FDTD simulation.
% The problem geometry is composed of two PEC plates extending to
% infinity in y and z dimensions, parallel to each other with 1 meter
% separation. The space between the PEC plates is filled with air.
% A sheet of current source paralle to the PEC plates is placed
% at the center of the problem space. The current source excites fields
% in the problem space due to a z-directed current density Jz,
% which has a Gaussian waveform in time.
% Define initial constants
eps_0 = 8.854187817e-12; % permittivity of free space
mu_0 = 4*pi*1e-7; % permeability of free space
c = 1/sqrt(mu_0*eps_0); % speed of light
% Define problem geometry and parameters
domain_size = 1; % 1D problem space length in meters
dx = 1e-3; % cell size in meters, Δx=0.001m
dt = 3e-12; % duration of time step in seconds
number_of_time_steps = 2000; % number of iterations
nx = round(domain_size/dx); % number of cells in 1D problem space
source_position = 0.5; % position of the current source Jz
% Initialize field and material arrays
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); % free space
sigma_e_z = zeros(nx+1, 1); % free space
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); % free space
sigma_m_y = zeros(nx, 1); % free space
% Calculate FDTD updating coefficients
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);
% Define the Gaussian source waveform
time = dt * [0:number_of_time_steps-1].';
Jz_waveform = exp(-((time-2e-10)/5e-11).^2)*1e-3/dx;
source_position_index = round(nx * source_position/domain_size)+1;
% Subroutine to initialize plotting
initialize_plotting_parameters;
% FDTD loop
for time_step = 1:number_of_time_steps
% Update Jz for the current time step
Jz(source_position_index) = Jz_waveform(time_step);
% Update magnetic field
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);
% Update electric field
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; % Apply PEC boundary condition at x = 0 m
Ez(nx+1) = 0; % Apply PEC boundary condition at x = 1 m
% Subroutine to plot the current state of the fields
plot_fields;
end
第2個是initialize_plotting_parameters，看名字就知道是初始化參數：
% Subroutine used to initialize 1D plot
Ez_positions = [0:nx]*dx;
Hy_positions = ([0:nx-1]+0.5)*dx;
v = [0 -0.1 -0.1; 0 -0.1 0.1; 0 0.1 0.1; 0 0.1 -0.1; ...
1 -0.1 -0.1; 1 -0.1 0.1; 1 0.1 0.1; 1 0.1 -0.1];
f = [1 2 3 4; 5 6 7 8];
axis([0 1 -0.2 0.2 -0.2 0.2]);
lez = line(Ez_positions, Ez*0, Ez, 'Color', 'b', 'linewidth', 1.5);
lhy = line(Hy_positions, 377*Hy, Hy*0, 'Color', 'r', 'LineWidth', 1.5, 'LineStyle','-.');
set(gca, 'fontsize', 12, 'fontweight', 'bold');
set(gcf,'Color','white');
axis square;
legend('E_{z}', 'H_{y} \times 377', 'location', 'northeast');
xlabel('x [m]');
ylabel('[A/m]');
zlabel('[V/m]');
grid on;
p = patch('vertices', v, 'faces', f, 'facecolor', 'g', 'facealpha', 0.2);
text(0, 1, 1.1, 'PEC', 'horizontalalignment', 'center', 'fontweight', 'bold');
text(1, 1, 1.1, 'PEC', 'horizontalalignment', 'center', 'fontweight', 'bold');
第3個就是畫圖：
% Subroutine used to plot 1D transient field delete(lez); delete(lhy); lez = line(Ez_positions, Ez*0, Ez, 'Color', 'b', 'LineWidth', 1.5); lhy = line(Hy_positions, 377*Hy, Hy*0, 'Color', 'r', 'LineWidth', 1.5, 'LineStyle', '-.'); ts = num2str(time_step); ti = num2str(dt*time_step*1e9); title(['time step = ' ts ' , time = ' ti ' ns']); drawnow;
運行結果：
上圖，從PEC板反射后，電場分量極性變反，磁場分量極性不變。
反射后，電場分量極性再次改變；
