声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 4479|回复: 9

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

[复制链接]
发表于 2006-3-27 09:11 | 显示全部楼层 |阅读模式

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

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

x
请教校长有没有在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([0 3 -1 1]);<BR>ylabel('EZ');</P>
<P>subplot(2,1,2),plot(x,hy,'b'),axis([0 3 -3.0e-3 3.0e-3]);<BR>xlabel('x (meters)');ylabel('HY');</P>
<P>rect=get(gcf,'Position');<BR>rect(1:2)=[0 0];</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([0 3 -1 1]);<BR>title(['time = ',rtime,' ns']);<BR>ylabel('EZ');</P>
<P>subplot(2,1,2),plot(x,hy,'b'),axis([0 3 -3.0e-3 3.0e-3]);<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>
发表于 2006-5-20 15:30 | 显示全部楼层

三维FDTD的Matlab编程

请问校长有没有关于三维FDTDMatlab的程序,谢谢校长了
发表于 2006-8-11 02:16 | 显示全部楼层
慢慢看
发表于 2006-8-12 07:50 | 显示全部楼层
3维的FDTD程序

  1. %***********************************************************************
  2. %     3-D FDTD code with PEC boundaries
  3. %***********************************************************************
  4. %
  5. %     Program author: Susan C. Hagness
  6. %                     Department of Electrical and Computer Engineering
  7. %                     University of Wisconsin-Madison
  8. %                     1415 Engineering Drive
  9. %                     Madison, WI 53706-1691
  10. %                     608-265-5739
  11. %                     hagness@engr.wisc.edu
  12. %
  13. %     Date of this version:  February 2000
  14. %
  15. %     This MATLAB M-file implements the finite-difference time-domain
  16. %     solution of Maxwell's curl equations over a three-dimensional
  17. %     Cartesian space lattice comprised of uniform cubic grid cells.
  18. %     
  19. %     To illustrate the algorithm, an air-filled rectangular cavity
  20. %     resonator is modeled.  The length, width, and height of the
  21. %     cavity are 10.0 cm (x-direction), 4.8 cm (y-direction), and
  22. %     2.0 cm (z-direction), respectively.
  23. %
  24. %     The computational domain is truncated using PEC boundary
  25. %     conditions:
  26. %          ex(i,j,k)=0 on the j=1, j=jb, k=1, and k=kb planes
  27. %          ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes
  28. %          ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes
  29. %     These PEC boundaries form the outer lossless walls of the cavity.
  30. %
  31. %     The cavity is excited by an additive current source oriented
  32. %     along the z-direction.  The source waveform is a differentiated
  33. %     Gaussian pulse given by
  34. %          J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2),
  35. %     where tau=50 ps.  The FWHM spectral bandwidth of this zero-dc-
  36. %     content pulse is approximately 7 GHz. The grid resolution
  37. %     (dx = 2 mm) was chosen to provide at least 10 samples per
  38. %     wavelength up through 15 GHz.
  39. %
  40. %     To execute this M-file, type "fdtd3D" at the MATLAB prompt.
  41. %     This M-file displays the FDTD-computed Ez fields at every other
  42. %     time step, and records those frames in a movie matrix, M, which
  43. %     is played at the end of the simulation using the "movie" command.
  44. %
  45. %***********************************************************************

  46. clear

  47. %***********************************************************************
  48. %     Fundamental constants
  49. %***********************************************************************

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

  53. %***********************************************************************
  54. %     Grid parameters
  55. %***********************************************************************

  56. ie=50;       %number of grid cells in x-direction
  57. je=24;       %number of grid cells in y-direction
  58. ke=10;       %number of grid cells in z-direction

  59. ib=ie+1;     
  60. jb=je+1;   
  61. kb=ke+1;   

  62. is=26;       %location of z-directed current source
  63. js=13;       %location of z-directed current source

  64. kobs=5;

  65. dx=0.002;          %space increment of cubic lattice
  66. dt=dx/(2.0*cc);    %time step

  67. nmax=500;          %total number of time steps

  68. %***********************************************************************
  69. %     Differentiated Gaussian pulse excitation
  70. %***********************************************************************

  71. rtau=50.0e-12;
  72. tau=rtau/dt;
  73. ndelay=3*tau;
  74. srcconst=-dt*3.0e+11;

  75. %***********************************************************************
  76. %     Material parameters
  77. %***********************************************************************

  78. eps=1.0;
  79. sig=0.0;        

  80. %***********************************************************************
  81. %     Updating coefficients
  82. %***********************************************************************

  83. ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps));
  84. cb=(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps));
  85. da=1.0;
  86. db=dt/muz/dx;

  87. %***********************************************************************
  88. %     Field arrays
  89. %***********************************************************************

  90. ex=zeros(ie,jb,kb);
  91. ey=zeros(ib,je,kb);
  92. ez=zeros(ib,jb,ke);
  93. hx=zeros(ib,je,ke);
  94. hy=zeros(ie,jb,ke);
  95. hz=zeros(ie,je,kb);

  96. %***********************************************************************
  97. %     Movie initialization
  98. %***********************************************************************

  99. tview(:,:)=ez(:,:,kobs);
  100. sview(:,:)=ez(:,js,:);

  101. subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview');
  102. shading flat;
  103. caxis([-1.0 1.0]);
  104. colorbar;
  105. axis image;
  106. title(['Ez(i,j,k=5), time step = 0']);
  107. xlabel('i coordinate');
  108. ylabel('j coordinate');

  109. subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview');
  110. shading flat;
  111. caxis([-1.0 1.0]);
  112. colorbar;
  113. axis image;
  114. title(['Ez(i,j=13,k), time step = 0']);
  115. xlabel('i coordinate');
  116. ylabel('k coordinate');

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

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

  120. %***********************************************************************
  121. %     BEGIN TIME-STEPPING LOOP
  122. %***********************************************************************

  123. for n=1:nmax
  124.    
  125. %***********************************************************************
  126. %     Update electric fields
  127. %***********************************************************************

  128. ex(1:ie,2:je,2:ke)=ca*ex(1:ie,2:je,2:ke)+...
  129.                    cb*(hz(1:ie,2:je,2:ke)-hz(1:ie,1:je-1,2:ke)+...
  130.                        hy(1:ie,2:je,1:ke-1)-hy(1:ie,2:je,2:ke));

  131. ey(2:ie,1:je,2:ke)=ca*ey(2:ie,1:je,2:ke)+...
  132.                    cb*(hx(2:ie,1:je,2:ke)-hx(2:ie,1:je,1:ke-1)+...
  133.                        hz(1:ie-1,1:je,2:ke)-hz(2:ie,1:je,2:ke));
  134.                     
  135. ez(2:ie,2:je,1:ke)=ca*ez(2:ie,2:je,1:ke)+...
  136.                    cb*(hx(2:ie,1:je-1,1:ke)-hx(2:ie,2:je,1:ke)+...
  137.                        hy(2:ie,2:je,1:ke)-hy(1:ie-1,2:je,1:ke));
  138.                     
  139. ez(is,js,1:ke)=ez(is,js,1:ke)+...
  140.                srcconst*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));

  141. %***********************************************************************
  142. %     Update magnetic fields
  143. %***********************************************************************

  144. hx(2:ie,1:je,1:ke)=hx(2:ie,1:je,1:ke)+...
  145.                    db*(ey(2:ie,1:je,2:kb)-ey(2:ie,1:je,1:ke)+...
  146.                        ez(2:ie,1:je,1:ke)-ez(2:ie,2:jb,1:ke));
  147.                
  148. hy(1:ie,2:je,1:ke)=hy(1:ie,2:je,1:ke)+...
  149.                    db*(ex(1:ie,2:je,1:ke)-ex(1:ie,2:je,2:kb)+...
  150.                        ez(2:ib,2:je,1:ke)-ez(1:ie,2:je,1:ke));
  151.                
  152. hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)+...
  153.                    db*(ex(1:ie,2:jb,2:ke)-ex(1:ie,1:je,2:ke)+...
  154.                        ey(1:ie,1:je,2:ke)-ey(2:ib,1:je,2:ke));
  155.                     
  156. %***********************************************************************
  157. %     Visualize fields
  158. %***********************************************************************

  159. if mod(n,2)==0;

  160. timestep=int2str(n);
  161. tview(:,:)=ez(:,:,kobs);
  162. sview(:,:)=ez(:,js,:);

  163. subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview');
  164. shading flat;
  165. caxis([-1.0 1.0]);
  166. colorbar;
  167. axis image;
  168. title(['Ez(i,j,k=5), time step = ',timestep]);
  169. xlabel('i coordinate');
  170. ylabel('j coordinate');

  171. subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview');
  172. shading flat;
  173. caxis([-1.0 1.0]);
  174. colorbar;
  175. axis image;
  176. title(['Ez(i,j=13,k), time step = ',timestep]);
  177. xlabel('i coordinate');
  178. ylabel('k coordinate');

  179. nn=n/2;
  180. M(:,nn)=getframe(gcf,rect);

  181. end;

  182. %***********************************************************************
  183. %     END TIME-STEPPING LOOP
  184. %***********************************************************************

  185. end

  186. movie(gcf,M,0,10,rect);
复制代码
发表于 2006-8-22 15:30 | 显示全部楼层
请问校长一个关于仿真天线的问题,: 怎样才能在fdtd程序中加入要仿真的器件呢? 有相关的仿真程序校长能给新手法一份么?非常感谢!!!我的邮箱:bingbing1341@163.com
发表于 2008-3-3 09:23 | 显示全部楼层
请问校长有没有二维的光子晶体电磁场传输特性的FTDT的matlab程序,非常非常感谢!
发表于 2008-12-4 21:18 | 显示全部楼层
正好需要可惜还是看不懂!
发表于 2008-12-4 21:21 | 显示全部楼层
相当的复杂,存下来慢慢消化
发表于 2011-5-8 12:09 | 显示全部楼层
请问;一维光子晶体传输矩阵模拟禁带Matlab编程原代码?
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-11 07:16 , Processed in 0.076326 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表