|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
<P>下面是小弟编的程,出问题,解决不了了,哪位好心的高手帮忙看一下.<BR><BR>function FDTD_debug</P>
<P>%constants<BR>c_0 = 3E8; % Speed of light in free space<BR>Nz = 100; % Number of cells in z-direction<BR>Nt = 200; % Number of time steps<BR>N = Nz; % Global number of cells</P>
<P>E = zeros(Nz,1); % E component<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>H1 = zeros(Nz,1); <BR>J = zeros(Nz,1);<BR>J1 = zeros(Nz,1);<BR>K = zeros(Nz,1);<BR>K1 = zeros(Nz,1);<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 = 10e9; % 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;<BR>omega_p=10e8;</P>
<P>tao=10e8;</P>
<P>lambda_ref = c_ref/f_ref; % Reference wavelength</P>
<P>Dx_ref = lambda_ref/20; % Reference cells width<BR>Dz = Dx_ref;<BR>nDz = Dz/Dx_ref;</P>
<P>Z = Nz*Dz;<BR>r = 1; % Normalized time step<BR>Dt_ref = r*Dx_ref/c_ref; % Reference time step<BR>Dt = Dt_ref;</P>
<P>% Source position<BR>Nz_Source = 0.5*Nz;<BR>N_Source = Nz_Source;</P>
<P>for n = 1:Nt-1<BR> t = Dt_ref*r*(n-1); % Actual time<BR> <BR> %Source function series<BR> Source_type = 1;<BR> switch Source_type <BR> case 1 % modified source function<BR> ncycles = 2;<BR> if t < ncycles*2.0*pi/(omega_0) <BR> pulse = -0.5*( 1.0 - cos(omega_0*t/ncycles) ) * sin(omega_0*t);<BR> else <BR> pulse = 0;<BR> end<BR> case 2 % sigle cos source function<BR> if t < 2.0*pi/(omega_0)<BR> pulse = 8*c_0^2*Dt_ref^2*mu_ref*omega_0*cos(omega_0*t);<BR> else <BR> pulse = 0;<BR> end<BR> case 3 % Gaussian pulse<BR> if t < Dt_ref*r*50<BR> pulse = -40*c_0*(t-t*25/(n-1))*exp(-(t-t*25/(n-1))^2/2/(50/2.3548)^2/(t/(n-1))^2);<BR> else<BR> pulse = 0;<BR> end <BR> end<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> J1(N_Source) = J1(N_Source) - r*pulse;<BR> J1 = J;<BR> <BR> K1(N_Source) = K1(N_Source) - r*pulse;<BR> K1 = K;<BR> <BR> b=2:Nz-1<BR> E(b)=E1(b)-(Dt/(eps_0*Dz))*(H(b+0.5)-H(b-0.5)+0.5*(J1(b+0.5)+J1(b-0.5))*Dz);<BR> H(b+0.5)=H1(b+0.5)-(Dt/(mu_0*Dz))*(E1(b+1)-E1(b)+K1(b+0.5)*Dz);<BR> K(b+0.5)=((1-0.5*tao*Dt)/(1+0.5*tao*Dt))*K1(b+0.5)+(mu_0*omega_p^2*Dt/(1+0.5*tao*Dt))*H(b+0.5);<BR> J(b+0.5)=((1-0.5*tao*Dt)/(1+0.5*tao*Dt))*J1(b+0.5)+(0.5*eps_0*omega_p^2*Dt/(1+0.5*tao*Dt))*(E(b)+E(b+1));<BR> <BR> %Boundary<BR> E(1) = 0; E(Nz) = 0; % Dirichlet<BR> %E(1) = E(2);E(Nz) = E(Nz-1); % Neumann<BR> %E(Nz) = E2(Nz-1) + (r-1)/(r+1)*(E(Nz-1)-E2(Nz)); % Mur ABC @ right side z = Nz<BR> %E(1) = E2(2) + (r-1)/(r+1)*(E(2)-E2(1)); % Mur ABC @ left side z = 0<BR> <BR> %display<BR> %choice***********************************************<BR> display = 3; % choice: '1' = line plot; '2' = imagesc; '3' = surf<BR> switch display<BR> case 1<BR> i = 1:Nz;<BR> plot(i, E(i));<BR> axis([0 Nz -4 4]);<BR> case 2<BR> A = rot90(E);<BR> imagesc(A);<BR> case 3<BR> i = 1:Nz;<BR> for j = 1:Nz<BR> Es(i,j) = E(i);<BR> end<BR> Es = rot90(Es);<BR> j = 1:Nz;<BR> surf(i,j,Es);<BR> axis([0 Nz 0 Nz -4 4]) <BR> otherwise<BR> A = rot90(E);<BR> imagesc(A);<BR> end <BR> pause(0.005); <BR>end<BR><BR>*************************************************************************************************************<BR><FONT color=#ff0033>??? Subscript indices must either be real positive integers or logicals.</FONT></P>
<P><FONT color=#ff0033>Error in ==> E:\matlab6.5\work\FDTD1Dme.m<BR>On line 91 ==> E(b)=E1(b)-(Dt/(eps_0*Dz))*(E(b+0.5)-E(b-0.5)+0.5*(E1(b+0.5)+E1(b-0.5))*Dz);</FONT></P>[em06] |
|