2265阅读
1回复

[讨论]我们一起来讨论1D FDTD （附代码）

UID ：13713

注册：2008-06-12
登录：2008-07-11
发帖：12
等级：仿真新人

0楼发表于: 2008-06-17 11:50:16

%***********************************************************************
% 1-D FDTD code with simple radiation boundary conditions
%***********************************************************************
%
% Program author: Susan C. Hagness
% Department of Electrical and Computer Engineering
% University of Wisconsin-Madison
% 1415 Engineering Drive
% Madison, WI 53706-1691
% 608-265-5739
% hagness@engr.wisc.edu
%
% Date of this version: February 2000
%
% This MATLAB M-file implements the finite-difference time-domain
% solution of Maxwell's curl equations over a one-dimensional space
% lattice comprised of uniform grid cells.
%
% To illustrate the algorithm, a sinusoidal wave (1GHz) propagating
% in a nonpermeable lossy medium (epsr=1.0, sigma=5.0e-3 S/m) is
% modeled. The simplified finite difference system for nonpermeable
% media (discussed in Section 3.6.6 of the text) is implemented.
%
% The grid resolution (dx = 1.5 cm) is chosen to provide 20
% samples per wavelength. The Courant factor S=c*dt/dx is set to
% the stability limit: S=1. In 1-D, this is the "magic time step."
%
% The computational domain is truncated using the simplest radiation
% boundary condition for wave propagation in free space:
%
% Ez(imax,n+1) = Ez(imax-1,n)
%
% To execute this M-file, type "fdtd1D" at the MATLAB prompt.
% This M-file displays the FDTD-computed Ez and Hy fields at every
% time step, and records those frames in a movie matrix, M, which is
% played at the end of the simulation using the "movie" command.
%
%***********************************************************************

clear

%***********************************************************************
% Fundamental constants
%***********************************************************************

cc=2.99792458e8; %speed of light in free space
muz=4.0*pi*1.0e-7; %permeability of free space
epsz=1.0/(cc*cc*muz); %permittivity of free space

freq=1.0e+9; %frequency of source excitation
lambda=cc/freq; %wavelength of source excitation
omega=2.0*pi*freq;

%***********************************************************************
% Grid parameters
%***********************************************************************

ie=200; %number of grid cells in x-direction

ib=ie+1;

dx=lambda/20.0; %space increment of 1-D lattice
dt=dx/cc; %time step
omegadt=omega*dt;

nmax=round(12.0e-9/dt); %total number of time steps

%***********************************************************************
% Material parameters
%***********************************************************************

eps=1.0;
sig=5.0e-3;

%***********************************************************************
% Updating coefficients for space region with nonpermeable media
%***********************************************************************

scfact=dt/muz/dx;

ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps));
cb=scfact*(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps));

%***********************************************************************
% Field arrays
%***********************************************************************

ez(1:ib)=0.0;
hy(1:ie)=0.0;

%***********************************************************************
% Movie initialization
%***********************************************************************

x=linspace(dx,ie*dx,ie);

subplot(2,1,1),plot(x,ez(1:ie)/scfact,'r'),axis([0 3 -1 1]);
ylabel('EZ');

subplot(2,1,2),plot(x,hy,'b'),axis([0 3 -3.0e-3 3.0e-3]);
xlabel('x (meters)');ylabel('HY');

rect=get(gcf,'Position');
rect(1:2)=[0 0];

M=moviein(nmax/2,gcf,rect);

%***********************************************************************
% BEGIN TIME-STEPPING LOOP
%***********************************************************************

for n=1:nmax

%***********************************************************************
% Update electric fields
%***********************************************************************

ez(1)=scfact*sin(omegadt*n);

rbc=ez(ie);
ez(2:ie)=ca*ez(2:ie)+cb*(hy(2:ie)-hy(1:ie-1));
ez(ib)=rbc;

%***********************************************************************
% Update magnetic fields
%***********************************************************************

hy(1:ie)=hy(1:ie)+ez(2:ib)-ez(1:ie);

%***********************************************************************
% Visualize fields
%***********************************************************************
7A6sSfPUy

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

共1条评分

cem-uestc

rf币 +5

你发了优秀贴子+2~3

2008-06-17

离线cem-uestc

UID ：9061

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

1楼发表于: 2008-06-17 21:41:49

Matlab代码值得学习，制作仿真动画的好参考matlab代码，其价值超过FDTD代码，奖励

欢迎光临
http://www.mwtee.com/home.php?mod=space&uid=13535

发帖回复