请高手们帮忙指点迷津,小弟在这里先谢谢~!
<P>这是我编的一个matlab程序,是用一维FDTD方法对正弦电磁波穿过DNG媒质的平行板时的传播过程模拟,出了问题,但是不知道在哪,希望您能帮指点一下错误.<br><br>还有,能不能提供一个PML边界吸收条件的matlab程序,万分谢谢了~!<br><br>%constants<br>c_0 = 3E8; % Speed of light in free space<br>Nz = 5000; % Number of cells in z-direction<br>Nt = 2000; % Number of time steps<br>N = Nz; % Global number of cells</P><P>E = zeros(Nz,1); % E component<br>E2 = zeros(Nz,1);<br>E1 = zeros(Nz,1);<br>Es = zeros(Nz,Nz); % Es is the output for surf function<br>%H = zeros(Nz,1); % H component<br>pulse = 0; % Pulse</P>
<P>%to be set<br>mu_0 = 4.0*pi*1.0e-7; % Permeability of free space<br>eps_0 = 1.0/(c_0*c_0*mu_0);% Permittivty of free space</P>
<P>%c_ref = c_0; % Reference velocity<br>%eps_ref = eps_0;<br>%mu_ref = mu_0;</P>
<P>f_0 = 3*10e10; % Excitation frequency<br>%f_ref = f_0; % Reference frequency </P>
<P>omega_0 = 2.0*pi*f_0; % Excitation circular frequency<br>%omega_ref = omega_0;</P>
<P>Dt=9.5*10e-14;<br>Dz=3.0*10e-5;</P>
<P>r = 1; </P>
<P>T_p=1/f_0;<br>omega_p=10e8;</P>
<P>tao=10e8;</P>
<P>m=5;<br>n=10;</P>
<P>% Source position<br>Nz_Source = 0.12*Nz;<br>N_Source = Nz_Source;</P>
<P>for n1 = 1:Nt-1<br> t=Dt*r*(n1-1);<br> <br> x_on=1-(m*T_p-t)/(m*Tp);<br> x_off=(t-(m+n)*T_p)/(m*T_p);<br> g_on=10*x_on^3-15*x_on^4+6*x_on^5;<br> g_off=1-(10*x_off^3-15*x_off^4+6*x_off^5);<br> <br> %Source function series<br> Source_type = 1;<br> switch Source_type <br> case 1<br> if t<m*T_p;<br> pulse =g_on*sin(omega_0*t);<br> else <br> pulse = 0;<br> end <br> case 2<br> if t<=(m+n)*T_p;<br> pulse=sin(omega_0*t);<br> else <br> pulse = 0;<br> end<br> case 3<br> if t<=(m+n+m)*T_p;<br> pulse=g_off*sin(omega_0*t);<br> else <br> pulse = 0;<br> end <br> case 4<br> if t>(m+n+m)*T_p;<br> pulse = 0;<br> else <br> pulse = 0;<br> end<br> end<br> <br> E1(N_Source) = E1(N_Source) - r*pulse;<br> E1 = E;<br> <br> H1(N_Source) = H1(N_Source) - r*pulse;<br> H1 = H;<br> <br> J=diff(E,t);<br> J1=diff(E1,t);<br> <br> K=diff(H,t);<br> K1=diff(H1,t);<br> <br> i=2:Nz-1<br> E(i)=E1(i)-(Dt/(eps_0*Dz))*(H(i+0.5)-H(i-0.5)+0.5*(J1(i+0.5)+J1(i-0.5))*Dz);<br> H(i+0.5)=H1(i+0.5)-(Dt/(mu_0*Dz))*(E1(i+1)-E1(i)+K1(i+0.5)*Dz);<br> K(i+0.5)=((1-0.5*tao*Dt)/(1+0.5*tao*Dt))*K1(i+0.5)+(mu_0*omega_p^2*Dt/(1+0.5*tao*Dt))*H(i+0.5);<br> J(i+0.5)=((1-0.5*tao*Dt)/(1+0.5*tao*Dt))*J1(i+0.5)+(0.5*eps_0*omega_p^2*Dt/(1+0.5*tao*Dt))*(E(i)+E(i+1));<br> <br> %Boundary<br> E(1) = 0; E(Nz) = 0; % Dirichlet<br> <br> fprintf(fid,'%10.4f%10.4f\n',E(i),t);<br>end<br> plot(E(i),t);<br> hold on<br></P>
[此贴子已经被作者于2006-5-18 15:31:27编辑过]
回复：（young20123）请Happy老师帮忙指点迷津,谢谢~!...
<P>两个个非常明显的错误</P><P>x_on=1-(m*T_p-t)/(m*Tp);<BR>其中Tp应该是T_p吧</P>
<P>H1(N_Source) = H1(N_Source) - r*pulse;<BR>H1没有初始值</P>
<P>这种错误自己看一下错误信息应该就能看出来吧</P>
<P>另外给你一个简易，问问题的时候最好不要指名问谁，否则有些人看到了就算知道也未必回答你</P>
回复：（young20123）请Happy老师帮忙指点迷津,谢谢~!...
其他错误自己在检查检查吧回复：（suffer）回复：（young20123）请Happy老师帮忙...
<DIV class=quote><B>以下是引用<I>suffer</I>在2006-5-18 11:12:32的发言：</B><BR><P>两个个非常明显的错误</P>
<P>x_on=1-(m*T_p-t)/(m*Tp);<BR>其中Tp应该是T_p吧</P>
<P>H1(N_Source) = H1(N_Source) - r*pulse;<BR>H1没有初始值</P>
<P>这种错误自己看一下错误信息应该就能看出来吧</P>
<P>另外给你一个简易，问问题的时候最好不要指名问谁，否则有些人看到了就算知道也未必回答你</P></DIV>
<br>suffer 说的很对啊~~~ 谢谢,谢谢,小弟知道了 运行了,就说是<BR>Error: Missing operator, comma, or semicolon.<BR><BR>所以求大家帮忙找错误了 缺了括号或标点符号什么的，用一用字典也行啊。一般有错的地方它都会告诉你再第几行第几列的！！
[此贴子已经被作者于2006-5-18 15:39:39编辑过]
回复：（ericlin）缺了括号或标点符号什么的，用一用字...
<DIV class=quote><B>以下是引用<I>ericlin</I>在2006-5-18 15:36:49的发言：</B><BR>缺了括号或标点符号什么的，用一用字典也行啊。</DIV><br>呵呵 弱弱的问一下,字典怎么用啊,小弟刚刚接触matlab <P>打开金山词霸，开启取词功能，把鼠标放到你不认识的英文字上。<br>可能你少了逗号或分号吧。<br><br>你初学MATLAB就编了这么大一串已经很厉害了，继续努力啊。</P>
[此贴子已经被作者于2006-5-18 15:56:27编辑过]
谢谢大家了
页:
[1]