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

请教校长在matlab上的FDTD程序

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

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

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

回复:(young20123)请教校长在matlab上的FDTD程序

<P>%***********************************************************************<BR>%   1-D FDTD code with simple radiation boundary conditions<BR>%***********************************************************************<BR>%<BR>%   Program author: Susan C. Hagness<BR>%                     Department of Electrical and Computer Engineering<BR>%                     University of Wisconsin-Madison<BR>%                     1415 Engineering Drive<BR>%                     Madison, WI 53706-1691<BR>%                     608-265-5739<BR>%                     <a href="mailthagness@engr.wisc.edu" target="_blank" >hagness@engr.wisc.edu</A><BR>%<BR>%   Date of this version:February 2000<BR>%<BR>%   This MATLAB M-file implements the finite-difference time-domain<BR>%   solution of Maxwell's curl equations over a one-dimensional space<BR>%   lattice comprised of uniform grid cells.<BR>%<BR>%   To illustrate the algorithm, a sinusoidal wave (1GHz) propagating <BR>%   in a nonpermeable lossy medium (epsr=1.0, sigma=5.0e-3 S/m) is <BR>%   modeled.The simplified finite difference system for nonpermeable<BR>%   media (discussed in Section 3.6.6 of the text) is implemented.<BR>%<BR>%   The grid resolution (dx = 1.5 cm) is chosen to provide 20<BR>%   samples per wavelength.The Courant factor S=c*dt/dx is set to<BR>%   the stability limit: S=1.In 1-D, this is the "magic time step."<BR>%<BR>%   The computational domain is truncated using the simplest radiation<BR>%   boundary condition for wave propagation in free space: <BR>%<BR>%                      Ez(imax,n+1) = Ez(imax-1,n)<BR>%<BR>%   To execute this M-file, type "fdtd1D" at the MATLAB prompt.<BR>%   This M-file displays the FDTD-computed Ez and Hy fields at every <BR>%   time step, and records those frames in a movie matrix, M, which is<BR>%   played at the end of the simulation using the "movie" command.<BR>%<BR>%***********************************************************************</P>
<P>clear</P>
<P>%***********************************************************************<BR>%   Fundamental constants<BR>%***********************************************************************</P>
<P>cc=2.99792458e8;            %speed of light in free space<BR>muz=4.0*pi*1.0e-7;          %permeability of free space<BR>epsz=1.0/(cc*cc*muz);       %permittivity of free space</P>
<P>freq=1.0e+9;                %frequency of source excitation<BR>lambda=cc/freq;             %wavelength of source excitation<BR>omega=2.0*pi*freq;</P>
<P>%***********************************************************************<BR>%   Grid parameters<BR>%***********************************************************************</P>
<P>ie=200;                     %number of grid cells in x-direction</P>
<P>ib=ie+1;</P>
<P>dx=lambda/20.0;             %space increment of 1-D lattice<BR>dt=dx/cc;                   %time step<BR>omegadt=omega*dt;</P>
<P>nmax=round(12.0e-9/dt);   %total number of time steps</P>
<P>%***********************************************************************<BR>%   Material parameters<BR>%***********************************************************************</P>
<P>eps=1.0;<BR>sig=5.0e-3;</P>
<P>%***********************************************************************<BR>%   Updating coefficients for space region with nonpermeable media<BR>%***********************************************************************</P>
<P>scfact=dt/muz/dx;</P>
<P>ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps));<BR>cb=scfact*(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps));</P>
<P>%***********************************************************************<BR>%   Field arrays<BR>%***********************************************************************</P>
<P>ez(1:ib)=0.0;<BR>hy(1:ie)=0.0;</P>
<P>%***********************************************************************<BR>%   Movie initialization<BR>%***********************************************************************</P>
<P>x=linspace(dx,ie*dx,ie);</P>
<P>subplot(2,1,1),plot(x,ez(1:ie)/scfact,'r'),axis();<BR>ylabel('EZ');</P>
<P>subplot(2,1,2),plot(x,hy,'b'),axis();<BR>xlabel('x (meters)');ylabel('HY');</P>
<P>rect=get(gcf,'Position');<BR>rect(1:2)=;</P>
<P>M=moviein(nmax/2,gcf,rect);</P>
<P>%***********************************************************************<BR>%   BEGIN TIME-STEPPING LOOP<BR>%***********************************************************************</P>
<P>for n=1:nmax</P>
<P>%***********************************************************************<BR>%   Update electric fields<BR>%***********************************************************************</P>
<P>ez(1)=scfact*sin(omegadt*n);</P>
<P>rbc=ez(ie);<BR>ez(2:ie)=ca*ez(2:ie)+cb*(hy(2:ie)-hy(1:ie-1));<BR>ez(ib)=rbc;</P>
<P>%***********************************************************************<BR>%   Update magnetic fields<BR>%***********************************************************************</P>
<P>hy(1:ie)=hy(1:ie)+ez(2:ib)-ez(1:ie);</P>
<P>%***********************************************************************<BR>%   Visualize fields<BR>%***********************************************************************</P>
<P>if mod(n,2)==0;</P>
<P>rtime=num2str(round(n*dt/1.0e-9));</P>
<P>subplot(2,1,1),plot(x,ez(1:ie)/scfact,'r'),axis();<BR>title(['time = ',rtime,' ns']);<BR>ylabel('EZ');</P>
<P>subplot(2,1,2),plot(x,hy,'b'),axis();<BR>title(['time = ',rtime,' ns']);<BR>xlabel('x (meters)');ylabel('HY');</P>
<P>M(:,n/2)=getframe(gcf,rect);</P>
<P>end</P>
<P>%***********************************************************************<BR>%   END TIME-STEPPING LOOP<BR>%***********************************************************************</P>
<P>end</P>
<P>movie(gcf,M,0,10,rect);<BR></P>

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

三维FDTD的Matlab编程

请问校长有没有关于三维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',),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',),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)=;

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',),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',),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编程原代码?
页: [1]
查看完整版本: 请教校长在matlab上的FDTD程序