我对程序白痴的很，这个是网络上找到的一个程序，如果侵犯了原作者的权益，请通知我，我将及时处理
现在把它分享出来，希望对有需要的人有帮助
-
- % 二维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
二维金属圆柱仿真.m
