[求助]一维程序求解读

%%Try changing the value of "profile". m*Zq3j  
close all; 4LEWOWF}  
clear all; (E{>L).~  
L = 5; %长度 ($Cy-p  
N = 505; %# 空间划分  n[vwwY  
Niter = 400; %迭代次数 /5Od:n  
fs = 300e6; %频率 TY."?` [FK  
ds = L/N; %步长 |fL|tkGEa  
dt = ds/300e6; %”单位时间“ jGg,)~)Y  
eps0 = 8.854e-12; %真空中的介电常数
F9ys.Bc  
mu0 = pi*4e-7; %真空中的磁导率 <EhOIN7@*D  
x = linspace(0,L,N); %x 轴划分 ~[_u@8l!mN  
showWKB=0; PykVXZ7j;  
ae = ones(N,1)*dt/(ds*eps0); ;6 ?a8t@  
am = ones(N,1)*dt/(ds*mu0); ;PS V3Zh  
as = ones(N,1); B0h|Y.S8%1  
epsr = ones(N,1); R[C+?qux  
mur= ones(N,1); k? <.y
r1  
sigma = zeros(N,1); jR1o<]
?  
%Try profile = 1,2,3,4,5,6 in sequence. FAkrM?0/  
%You can define epsr(i)介电常数ε=εr*ε0, mur(i)磁导率μ  %and sigma(i)导电性 ?/M:  
profile =1;
)8cb @N  
for i=1:N/2 q,i&%  
   epsr(i) = 1; ~.Wlv;  
   mur(i) = 1; VjI=5)+~  
   w1 = 0.5; *~D|M  
   w2 = 1.5; "QKCZ8_C  
   if (profile==1) %介电窗口 \;!}z3Ww  
      if (abs(x(i)-L/2)<0.5) epsr(i)=4; end ; _ziRy  
   end 6'{/Ote  
   if (profile==2)%平稳过渡 GI se|[p  
    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 Fy; sVB  
      if (abs(x(i)-L/2)<0.5) epsr(i)=4; end -w dbH`2Z"  
   end IP LKOT~  
   if (profile==3) %介质不连续 ^n|yfvR  
      if (x(i)>L/2) epsr(i) = 9; end S/itK3  
   end XIGz_g;#'w  
   if (profile==4) V-{3)6I$hG  
      if (x(i)>(L/2-0.1443)) epsr(i) = 3; end ~+A(zlYr~  
      if (x(i)>L/2) epsr(i) = 9; end i,H(6NL.  
   end h
SeXxSb:  
  %ER"Udh  
end 4.=jKj9j  
ae = ae./epsr; uPT2ga]  
am = am./mur; GM~Ek]9C%  
ae = ae./(1+dt*(sigma./epsr)/(2*eps0)); URd0|?t9^L  
as = (1-dt*(sigma./epsr)/(2*eps0))./(1+dt*(sigma./epsr)/(2*eps0)); kxanzsSr9  
Hz = zeros(N,1); @A5'vf|2;.  
Ey = zeros(N,1); t)4><22of  
figure(2); =bBV A0y  
set(gcf,'doublebuffer','on'); 5}3#l/  
plot(Ey); ^xyU*A}D  
grid on; rD\)ndPv  
for iter=1:Niter   _o52#Q4   
   Ey(3) = Ey(3)+2*(1-exp(-((iter-1)/50)^2))*sin(2*pi*fs*dt*iter); x>cl$41!W  
   Hz(1) = Hz(2); +]3kcm7B  
   for i=2:N-1 9-W3}4'e  
      Hz(i) = Hz(i)-am(i)*(Ey(i+1)-Ey(i)); )+ V)]dS@%  
   end o=nF.y  
   Ey(N) = Ey(N-1); K7 J RCLA  
   for i=2:N-1 }T5 E^  
      Ey(i) = as(i)*Ey(i)-ae(i)*(Hz(i)-Hz(i-1)); S-f .NC}:i  
   end [u $X.=(  
   figure(2) gW5yLb_Vz$  
   hold off h-f`as"d  
   plot(x,Ey,'b'); _qxBjB4t"a  
   axis([3*ds L -2 2]); qkM)zOZ^  
   grid on; (GW"iL#.  
   title('E (blue) and 377*H (red)'); :>|dE%/e$  
   hold on oH,{'S@q  
   plot(x,377*Hz,'r'); EV'i/*v}\  
   xlabel('x (m)'); O"GuVC}B  
   pause(0); S4_C8  
   iter 57K\sT4[  
end $} @gR] Z  
if (showWKB==1) b9xvLR8  
 phase = cumsum((epsr).^0.5)*ds; EKD?
j  
 beta0 = 2*pi*fs/(300e6); 8.!+Hm4  
   theory = sin(2*pi*fs*(Niter+4)*dt-beta0*phase)./(epsr.^0.25); /szwVA  
   .. n)Zu>  
j* ZU}Ss  

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