982阅读
0回复

[讨论]啊，受不了了，我的二维TEfdtd老出错，请各位抽出宝贵的几分钟给我看看行吗？教我一下，谢谢

自强不息厚德载物

UID ：58301

注册：2010-04-28
登录：2022-12-16
发帖：332
等级：仿真三级

0楼发表于: 2011-10-24 11:28:42

%常数定义
ie=80;
je=80;
i0=25;
j0=25;
epsr=8.854e-12;
miu=4*pi*1e-7;
sigma=0;
sigmam=0;
Z0=sqrt(miu/epsr);
ddx=0.01;
dy=ddx;
dt=ddx/(6e8);
ia=4;
ib=ie-4+1;
ja=4;
jb=je-4+1;
%初始化数组
sigmax=0.001;
sigmamx=Z0^2*sigmax;
sigmay=0.001;
sigmamy=Z0^2*sigmay;

ex=zeros(ie,je);
ey=zeros(ie,je);
hz=zeros(ie+1,je+1);
hzx=zeros(ie+1,je+1);
hzy=zeros(ie+1,je+1);
nsteps=300;
ee=zeros(nsteps,1);
t0=12;
spread=25;
T=0;
%主程序
for n=1:nsteps
T=T+1;
%真空中,fdtd迭代
ca=(1.-sigma*dt/(2*epsr))/(1.+sigma*dt/(2*epsr));
cb=dt/epsr/(1.+sigma*dt/(2*epsr));
cp=(1.-sigmam*dt/(2*miu))/(1.+sigmam*dt/(2*miu));
cq=dt/epsr/(1.+sigmam*dt/(2*miu));
cb=dt/epsr;
cq=dt/miu;
ex(ia:ib,ja:jb)=ex(ia:ib,ja:jb)+cb.*(hz(ia:ib,ja:jb)-hz(ia:ib,ja-1:jb-1)/dy);
ey(ia:ib,ja:jb)=ey(ia:ib,ja:jb)+cb.*(hz(ia:ib,ja:jb)-hz(ia-1:ib-1,ja:jb)/ddx);
%source
hz(ia-1:ib-1,ja-1:jb-1)=hz(ia-1:ib-1,ja-1:jb-1)+cq.*((ey(ia:ib,ja:jb)-ey(ia-1:ib-1,ja-1:jb-1))/ddx-(ex(ia:ib,ja:jb)-ex(ia-1:ib-1,ja-1:jb-1))/dy);
pulse=exp(-0.5*(T*dt-t0)^2/spread^2)/10000;
ex(i0,j0)=ex(i0,j0)+pulse;
%pml层顶点计算
e1=(1-exp(-sigmay*dt/epsr))/(dy*sigmay);
% ex在四个角顶点的计算
%左下角，hz ->1:ja
ex(1:ia-1,1:ja-1)=exp(-sigmay*dt/epsr).*ex(1:ia-1,1:ja-1)-e1.*(hz(1:ia-1,1:ja-1)-hz(1:ia-1,2:ja));
%左上角
ex(1:ia-1,jb-1:je-1)=exp(-sigmay*dt/epsr).*ex(1:ia-1,jb-1:je-1)-e1.*(hz(1:ia-1,jb-1:je-1)-hz(1:ia-1,jb:je));
%右下角
ex(ib-1:ie-1,1:ja-1)=exp(-sigmay*dt/epsr).*ex(ib-1:ie-1,1:ja-1)-e1.*(hz(ib-1:ie-1,1:ja-1)-hz(ib-1:ie-1,2:ja));
%右上角
ex(ib-1:ie-1,jb-1:je-1)=exp(-sigmay*dt/epsr).*ex(ib-1:ie-1,jb-1:je-1)-e1.*(hz(ib-1:ie-1,jb-1:je-1)-hz(ib-1:ie-1,jb:je));
% ey在四个角顶点的计算
e11=(1-exp(-sigmax*dt/epsr))/(ddx*sigmax);
%左下角，hz ->1:ia
ey(1:ia-1,1:ja-1)=exp(-sigmax*dt/epsr).*ey(1:ia-1,1:ja-1)-e11.*(hz(2:ia,1:ja-1)-hz(1:ia-1,1:ja-1));
%左上角
ey(1:ia-1,jb-1:je-1)=exp(-sigmax*dt/epsr).*ey(1:ia-1,jb-1:je-1)-e11.*(hz(2:ia,jb-1:je-1)-hz(1:ia-1,jb-1:je-1));
%右下角
ey(ib-1:ie-1,1:ja-1)=exp(-sigmax*dt/epsr).*ey(ib-1:ie-1,1:ja-1)-e11.*(hz(ib:ie,1:ja-1)-hz(ib-1:ie-1,1:ja-1));
%右上角
ey(ib-1:ie-1,jb-1:je-1)=exp(-sigmax*dt/epsr).*ey(ib-1:ie-1,jb-1:je-1)-e11.*(hz(ib:ie,jb-1:je-1)-hz(ib-1:ie-1,jb-1:je-1));

%%%hzx在四个角顶点的计算
hhzx1=(1-exp(-sigmamx*dt/epsr))/(ddx*sigmamx);
%左下角
hzx(2:ia,2:ja)=exp(-sigmamx*dt/miu).*hzx(2:ia,2:ja)-hhzx1.*(ey(2:ia,2:ja)-ey(1:ia-1,2:ja));
%右下角
hzx(ib+1:ie,2:ja)=exp(-sigmamx*dt/miu).*hzx(ib+1:ie,2:ja)-hhzx1.*(ey(ib+1:ie,2:ja)-ey(ib:ie-1,2:ja));
%左上角
hzx(2:ia,jb+1:je)=exp(-sigmamx*dt/miu).*hzx(2:ia,jb+1:je)-hhzx1.*(ey(2:ia,jb+1:je)-ey(1:ia-1,jb+1:je));
%右上角
hzx(ib+1:ie,jb+1:je)=exp(-sigmamx*dt/miu).*hzx(ib+1:ie,jb+1:je)-hhzx1.*(ey(ib+1:ie,jb+1:je)-ey(ib:ie-1,jb+1:je));
%%%hzy在四个角顶点的计算
hhzy1=(1-exp(-sigmamy*dt/epsr))/(ddx*sigmamy);
%左下角
hzy(2:ia,2:ja)=exp(-sigmay*dt/epsr).*hzy(2:ia,2:ja)-hhzy1.*(ex(2:ia,2:ja)-ex(2:ia,1:ja-1));
%左上角
hzy(2:ia,jb+1:je)=exp(-sigmay*dt/epsr).*hzy(2:ia,jb+1:je)-hhzy1.*(ex(2:ia,jb+1:je)-ex(2:ia,jb:je-1));
%右下角
hzy(ib+1:ie,2:ja)=exp(-sigmay*dt/epsr).*hzy(ib+1:ie,2:ja)-hhzy1.*(ex(ib+1:ie,2:ja)-ex(ib+1:ie,1:ja-1));
%右上角
hzy(ib+1:ie,jb+1:je)=exp(-sigmay*dt/epsr).*hzy(ib+1:ie,jb+1:je)-hhzy1.*(ex(ib+1:ie,jb+1:je)-ex(ib+1:ie,jb:je-1));
%pml 左右两边
%%%% the left side
ex(1:ia-1,ja:jb-1)=ex(1:ia-1,ja:jb-1)-Z0/2.*(hz(1:ia-1,ja:jb-1)-hz(1:ia-1,ja+1:jb));
ey(1:ia-1,ja:jb-1)=exp(-sigmax*dt/epsr).*ey(1:ia-1,ja:jb-1)-e11.*(hz(2:ia,ja:jb-1)-hz(1:ia-1,ja:jb-1));
hzx(2:ia,ja+1:jb)=exp(-sigmamx*dt/miu).*hzx(2:ia,ja+1:jb)-hhzx1.*(ey(2:ia,ja+1:jb)-ey(1:ia-1,ja+1:jb));
hzy(2:ia,ja+1:jb)=hzy(2:ia,ja+1:jb)-1/(2*Z0).*(ex(2:ia,ja:jb-1)-ex(2:ia,ja+1:jb));
%%% the right side
ex(ib:ie-1,ja:jb-1)=ex(ib:ie-1,ja:jb-1)-Z0/2.*(hz(ib:ie-1,ja:jb-1)-hz(ib:ie-1,ja+1:jb));
ey(ib:ie-1,ja:jb-1)=exp(-sigmax*dt/epsr).*ey(ib:ie-1,ja:jb-1)-e11.*(hz(ib+1:ie,ja:jb-1)-hz(ib:ie-1,ja:jb-1));
hzx(ib+1:ie,ja+1:jb)=exp(-sigmamx*dt/miu).*hzx(ib+1:ie,ja+1:jb)-hhzx1.*(ey(ib+1:ie,ja+1:jb)-ey(ib:ie-1,ja+1:jb));
hzy(ib+1:ie,ja+1:jb)=hzy(ib+1:ie,ja+1:jb)-1/(2*Z0).*(ex(ib+1:ie,ja:jb-1)-ex(ib+1:ie,ja+1:jb));
%pml 上下两边
%%the bottom side
ex(ia:ib-1,1:ja-1)=exp(-sigmay*dt/epsr).*ex(ia:ib-1,1:ja-1)-e1.*(hz(ia:ib-1,1:ja-1)-hz(ia:ib-1,2:ja));
ey(ia:ib-1,1:ja-1)=ey(ia:ib-1,1:ja-1)-(Z0/2).*(hz(ia+1:ib,1:ja-1)-hz(ia:ib-1,1:ja-1));
hzx(ia+1:ib,2:ja)=hzx(ia+1:ib,2:ja)-1/(2*Z0).*(ey(ia+1:ib,2:ja)-ey(ia:ib-1,2:ja));
hzy(ia+1:ib,2:ja)=exp(-sigmay*dt/epsr).*hzy(ia+1:ib,2:ja)-hhzy1.*(ex(ia+1:ib,2:ja)-ex(ia+1:ib,1:ja-1));
%%the upper side
ex(ia:ib-1,jb:je-1)=exp(-sigmay*dt/epsr).*ex(ia:ib-1,jb:je-1)-e1.*(hz(ia:ib-1,jb:je-1)-hz(ia:ib-1,jb+1:je));
ey(ia:ib-1,jb:je-1)=ey(ia:ib-1,jb:je-1)-(Z0/2).*(hz(ia+1:ib,jb:je-1)-hz(ia:ib-1,jb:je-1));
hzx(ia+1:ib,jb+1:je)=hzx(ia+1:ib,jb+1:je)-1/(2*Z0).*(ey(ia+1:ib,jb+1:je)-ey(ia:ib-1,jb+1:je));
hzy(ia+1:ib,jb+1:je)=exp(-sigmay*dt/epsr).*hzy(ia+1:ib,jb+1:je)-hhzy1.*(ex(ia+1:ib,jb+1:je)-ex(ia+1:ib,jb:je-1));
%真空与pml交界面
ey(ib,ja:jb)=ey(ib,ja:jb)-Z0/2.*(hz(ib+1,ja:jb)-hz(ib,ja:jb));
ex(ia:ib,ja)=ex(ia:ib,ja)-Z0/2.*(hz(ia:ib,ja)-hz(ia:jb,ja+1));

发帖回复