声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1244|回复: 10

请高手们帮忙指点迷津,小弟在这里先谢谢~!

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

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

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

x
<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&lt;m*T_p;<br>                pulse =g_on*sin(omega_0*t);<br>            else <br>                pulse = 0;<br>            end    <br>        case 2<br>            if t&lt;=(m+n)*T_p;<br>                 pulse=sin(omega_0*t);<br>            else <br>                pulse = 0;<br>            end<br>        case 3<br>            if t&lt;=(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&gt;(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编辑过]

回复
分享到:

使用道具 举报

发表于 2006-5-18 11:12 | 显示全部楼层

回复:(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>
发表于 2006-5-18 11:12 | 显示全部楼层

回复:(young20123)请Happy老师帮忙指点迷津,谢谢~!...

其他错误自己在检查检查吧
发表于 2006-5-18 11:13 | 显示全部楼层

回复:(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 说的很对啊~~~
 楼主| 发表于 2006-5-18 15:27 | 显示全部楼层
谢谢,谢谢,小弟知道了
 楼主| 发表于 2006-5-18 15:30 | 显示全部楼层
运行了,就说是<BR>Error: Missing operator, comma, or semicolon.<BR><BR>所以求大家帮忙找错误了
发表于 2006-5-18 15:36 | 显示全部楼层
缺了括号或标点符号什么的,用一用字典也行啊。一般有错的地方它都会告诉你再第几行第几列的!!
[此贴子已经被作者于2006-5-18 15:39:39编辑过]

发表于 2006-5-18 15:38 | 显示全部楼层

回复:(ericlin)缺了括号或标点符号什么的,用一用字...

<DIV class=quote><B>以下是引用<I>ericlin</I>在2006-5-18 15:36:49的发言:</B><BR>缺了括号或标点符号什么的,用一用字典也行啊。</DIV>
<br>呵呵[em17]
 楼主| 发表于 2006-5-18 15:42 | 显示全部楼层
弱弱的问一下,字典怎么用啊,小弟刚刚接触matlab
发表于 2006-5-18 15:50 | 显示全部楼层
<P>打开金山词霸,开启取词功能,把鼠标放到你不认识的英文字上。<br>可能你少了逗号或分号吧。<br><br>你初学MATLAB就编了这么大一串已经很厉害了,继续努力啊。</P>
[此贴子已经被作者于2006-5-18 15:56:27编辑过]

 楼主| 发表于 2006-5-18 15:52 | 显示全部楼层
谢谢大家了
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-11-25 19:53 , Processed in 0.067208 second(s), 17 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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