[求助]新手请教taflove书中的一维程序有些疑问

%*********************************************************************** W%4=x>J-  
%     1-D FDTD code with simple radiation boundary conditions J91[w?,  
%*********************************************************************** < Hkq  
% E7t;p)x  
%     Program author: Susan C. Hagness I]^>>>p$  
%                     Department of Electrical and Computer Engineering kk%32(By  
%                     University of Wisconsin-Madison <XIIT-b[  
%                     1415 Engineering Drive #obRr#8  
%                     Madison, WI 53706-1691 @(Y!$><Is  
%                     608-265-5739 =}vT>b
  
%                     [url=mailto:hagness@engr.wisc.edu]hagness@engr.wisc.edu[/url] )\|+G5#`  
% MC* Hl`C  
%     Date of this version:  February 2000 !bP%\)5  
% &'SD1m1P  
%     This MATLAB M-file implements the finite-difference time-domain z@yTkH_  
%     solution of Maxwell's curl equations over a one-dimensional space -~mgct5  
%     lattice comprised of uniform grid cells. 0?<#!  
% tO&ffZP8$  
%     To illustrate the algorithm, a sinusoidal wave (1GHz) propagating *!%lBt{2  
%     in a nonpermeable lossy medium (epsr=1.0, sigma=5.0e-3 S/m) is 0V4B Q:v  
%     modeled.  The simplified finite difference system for nonpermeable IdQ./@?  
%     media (discussed in Section 3.6.6 of the text) is implemented.97081162 &^r>Q`u  
% b~J)LXj]w  
%     The grid resolution (dx = 1.5 cm) is chosen to provide 20 ;$QC_l''b  
%     samples per wavelength.  The Courant factor S=c*dt/dx is set to N ~{N Nf Y  
%     the stability limit: S=1.  In 1-D, this is the "magic time step." J0?kEr  
% f#=c=e-A  
%     The computational domain is truncated using the simplest radiation  .qgUD  
%     boundary condition for wave propagation in free space: 1Y|a:){G  
% ")T\_ME  
%                      Ez(imax,n+1) = Ez(imax-1,n) LWyr  
% |5bLV^mv]i  
%     To execute this M-file, type "fdtd1D" at the MATLAB prompt. nN\H'{Wzd  
%     This M-file displays the FDTD-computed Ez and Hy fields at every 'v"=  
%     time step, and records those frames in a movie matrix, M, which is 5J10S  
%     played at the end of the simulation using the "movie" command. /]_t->  
% gv''A"  
%*********************************************************************** ;f=m+QXU  
clear 8NWo)y49H  
%*********************************************************************** Ix5&B6L8  
%     Fundamental constants W6&vyOc  
%*********************************************************************** EU$.{C_O(  
cc=2.99792458e8;            %speed of light in free space HeOdCr-PN  
muz=4.0*pi*1.0e-7;          %permeability of free space s V_(9@b  
epsz=1.0/(cc*cc*muz);       %permittivity of free space
WwDM^}e  
freq=1.0e+9;                %frequency of source excitation yOq@w!xz  
lambda=cc/freq;             %wavelength of source excitation .\n` 4A1z  
omega=2.0*pi*freq; 4f([EV[6dK  
%*********************************************************************** Fl-\{vOn  
%     Grid parameters }d<R 5  
%*********************************************************************** {'5"i?>s0>  
ie=200;                     %number of grid cells in x-direction "9wD|wsz  
ib=ie+1; 2;8m0+tl  
dx=lambda/20.0;             %space increment of 1-D lattice {_QdB;VwH  
dt=dx/cc;                   %time step 7l D-|yx  
omegadt=omega*dt; !y= R)k  
nmax=round(12.0e-9/dt);     %total number of time steps zaqX};b  
%*********************************************************************** ;V xRaj?  
%     Material parameters <s9?9^!!V^  
%*********************************************************************** i"WYcF|  
eps=1.0; zxbfh/=  
sig=5.0e-3; wI$a1H  
%*********************************************************************** X2z<cJG|d@  
%     Updating coefficients for space region with nonpermeable media xT%`"eM}  
%*********************************************************************** ']r8q 
%  
scfact=dt/muz/dx; DN*5q9.  
ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps)); 9 r!zYZ`)  
cb=scfact*(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps));        这个式子的我认为他表达的不对分子应该是:(dt/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps)) mQj=-\p  
%*********************************************************************** uu9M}]mDl  
%     Field arrays  Y{p$%  
%*********************************************************************** yD7BZI xW  
ez(1:ib)=0.0; FACw;/rW  
hy(1:ie)=0.0; ;2p+i/sVj  
%*********************************************************************** X\AH^I6S  
%     Movie initialization nlwqSXw  
%*********************************************************************** ^7-zwl(>?N  
x=linspace(dx,ie*dx,ie); =dmr
,WE  
subplot(2,1,1),plot(x,ez(1:ie)/scfact,'r'),axis([0 3 -1 1]); ;mkkaW,D*  
ylabel('EZ'); Y;"k5+ q  
subplot(2,1,2),plot(x,hy,'b'),axis([0 3 -3.0e-3 3.0e-3]); S#7YJ7 K"N  
xlabel('x (meters)');ylabel('HY'); *l+#<5x  
rect=get(gcf,'Position'); 这个式子是往下的3行是啥意思 j/FLEsU!R  
rect(1:2)=[0 0]; [9 W@<p  
M=moviein(nmax/2,gcf,rect); F5h/>  
%*********************************************************************** c0qp-=^&.  
%     BEGIN TIME-STEPPING LOOP D 2X_Yv
 
%*********************************************************************** rtV`Q[E  
for n=1:nmax fw%`[(hK  
%*********************************************************************** 7>FXsUt_  
%     Update electric fields  =<HDek  
%*********************************************************************** ,NSf  
ez(1)=scfact*sin(omegadt*n); 9iA rBL"  
rbc=ez(ie); i
%hCV o  
ez(2:ie)=ca*ez(2:ie)+cb*(hy(2:ie)-hy(1:ie-1)); ZJYn[\]  
ez(ib)=rbc; Lo%n{*if  
%*********************************************************************** Cn{Hk)6  
%     Update magnetic fields l8\UO<^fY  
%*********************************************************************** gcJ!_KZK  
hy(1:ie)=hy(1:ie)+ez(2:ib)-ez(1:ie);  这个式子?不是公式不会推而是这个式子的表达跟葛老师书上的公式表达不一致,希望能给解答 Zxa.x?:?n  
%*********************************************************************** N
ep4J;  
%     Visualize fields LeKovt%  
%*********************************************************************** >f(?Mxh2  
if mod(n,2)==0; szI7I$Qb  
rtime=num2str(round(n*dt/1.0e-9)); b/wpk~qi  
subplot(2,1,1),plot( .. x:|Y)Dn\  
Zf'*pp T&q  

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

rect=get(gcf,'Position'); 这个式子是往下的3行是啥意思 aSL`yuXu  
rect(1:2)=[0 0]; HU3:6R&  
M=moviein(nmax/2,gcf,rect); +7Ws`qhEe  
M(:,n/2)=getframe(gcf,rect);这个式子的意思? s#2t\}/  
_;lw,;ftA  
这是matlab动画显示的，多少帧显示。随便找本matlab的书上都有的。