二维FDTD TE波圆柱仿真

% 二维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(0,9*dev,1000);
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(0,s_range,nodes); %在空间划分采样点数
taper = 1/sqrt(2*pi)/dev_s*exp(-1*(spacial - p).^2/2/dev_s^2);
taper = taper./max(taper);

figure(1);
plot(spacial,taper,'c');
title('脉冲源');
xlabel('x轴');
ylabel('Pluse');

% TE波的分量初始化
%----------------------------------------
Ex = zeros(nodes,nodes);
Ey = zeros(nodes,nodes);
Hz = zeros(nodes,nodes);

% 加入金属圆柱
%--------------------------------
center_s = round(nodes/2);
cn = center_s*del_s;
PEC = ones(nodes,nodes); %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(k,j) = 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/eps/del_s;
c_MUR = (c*del_t - del_s)/(c*del_t + del_s);

frames=80;
figure(3);
axis([0 nodes 0 nodes 0 1]);
FDTD_M = moviein(frames);

while (done~=0),
% 初始化
Ex_o = Ex;
Ey_

o = Ey;
Hz_o = Hz;

% 电场值的中心差分

Hterm = c_eps*(Hz(2:nodes-1,1:nodes) - Hz(1:nodes-2,1:nodes));
Ex(2:nodes-1,1:nodes) = Ex(2:nodes-1,1:nodes) + Hterm;

Hterm = c_eps*(Hz(1:nodes,1:nodes-2) - Hz(1:nodes,2:nodes-1));
Ey(1:nodes,2:nodes-1) = Ey(1:nodes,2:nodes-1) + Hterm;

%加入脉冲源
t = n*del_t;
term = (t-mean);
pulse = (-1/sqrt(2*pi)/dev^3)*term;
pulse = pulse* exp((-1/2/dev^2)*term^2) / pulsenorm;
source = pulse * ones(1,nodes) .* taper; % DBD correction

%t = 80; %可选择高斯脉冲
%term = (n-t);
%pulse =exp(-4*pi*term^2/(10^2)) ;
%source = pulse * ones(1,nodes) ;

Ex(2,1:nodes) = Ex(2,1:nodes) + source;

% 4条边界线的处理

% 左边
Ex(1,1:nodes) = Ex_o(2,1:nodes) + c_MUR*(Ex(2,1:nodes) - Ex(1,1:nodes));

% 上边
Ey(1:nodes,nodes) = Ey_o(1:nodes,nodes-1) + c_MUR*(Ey(1:nodes,nodes-1) - Ey(1:nodes,nodes));

% 右边
Ex(nodes,1:nodes) = Ex_o(nodes-1,1:nodes) + c_MUR*(Ex(nodes-1,1:nodes) - Ex(nodes,1:nodes));

% 下边
Ey(1:nodes,1) = Ey_o(1:nodes,2) + c_MUR*(Ey(1:nodes,2) - Ey(1:nodes,1));

% 加入金属圆柱
Ex = Ex.*PEC;
Ey = Ey.*PEC;
Hz = Hz.*PEC;

% 磁场值的中心差分
Eterm1 = c_mu*(Ex(2:nodes,1:nodes-1) - Ex(1:nodes-1,1:nodes-1));
Eterm2 = c_mu*(Ey(1:nodes-1,1:nodes-1) - Ey(1:nodes-1,2:nodes));
Hz(1:nodes-1,1:nodes-1) = Hz(1:nodes-1,1:nodes-1) + Eterm1 + Eterm2;

% 描绘场分布
cut = 3;
figure(3);
clf;
if (abs(round(n/cut) - (n/cut)) == 0), %每5个时间步进行绘图
F = F + 1; %记录绘图的步数
mesh(abs(Ex + i*Ey)); %电场的幅值
axis([0 nodes 0 nodes 0 1]);
FDTD_M(:,F) = getframe;
end;

% 进行时间迭代，done=0时停止计算
if (done ~= 0),
n = n+1; %记录迭代的时间步数
end;
if (n == cut*frames),
done = 0;
end;
end;

figure(4);
mesh(abs(Ex + i*Ey));

% figure(5);
%title('幅值分布图');
%surface(abs(Ex + i*Ey));
%shading flat;

%figure(6);
%title('相位分布图');
%surface(angle(Ex + i*Ey));
%shading interp;
clear all