请教校长在matlab上的FDTD程序

young20123 · 发表于 2006-3-27 09:11

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

请教校长有没有在matlab上的FDTD程序关于一维时域有限差分法（分层介质，反射，透过系数-FFT），先谢谢校长了

[ 本帖最后由 cdwxg 于 2006-8-22 21:01 编辑 ]

风花雪月 · 发表于 2006-4-1 06:57

%*********************************************************************** % 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 % <a href="mailthagness@engr.wisc.edu" target="_blank" >hagness@engr.wisc.edu</A> % % 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 %***********************************************************************
if mod(n,2)==0;
rtime=num2str(round(n*dt/1.0e-9));
subplot(2,1,1),plot(x,ez(1:ie)/scfact,'r'),axis([0 3 -1 1]); title(['time = ',rtime,' ns']); ylabel('EZ');
subplot(2,1,2),plot(x,hy,'b'),axis([0 3 -3.0e-3 3.0e-3]); title(['time = ',rtime,' ns']); xlabel('x (meters)');ylabel('HY');
M(:,n/2)=getframe(gcf,rect);
end
%*********************************************************************** % END TIME-STEPPING LOOP %***********************************************************************
end
movie(gcf,M,0,10,rect);

sanchabao · 发表于 2006-5-20 15:30

请问校长有没有关于三维FDTDMatlab的程序,谢谢校长了

pheigenbau · 发表于 2006-8-11 02:16

慢慢看

VibInfo · 发表于 2006-8-12 07:50

3维的FDTD程序

%***********************************************************************
% 3-D FDTD code with PEC boundaries
%***********************************************************************
%
% 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 three-dimensional
% Cartesian space lattice comprised of uniform cubic grid cells.
%
% To illustrate the algorithm, an air-filled rectangular cavity
% resonator is modeled. The length, width, and height of the
% cavity are 10.0 cm (x-direction), 4.8 cm (y-direction), and
% 2.0 cm (z-direction), respectively.
%
% The computational domain is truncated using PEC boundary
% conditions:
% ex(i,j,k)=0 on the j=1, j=jb, k=1, and k=kb planes
% ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes
% ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes
% These PEC boundaries form the outer lossless walls of the cavity.
%
% The cavity is excited by an additive current source oriented
% along the z-direction. The source waveform is a differentiated
% Gaussian pulse given by
% J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2),
% where tau=50 ps. The FWHM spectral bandwidth of this zero-dc-
% content pulse is approximately 7 GHz. The grid resolution
% (dx = 2 mm) was chosen to provide at least 10 samples per
% wavelength up through 15 GHz.
%
% To execute this M-file, type "fdtd3D" at the MATLAB prompt.
% This M-file displays the FDTD-computed Ez fields at every other
% 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
%***********************************************************************
% Grid parameters
%***********************************************************************
ie=50; %number of grid cells in x-direction
je=24; %number of grid cells in y-direction
ke=10; %number of grid cells in z-direction
ib=ie+1;
jb=je+1;
kb=ke+1;
is=26; %location of z-directed current source
js=13; %location of z-directed current source
kobs=5;
dx=0.002; %space increment of cubic lattice
dt=dx/(2.0*cc); %time step
nmax=500; %total number of time steps
%***********************************************************************
% Differentiated Gaussian pulse excitation
%***********************************************************************
rtau=50.0e-12;
tau=rtau/dt;
ndelay=3*tau;
srcconst=-dt*3.0e+11;
%***********************************************************************
% Material parameters
%***********************************************************************
eps=1.0;
sig=0.0;
%***********************************************************************
% Updating coefficients
%***********************************************************************
ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps));
cb=(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps));
da=1.0;
db=dt/muz/dx;
%***********************************************************************
% Field arrays
%***********************************************************************
ex=zeros(ie,jb,kb);
ey=zeros(ib,je,kb);
ez=zeros(ib,jb,ke);
hx=zeros(ib,je,ke);
hy=zeros(ie,jb,ke);
hz=zeros(ie,je,kb);
%***********************************************************************
% Movie initialization
%***********************************************************************
tview(:,:)=ez(:,:,kobs);
sview(:,:)=ez(:,js,:);
subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview');
shading flat;
caxis([-1.0 1.0]);
colorbar;
axis image;
title(['Ez(i,j,k=5), time step = 0']);
xlabel('i coordinate');
ylabel('j coordinate');
subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview');
shading flat;
caxis([-1.0 1.0]);
colorbar;
axis image;
title(['Ez(i,j=13,k), time step = 0']);
xlabel('i coordinate');
ylabel('k coordinate');
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
%***********************************************************************
ex(1:ie,2:je,2:ke)=ca*ex(1:ie,2:je,2:ke)+...
cb*(hz(1:ie,2:je,2:ke)-hz(1:ie,1:je-1,2:ke)+...
hy(1:ie,2:je,1:ke-1)-hy(1:ie,2:je,2:ke));
ey(2:ie,1:je,2:ke)=ca*ey(2:ie,1:je,2:ke)+...
cb*(hx(2:ie,1:je,2:ke)-hx(2:ie,1:je,1:ke-1)+...
hz(1:ie-1,1:je,2:ke)-hz(2:ie,1:je,2:ke));
ez(2:ie,2:je,1:ke)=ca*ez(2:ie,2:je,1:ke)+...
cb*(hx(2:ie,1:je-1,1:ke)-hx(2:ie,2:je,1:ke)+...
hy(2:ie,2:je,1:ke)-hy(1:ie-1,2:je,1:ke));
ez(is,js,1:ke)=ez(is,js,1:ke)+...
srcconst*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));
%***********************************************************************
% Update magnetic fields
%***********************************************************************
hx(2:ie,1:je,1:ke)=hx(2:ie,1:je,1:ke)+...
db*(ey(2:ie,1:je,2:kb)-ey(2:ie,1:je,1:ke)+...
ez(2:ie,1:je,1:ke)-ez(2:ie,2:jb,1:ke));
hy(1:ie,2:je,1:ke)=hy(1:ie,2:je,1:ke)+...
db*(ex(1:ie,2:je,1:ke)-ex(1:ie,2:je,2:kb)+...
ez(2:ib,2:je,1:ke)-ez(1:ie,2:je,1:ke));
hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)+...
db*(ex(1:ie,2:jb,2:ke)-ex(1:ie,1:je,2:ke)+...
ey(1:ie,1:je,2:ke)-ey(2:ib,1:je,2:ke));
%***********************************************************************
% Visualize fields
%***********************************************************************
if mod(n,2)==0;
timestep=int2str(n);
tview(:,:)=ez(:,:,kobs);
sview(:,:)=ez(:,js,:);
subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview');
shading flat;
caxis([-1.0 1.0]);
colorbar;
axis image;
title(['Ez(i,j,k=5), time step = ',timestep]);
xlabel('i coordinate');
ylabel('j coordinate');
subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview');
shading flat;
caxis([-1.0 1.0]);
colorbar;
axis image;
title(['Ez(i,j=13,k), time step = ',timestep]);
xlabel('i coordinate');
ylabel('k coordinate');
nn=n/2;
M(:,nn)=getframe(gcf,rect);
end;
%***********************************************************************
% END TIME-STEPPING LOOP
%***********************************************************************
end
movie(gcf,M,0,10,rect);

复制代码

bingbing1341 · 发表于 2006-8-22 15:30

请问校长一个关于仿真天线的问题,: 怎样才能在fdtd程序中加入要仿真的器件呢？有相关的仿真程序校长能给新手法一份么？非常感谢！！！我的邮箱：bingbing1341@163.com

mermaid · 发表于 2008-3-3 09:23

请问校长有没有二维的光子晶体电磁场传输特性的FTDT的matlab程序，非常非常感谢！

nanomod · 发表于 2008-12-4 21:18

正好需要可惜还是看不懂！

科技在线 · 发表于 2008-12-4 21:21

相当的复杂，存下来慢慢消化

penglilong008 · 发表于 2011-5-8 12:09

请问;一维光子晶体传输矩阵模拟禁带Matlab编程原代码?

账号		自动登录	找回密码
密码			我要加入

[综合讨论] 请教校长在matlab上的FDTD程序

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

回复：（young20123）请教校长在matlab上的FDTD程序

三维FDTD的Matlab编程