8141阅读
14回复

[求助]发个FDTD程序大家共享

UID ：4992

注册：2007-09-17
登录：2007-09-26
发帖：19
等级：仿真新人

0楼发表于: 2007-09-24 13:26:10

1-D:
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);
end
Es = rot90(Es);
j = 1:Nz;
surf(i,j,Es);
axis([0 Nz 0 Nz -4 4])
otherwise
A = rot90(E);
imagesc(A);
end
pause(0.005);
end

2D:
function FDTD_2D_Yee_2D_TM_final

%constants
c_0 = 3E8; % Speed of light in free space
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

%to be set
Nx = 100; % Number of cells in x-direction
Ny = 100; % Number of cells in y-direction
Nt = 500; % Number of time steps

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
Dy_ref = lambda_ref/20;

X = Nx*Dx_ref;
Y = Ny*Dy_ref;

r = 0.5; % Normalized factor
Dt_ref = r/c_ref*Dx_ref; % Reference time step
Dt = Dt_ref;

% initialization
Ez0 = zeros(Nx, Ny);
Ez1 = zeros(Nx, Ny);
Hx0 = zeros(Nx, Ny);
Hx1 = zeros(Nx, Ny);
Hy0 = zeros(Nx, Ny);
Hy1 = zeros(Nx, Ny);

% Source position
Nx_Source = int16(0.5*(Nx+1));
Ny_Source = int16(0.5*(Nx+1));

pulse = 0;

for n = 1:Nt % Sources' definitions
t = Dt_ref*r*n; % Actual time
Source_type = 1; % Choice of source type

switch Source_type
case 1 % Modified source function
ncycles = 1;
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
otherwise % For debug
pulse = 1;
end
Ez0(Nx_Source,Ny_Source) = Ez0(Nx_Source,Ny_Source) - r*pulse;

CHy = Dt_ref/mu_ref/Dx_ref; % Coefficients used below
CHx = Dt_ref/mu_ref/Dy_ref;
CEzHy = Dt_ref/eps_ref/Dx_ref;
CEzHx = Dt_ref/eps_ref/Dy_ref;

for i = 2:Nx
% H update
Hx1(i,1:Ny-1) = Hx0(i,1:Ny-1) - CHx.*(Ez0(i,2:Ny)-Ez0(i,1:Ny-1));
end
for i = 1:Nx-1
Hy1(i,2:Ny) = Hy0(i,2:Ny) + CHy.*(Ez0(i+1,2:Ny)-Ez0(i,2:Ny));
end

% Boundary conditions **************************************
boundary = 3; % Choice: '1'=Mur ABC; '2'=Dirichlet; '3'=Neumann
switch boundary
case 1 % For H Mur ABC
Hx1(1,1:Ny-1) = Hx0(2, 1:Ny-1) + (r-1)/(r+1).*(Hx1(2, 1:Ny-1)-Hx0(1,1:Ny-1)); % Mur ABC @ left side x = 0
Hx1(1:Nx, Ny) = Hx0(1:Nx,Ny-1) + (r-1)/(r+1).*(Hx1(1:Nx,Ny-1)-Hx0(1:Nx, Ny)); % Mur ABC @ left side y = Ny
Hy1(1:Nx-1,1) = Hy0(1:Nx-1, 2) + (r-1)/(r+1).*(Hy1(1:Nx-1, 2)-Hy0(1:Nx-1,1)); % Mur ABC @ left side y = 0
Hy1(Nx, 1:Ny) = Hy0(Nx-1,1:Ny) + (r-1)/(r+1).*(Hy1(Nx-1,1:Ny)-Hy0(Nx, 1:Ny)); % Mur ABC @ left side x = Nx
case 2 % Dirichlet

Hx1(1:Nx, 1) = 0;
Hx1(1:Nx,Ny) = 0;
Hx1(1, 1:Ny) = 0;
Hx1(Nx,1:Ny) = 0;

Hy1(1:Nx, 1) = 0;
Hy1(1:Nx,Ny) = 0;
Hy1(1, 1:Ny) = 0;
Hy1(Nx,1:Ny) = 0;
case 3 % Neumann
Hx1(1, 1:Ny-1) = Hx0(1, 1:Ny-1);
Hx1(Nx,1:Ny-1) = Hx1(Nx,1:Ny-1);
Hx1(1:Nx, Ny) = Hx0(1:Nx, Ny);
Hx1(1:Nx, 1) = Hx1(1:Nx, 1);

Hy1(1:Nx-1, 1) = Hy0(1:Nx-1, 1);
Hy1(1:Nx-1,Ny) = Hy0(1:Nx-1,Ny);
Hy1(Nx, 1:Ny) = Hy0(Nx, 1:Ny);
Hy1(1, 1:Ny) = Hy0(1, 1:Ny);

end
for i = 2:Nx
% E update
Ez1(i,2:Ny) = Ez0(i,2:Ny) + CEzHy.*(Hy1(i,2:Ny)-Hy1(i-1,2:Ny)) - CEzHx.*(Hx1(i,2:Ny)-Hx1(i,1:Ny-1));
end
switch boundary
case 1 % For E Mur ABC
Ez1(1,1:Ny) = Ez0(2,1:Ny) + (r-1)/(r+1).*(Ez1(2,1:Ny)-Ez0(1,1:Ny)); % Mur ABC @ right side x = 0
Ez1(1:Nx,1) = Ez0(1:Nx,2) + (r-1)/(r+1).*(Ez1(1:Nx,2)-Ez0(1:Nx,1)); % Mur ABC @ right side y = 0
case 2 % Dirichlet
Ez1(1, 1:Ny) = 0;
Ez1(1:Nx, 1) = 0;
Ez1(Nx,1:Ny) = 0;
Ez1(1:Nx,Ny) = 0;
case 3 % Neumann
Ez1(1, 1:Ny) = Ez0(1, 1:Ny);
Ez1(Nx,1:Ny) = Ez0(Nx,1:Ny);
Ez1(1:Nx, 1) = Ez0(1:Nx, 1);
Ez1(1:Nx,Ny) = Ez0(1:Nx,Ny);
end
Hx0 = Hx1;
Hy0 = Hy1;
Ez0 = Ez1;

% Display*********************************************
i = 1:Nx;
j = 1:Ny;
display = 1; % Choice: '1' = Ez, '2' = Hx, '3' = Hy, 'Otherwise' = three components
switch display
case 1
surf(i,j,Ez0);
axis([0 Nx 0 Ny -0.03 0.03]);
set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
xlabel('x in m');
set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
ylabel('y in m');
zlabel('Amplitude of Ez');
case 2
surf(i,j,Hx0);
axis([0 Nx 0 Ny -1e-4 1e-4]);
set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
xlabel('x in m');
set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
ylabel('y in m');
zlabel('Amplitude of Hx');
case 3
surf(i,j,Hy0);
axis([0 Nx 0 Ny -1e-4 1e-4]);
set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
xlabel('x in m');
set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
ylabel('y in m');
zlabel('Amplitude of Hy');
otherwise
subplot(2,2,1);
surf(i,j,Ez0);
title('Ex');
axis([0 Nx 0 Ny -0.03 0.03]);
set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
xlabel('x in m');
set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
ylabel('y in m');
zlabel('Amplitude of Ez');

subplot(2,2,2);
surf(i,j,Hx0);
title('Hx');
axis([0 Nx 0 Ny -1e-4 1e-4]);
set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
xlabel('x in m');
set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
ylabel('y in m');
zlabel('Amplitude of Hx');

subplot(2,2,3);
surf(i,j,Hy0);
title('Hy');
axis([0 Nx 0 Ny -1e-4 1e-4]);
set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
xlabel('x in m');
set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
ylabel('y in m');
zlabel('Amplitude of Hy');
end

%surf(i,j,Ez0);
%axis([0 Nx 0 Ny -0.03 0.03])
%A = rot90(Ez0);
%imagesc(A);
pause(0.005);
end

3D FDTD for hexahedral cavity with conducting walls

% Physical constants
eps0 = 8.8541878e-12; % Permittivity of vacuum
mu0 = 4e-7 * pi; % Permeability of vacuum
c0 = 299792458; % Speed of light in vacuum

% Parameter initiation
Lx = .05; Ly = .04; Lz = .03; % Cavity dimensions in meters
Nx = 25; Ny = 20; Nz = 15; % Number of cells in each direction
Cx = Nx / Lx; % Inverse cell dimensions
Cy = Ny / Ly;
Cz = Nz / Lz;
Nt = 1024; % Number of time steps
Dt = 1/(c0*norm([Cx Cy Cz])); % Time step

% Allocate field matrices
Ex = zeros(Nx , Ny+1, Nz+1);
Ey = zeros(Nx+1, Ny , Nz+1);
Ez = zeros(Nx+1, Ny+1, Nz );
Hx = zeros(Nx+1, Ny , Nz );
Hy = zeros(Nx , Ny+1, Nz );
Hz = zeros(Nx , Ny , Nz+1);

% Allocate time signals
Et = zeros(Nt,3);

% Initiate fields (near but not on the boundary)
Ex(1,2,2) = 1;
Ey(2,1,2) = 2;
Ez(2,2,1) = 3;

% Time stepping
for n = 1:Nt;

% Update H everywhere
Hx = Hx + (Dt/mu0)*(diff(Ey,1,3)*Cz - diff(Ez,1,2)*Cy);
Hy = Hy + (Dt/mu0)*(diff(Ez,1,1)*Cx - diff(Ex,1,3)*Cz);
Hz = Hz + (Dt/mu0)*(diff(Ex,1,2)*Cy - diff(Ey,1,1)*Cx);

% Update E everywhere except on boundary
Ex(:,2:Ny,2:Nz) = Ex(:,2:Ny,2:Nz) + (Dt /eps0) * ...
(diff(Hz(:,:,2:Nz),1,2)*Cy - diff(Hy(:,2:Ny,:),1,3)*Cz);
Ey(2:Nx,:,2:Nz) = Ey(2:Nx,:,2:Nz) + (Dt /eps0) * ...
(diff(Hx(2:Nx,:,:),1,3)*Cz - diff(Hz(:,:,2:Nz),1,1)*Cx);
Ez(2:Nx,2:Ny,:) = Ez(2:Nx,2:Ny,:) + (Dt /eps0) * ...
(diff(Hy(:,2:Ny,:),1,1)*Cx - diff(Hx(2:Nx,:,:),1,2)*Cy);

% Sample the electric field at chosen points
Eto(n,:) = [Ex(4,4,4) Ey(4,4,4) Ez(4,4,4)];
end

Version with for-loops for update of Ex and Hx.

% Physical constants
eps0 = 8.8541878e-12; % Permittivity of vacuum
mu0 = 4e-7 * pi; % Permeability of vacuum
c0 = 299792458; % Speed of light in vacuum

% Parameter initiation
Lx = .05; Ly = .04; Lz = .03; % Cavity dimensions in meters
Nx = 25; Ny = 20; Nz = 15; % Number of cells in each direction
Cx = Nx / Lx; % Inverse cell dimensions
Cy = Ny / Ly;
Cz = Nz / Lz;
Nt = 32; % Number of time steps
Dt = 1/(c0*norm([Cx Cy Cz])); % Time step

% Allocate field matrices
Ex = zeros(Nx , Ny+1, Nz+1);
Ey = zeros(Nx+1, Ny , Nz+1);
Ez = zeros(Nx+1, Ny+1, Nz );
Hx = zeros(Nx+1, Ny , Nz );
Hy = zeros(Nx , Ny+1, Nz );
Hz = zeros(Nx , Ny , Nz+1);

% Allocate time signals
Et = zeros(Nt,3);

% Initiate fields (near but not on the boundary)
Ex(1,2,2) = 1;
Ey(2,1,2) = 2;
Ez(2,2,1) = 3;

% Time stepping
for n = 1:Nt;
% Update H everywhere

for i = 1:Nx+1
for j = 1:Ny
for k = 1:Nz
Hx(i,j,k) = Hx(i,j,k) + (Dt/mu0)* ...
((Ey(i,j,k+1)-Ey(i,j,k))*Cz - (Ez(i,j+1,k)-Ez(i,j,k))*Cy);
end
end
end

Hy = Hy + (Dt/mu0)*((Ez(2:Nx+1,:,:)-Ez(1:Nx,:,:))*Cx ...
- (Ex(:,:,2:Nz+1)-Ex(:,:,1:Nz))*Cz);
Hz = Hz + (Dt/mu0)*((Ex(:,2:Ny+1,:)-Ex(:,1:Ny,:))*Cy ...
- (Ey(2:Nx+1,:,:)-Ey(1:Nx,:,:))*Cx);

% Update E everywhere except on boundary

for i = 1:Nx
for j = 2:Ny
for k = 2:Nz
Ex(i,j,k) = Ex(i,j,k) + (Dt /eps0) * ...
((Hz(i,j,k)-Hz(i,j-1,k))*Cy-(Hy(i,j,k)-Hy(i,j,k-1))*Cz);
end
end
end

Ey(2:Nx,:,2:Nz) = Ey(2:Nx,:,2:Nz) + (Dt /eps0) * ...
((Hx(2:Nx,:,2:Nz)-Hx(2:Nx,:,1:Nz-1))*Cz ...
- (Hz(2:Nx,:,2:Nz)-Hz(1:Nx-1,:,2:Nz))*Cx);
Ez(2:Nx,2:Ny,:) = Ez(2:Nx,2:Ny,:) + (Dt /eps0) * ...
((Hy(2:Nx,2:Ny,:)-Hy(1:Nx-1,2:Ny,:))*Cx ...
- (Hx(2:N ..

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

共1条评分

wwwjjjsss

RFEDA币 +3

优秀资料+1~3

2007-09-24

离线seapeng

淡定，随缘，感恩

UID ：5004