时域有限差分FDTD分析二维圆柱散射问题

% 二维FDTD  TE波圆柱仿真 

clear all;
close all;
clc;

% 定义常数
%-------------------------
c = 3.0E8;
mu = 1.2566E-6;
eps = 8.8542E-12;

f = 1E9;                                                                      %频率
lambda = c/f;                                                                 %波长
w_max = 2*pi*2*f;

cylinder_circ = lambda*8;                                                    %圆柱的周长，改变圆柱的尺寸 
cylinder_rad = cylinder_circ/2/pi;                                            %圆柱的半径
     
% 定义FDTD网格
%----------------------
del_s = lambda/20;                                                            %每最小波长20个采样点          
del_t = 0.5*del_s/c;                                                          %迭代时间步长

dim_s = 5;
s_range = dim_s * lambda;                     
s_range = ceil（s_range / del_s） * del_s;                                      %计算区域的长度
s_cells = s_range / del_s;
cells = s_cells;                                                              %划分的网格数
nodes = cells+1;                                                              %采样点数

% 定义时间脉冲源
%----------------------------------------
dev_larger = 2;
dev = 1/w_max*dev_larger;

dead = 4;
mean = dev*dead;

t = linspace（09*dev1000）;
term = （t-mean）;
pulse = （-1/sqrt（2*pi）/dev^3）.*term;
pulse = pulse.* exp（（-1/2/dev^2）.*term.^2）;
pulsenorm = max（pulse）;

p = s_range / 2;
dead_s = 1.5;
dev_s = s_range / 4 / dead_s;

spacial = linspace（0s_rangenodes）;                                          %在空间划分采样点数
taper = 1/sqrt（2*pi）/dev_s*exp（-1*（spacial - p）.^2/2/dev_s^2）;
taper = taper./max（taper）;

figure（1）;
plot（spacialtaper‘c‘）;
title（‘脉冲源‘）;
xlabel（‘x轴‘）;
ylabel（‘Pluse‘）;

% TE波的分量初始化
%----------------------------------------
Ex = zeros（nodesnodes）;
Ey = zeros（nodesnodes）;
Hz = zeros（nodesnodes）;

% 加入金属圆柱
%--------------------------------
center_s = round（nodes/2）;
cn = center_s*del_s;
PEC = ones（nodesnodes）;                                                       %nodes*nodes的全1矩阵
for k=1:nodes
  for j=1:nodes
    rad = sqrt（（k*del_s - cn）^2 + （j*del_s - cn）^2）;                           %加入圆柱也可以改为方柱
    if （rad <= cylinder_rad）  
    PEC（kj） = 0;
    end;
  end;
end;

figure（2）;
contour（PEC）;                                                                 %等高绘图
title（‘金属圆柱‘）;

% 计算参数设置
%--------------------------------
done = 1;
n = 0;
F = 0;

c_mu = del_t/mu/del_s;
c_eps = del_t