有限时域差分

clear
clc

% Fundamental constants

cc=2.99792458e8; %speed of light in free space【光速】
muz=4.0*pi*1.0e-7; %permeability of free space【磁导率】
epsz=1.0/(cc*cc*muz); %permittivity of free space【介电常数】

% Grid parameters

XGridNum=300; %number of grid cells in x-direction【x方向单元格数目】
YGridNum=200; %number of grid cells in y-direction【y方向单元格数目】
ZGridNum=10; %number of grid cells in z-direction【z方向单元格数目】

XB=XGridNum+1;
YB=YGridNum+1;
ZB=ZGridNum+1;

XS=150; %location of z-directed current source 【z方向电流激励源的位置x】
YS=100; %location of z-directed current source 【z方向电流激励源的位置y】

kobs=5;

dx=0.002; %space increment of cubic lattice【空间立方格步进】
dt=dx/(2.0*cc); %time step【时间步进】

nmax=500; %total number of time steps【总步数】

% Differentiated Gaussian pulse excitation【不同的高斯脉冲激励】
% J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2)

rtau=50.0e-12;
tau=rtau/dt;
ndelay=3*tau;
srcconst=-dt*3.0e+11;

% Material parameters【材料参数】

eps=1.0;%真空中的介电常数
sig=0.0;%真空中的磁导率

% Updating coefficients【更新系数】

ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps));
cb=(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps));
da=1.0;
db=dt/muz/dx;

% Field arrays

ex=zeros(XGridNum,YB,ZB);
ey=zeros(XB,YGridNum,ZB);
ez=zeros(XB,YB,ZGridNum);
hx=zeros(XB,YGridNum,ZGridNum);
hy=zeros(XGridNum,YB,ZGridNum);
hz=zeros(XGridNum,YGridNum,ZB);

% Movie initialization

tview(:,:)=ez(:,:,kobs);
sview(:,:)=ez(:,YS,:);
subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview');
shading flat;
caxis([-1.0 1.0]);
colorbar;%【显示右侧颜色示意条】
axis image;%【添加图形坐标】？
title(['Ez(i,j,k=5), time step = 0']);%【添加图形标题】
xlabel('i coordinate');%【添加i坐标显示】
ylabel('j coordinate');%【添加j坐标显示】

subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview');
shading flat;%【删除网格】
caxis([-1.0 1.0]);
colorbar;
axis image;
title(['Ez(i,j=13,k), time step = 0']);
xlabel('i coordinate');
ylabel('k coordinate');

rect=get(gcf,'Position');
rect(1:2)=[0 0];%画图
M=moviein(nmax/2,gcf,rect);

% BEGIN TIME-STEPPING LOOP【时域循环】

for n=1:nmax

% Update electric fields【更新电场】

ex(1:XGridNum,2:YGridNum,2:ZGridNum)=ca*ex(1:XGridNum,2:YGridNum,2:ZGridNum)+...
cb*(hz(1:XGridNum,2:YGridNum,2:ZGridNum)-hz(1:XGridNum,1:YGridNum-1,2:ZGridNum)+...
hy(1:XGridNum,2:YGridNum,1:ZGridNum-1)-hy(1:XGridNum,2:YGridNum,2:ZGridNum));



ey(2:XGridNum,1:YGridNum,

2:ZGridNum)=ca*ey(2:XGridNum,1:YGridNum,2:ZGridNum)+...
cb*(hx(2:XGridNum,1:YGridNum,2:ZGridNum)-hx(2:XGridNum,1:YGridNum,1:ZGridNum-1)+...
hz(1:XGridNum-1,1:YGridNum,2:ZGridNum)-hz(2:XGridNum,1:YGridNum,2:ZGridNum));



ez(2:XGridNum,2:YGridNum,1:ZGridNum)=ca*ez(2:XGridNum,2:YGridNum,1:ZGridNum)+...
cb*(hx(2:XGridNum,1:YGridNum-1,1:ZGridNum)-hx(2:XGridNum,2:YGridNum,1:ZGridNum)+...
hy(2:XGridNum,2:YGridNum,1:ZGridNum)-hy(1:XGridNum-1,2:YGridNum,1:ZGridNum));


ez(XS,YS,1:ZGridNum)=ez(XS,YS,1:ZGridNum)+...
srcconst*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));%加点源

fprintf('/n/n sourse = %d \n',srcconst*(n-ndelay)*exp(-((n-ndelay)^2/tau^2)));
for tt=1:ZGridNum
fprintf('%d \t',ez(24,10,tt));%每一层这一点的值
end

% Update magnetic fields【更新磁场】

hx(2:XGridNum,1:YGridNum,1:ZGridNum)=hx(2:XGridNum,1:YGridNum,1:ZGridNum)+...
db*(ey(2:XGridNum,1:YGridNum,2:ZB)-ey(2:XGridNum,1:YGridNum,1:ZGridNum)+...
ez(2:XGridNum,1:YGridNum,1:ZGridNum)-ez(2:XGridNum,2:YB,1:ZGridNum));

hy(1:XGridNum,2:YGridNum,1:ZGridNum)=hy(1:XGridNum,2:YGridNum,1:ZGridNum)+...
db*(ex(1:XGridNum,2:YGridNum,1:ZGridNum)-ex(1:XGridNum,2:YGridNum,2:ZB)+...
ez(2:XB,2:YGridNum,1:ZGridNum)-ez(1:XGridNum,2:YGridNum,1:ZGridNum));

hz(1:XGridNum,1:YGridNum,2:ZGridNum)=hz(1:XGridNum,1:YGridNum,2:ZGridNum)+...
db*(ex(1:XGridNum,2:YB,2:ZGridNum)-ex(1:XGridNum,1:YGridNum,2:ZGridNum)+...
ey(1:XGridNum,1:YGridNum,2:ZGridNum)-ey(2:XB,1:YGridNum,2:ZGridNum));

%【注意】本程序没有吸收边界条件；没有数值稳定性分析；
% Visualize fields【观察场分布】

if mod(n,2)==0;%求余

timestep=int2str(n);%将数字转换为字符
tview(:,:)=ez(:,:,kobs);
sview(:,:)=ez(:,YS,:);

TT=ez(5,5,5);%无用

subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview');%为指定位置开辟子图，使之成为当前图，并且归一化
shading flat;
caxis([-1.0 1.0]);
colorbar;
axis image;
title(['Ez(i,j,k=5), time step = ',timestep]);
xlabel('i coordinate');
ylabel('j coordinate');

subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview');
shading flat;
caxis([-1.0 1.0]);
colorbar;
axis image;
title(['Ez(i,j=13,k), time step = ',timestep]);
xlabel('i coordinate');
ylabel('k coordinate');

nn=n/2;%取中间进行设置
M(:,nn)=getframe(gcf,rect);%绘图设定颜色

end;

% END TIME-STEPPING LOOP

end

movie(gcf,M,0,10,rect);