2290阅读
2回复

[讨论]高手来看下这个一维的FDTD的MATLAB仿真

UID ：12369

注册：2008-05-13
登录：2008-06-12
发帖：4
等级：旁观者

0楼发表于: 2008-05-30 15:28:11

function FDTD_debug

%constants
c_0 = 3E8; % Speed of light in free space
Nz = 100; % Number of cells in z-direction
Nt = 200; % Number of time steps
N = Nz; % Global number of cells

E = zeros(Nz,1); % E component
E2 = zeros(Nz,1);
E1 = zeros(Nz,1);
Es = zeros(Nz,Nz); % Es is the output for surf function
%H = zeros(Nz,1); % H component
pulse = 0; % Pulse

%to be set
mu_0 = 4.0*pi*1.0e-7; % Permeability of free space
eps_0 = 1.0/(c_0*c_0*mu_0); % Permittivty of free space

c_ref = c_0; % Reference velocity
eps_ref = eps_0;
mu_ref = mu_0;

f_0 = 10e9; % Excitation frequency
f_ref = f_0; % Reference frequency

omega_0 = 2.0*pi*f_0; % Excitation circular frequency
omega_ref = omega_0;

lambda_ref = c_ref/f_ref; % Reference wavelength

Dx_ref = lambda_ref/20; % Reference cells width
Dz = Dx_ref;
nDz = Dz/Dx_ref;

Z = Nz*Dz;
r = 1; % Normalized time step
Dt_ref = r*Dx_ref/c_ref; % Reference time step
Dt = Dt_ref;

% Source position
Nz_Source = 0.5*Nz;
N_Source = Nz_Source;

for n = 1:Nt-1
t = Dt_ref*r*(n-1); % Actual time

%Source function series
Source_type = 1;
switch Source_type
case 1 % modified source function
ncycles = 2;
if t < ncycles*2.0*pi/(omega_0)
pulse = -0.5*( 1.0 - cos(omega_0*t/ncycles) ) * sin(omega_0*t);
else
pulse = 0;
end
case 2 % sigle cos source function
if t < 2.0*pi/(omega_0)
pulse = 8*c_0^2*Dt_ref^2*mu_ref*omega_0*cos(omega_0*t);
else
pulse = 0;
end
case 3 % Gaussian pulse
if t < Dt_ref*r*50
pulse = -40*c_0*(t-t*25/(n-1))*exp(-(t-t*25/(n-1))^2/2/(50/2.3548)^2/(t/(n-1))^2);
else
pulse = 0;
end
end
E2(N_Source) = E2(N_Source) - r*pulse;
E1 = E2;
E2 = E;

m = 2:Nz-1;
E(m) = r^2*(E2(m+1) - 2*E2(m) + E2(m-1)) + 2*E2(m) - E1(m);

%Boundary
E(1) = 0; E(Nz) = 0; % Dirichlet
%E(1) = E(2);E(Nz) = E(Nz-1); % Neumann
%E(Nz) = E2(Nz-1) + (r-1)/(r+1)*(E(Nz-1)-E2(Nz)); % Mur ABC @ right side z = Nz
%E(1) = E2(2) + (r-1)/(r+1)*(E(2)-E2(1)); % Mur ABC @ left side z = 0

%display
%choice***********************************************
display = 3; % choice: '1' = line plot; '2' = imagesc; '3' = surf
switch display
case 1
i = 1:Nz;
plot(i, E(i));
axis([0 Nz -4 4]);
case 2
A = rot90(E);
imagesc(A);
case 3
i = 1:Nz;
for j = 1:Nz
Es(i,j) = E(i);
&nbs ..

未注册仅能浏览部分内容，查看全部内容及附件请先登录或注册

离线cem-uestc

UID ：9061

注册：2008-03-07
登录：2019-01-05
发帖：2575
等级：荣誉管理员

1楼发表于: 2008-06-01 21:18:24

一维FDTD的最多数据是一维的，还真不能绘制3D的图。
不过可以同时绘制电场、磁场分量的历史曲线，他们在垂直的两个平面，这样可以是3D的

离线cem-uestc

UID ：9061

注册：2008-03-07
登录：2019-01-05
发帖：2575
等级：荣誉管理员

2楼发表于: 2008-06-01 21:20:04

不像1D的FDTD啊，像一维时域波动方程

发帖回复