clear
IE=60;JE=60;
ga=ones(IE,JE);
dz=zeros(IE,JE);% D
ex=zeros(IE,JE);
hx=zeros(IE,JE);hy=hx;
ic=IE/2;jc=JE/2;%观测点
ddx=0.01;dt=ddx/6e8;
epsz=8.851*1e-12;
%%%%%%%%%%表明谐振子
t0=20;
spread=6;
T=0;
nsteps=1;
%%%%%%%%%%%%%%%电磁场时间步循环
while(nsteps>0)
nsteps = input('输入循环的时间步，小于0则终止: ');
for n=1:nsteps
T=T+1;
%%%%%%%%%%%dz
for j=2:JE
for i=2:IE
dz(i,j)=dz(i,j)+0.5*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1));
end
end
%%%%%%加入激励源
pulse=sin(2*pi*1500*1e6*dt*T);%exp(-0.5*((t0-T)/spread)^2);
dz(ic,jc)=pulse;
%%%%%%%%%%%ex,ey,ez
for j=2:JE
for i=2:IE
ez(i,j)=ga(i,j)*dz(i,j);
end
end
%%%%%%%%%%%hx
for j=1:JE-1
for i=1:IE-1
hx(i,j)=hx(i,j)+0.5*(ez(i,j)-ez(i,j+1));
end
end
%%%%%%%%%%%%%hy
for j=1:JE-1
for i=1:IE-1
hy(i,j)=hy(i,j)+0.5*(ez(i+1,j)-ez(i,j));
end
end
figure(1)
for i=1:IE
for j=1:JE
ezkout(i,j)=(ez(i,j));
end
end
[i,j]=meshgrid(1:IE,1:JE);
mesh(i,j,(ezkout))
surf(i,j,(ezkout))
axis([1,IE,1,JE,-0.5,1]);
set(gca,'FontName','TimesNewRoman');
set(gca,'FontSize',14);
title(['时间步为:',num2str(n),'步']);
pause
end
end
