743阅读
0回复

[求助]一维程序求解读

离线422204942

UID ：74635

注册：2011-03-28
登录：2011-05-20
发帖：25
等级：仿真新人

0楼发表于: 2011-05-09 13:20:32

%%Try changing the value of "profile".
close all;
clear all;
L = 5; %长度
N = 505; %# 空间划分
Niter = 400; %迭代次数
fs = 300e6; %频率
ds = L/N; %步长
dt = ds/300e6; %”单位时间“
eps0 = 8.854e-12; %真空中的介电常数
mu0 = pi*4e-7; %真空中的磁导率
x = linspace(0,L,N); %x 轴划分
showWKB=0;
ae = ones(N,1)*dt/(ds*eps0);
am = ones(N,1)*dt/(ds*mu0);
as = ones(N,1);
epsr = ones(N,1);
mur= ones(N,1);
sigma = zeros(N,1);
%Try profile = 1,2,3,4,5,6 in sequence.
%You can define epsr(i)介电常数ε=εr*ε0, mur(i)磁导率μ %and sigma(i)导电性
profile =1;
for i=1:N/2
   epsr(i) = 1;
   mur(i) = 1;
   w1 = 0.5;
   w2 = 1.5;
   if (profile==1) %介电窗口
      if (abs(x(i)-L/2)<0.5) epsr(i)=4; end
   end
   if (profile==2)%平稳过渡
   if (abs(x(i)-L/2)<1.5) epsr(i)=1+3*(1+cos(pi*(abs(x(i)-L/2)-w1)/(w2-w1)))/2; end
      if (abs(x(i)-L/2)<0.5) epsr(i)=4; end
   end
   if (profile==3) %介质不连续
      if (x(i)>L/2) epsr(i) = 9; end
   end
if (profile==4)
      if (x(i)>(L/2-0.1443)) epsr(i) = 3; end
      if (x(i)>L/2) epsr(i) = 9; end
   end

end
ae = ae./epsr;
am = am./mur;
ae = ae./(1+dt*(sigma./epsr)/(2*eps0));
as = (1-dt*(sigma./epsr)/(2*eps0))./(1+dt*(sigma./epsr)/(2*eps0));
Hz = zeros(N,1);
Ey = zeros(N,1);
figure(2);
set(gcf,'doublebuffer','on');
plot(Ey);
grid on;
for iter=1:Niter
   Ey(3) = Ey(3)+2*(1-exp(-((iter-1)/50)^2))*sin(2*pi*fs*dt*iter);
   Hz(1) = Hz(2);
   for i=2:N-1
      Hz(i) = Hz(i)-am(i)*(Ey(i+1)-Ey(i));
   end
   Ey(N) = Ey(N-1);
   for i=2:N-1
      Ey(i) = as(i)*Ey(i)-ae(i)*(Hz(i)-Hz(i-1));
   end
   figure(2)
   hold off
   plot(x,Ey,'b');
   axis([3*ds L -2 2]);
   grid on;
   title('E (blue) and 377*H (red)');
   hold on
   plot(x,377*Hz,'r');
   xlabel('x (m)');
   pause(0);
   iter
end
if (showWKB==1)
phase = cumsum((epsr).^0.5)*ds;
beta0 = 2*pi*fs/(300e6);
   theory = sin(2*pi*fs*(Niter+4)*dt-beta0*phase)./(epsr.^0.25);
   ..

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

发帖回复