按照sheen经典文章分析矩形贴片天线

我的结果一点都不对.
初学,请大侠们指点:)
源程序如下:
%  fundamental constants
cc=3e8;
muz=4.0*pi*1.0e-7;
epsz=1.0/(cc*cc*muz);
% grid parameters
ie=60;   %number of grid cells in x-direction
je=100;  %number of grid cells in y-direction
ke=16;   %number of grid cells in z-direction
ib=ie+1;
jb=je+1;
kb=ke+1;
dx=0.389e-3;
dy=0.4e-3;
dz=0.265e-3;
dt=0.441e-12;
nmax=1000;
% gaussian pulse excition
T=15e-12;
tau=T/dt;
ndelay=3*tau;
%material parameters
eps=[2.2  (2.2+1)/2  1.0];
%updating coefficients
tux=dt/(muz*dx);
tuy=dt/(muz*dy);
tuz=dt/(muz*dz);
tex1=dt/(epsz*eps(1)*dx);
tey1=dt/(epsz*eps(1)*dy);
tez1=dt/(epsz*eps(1)*dz);
tex2=dt/(epsz*eps(2)*dx);
tey2=dt/(epsz*eps(2)*dy);
tez2=dt/(epsz*eps(2)*dz);
tex3=dt/(epsz*eps(3)*dx);
tey3=dt/(epsz*eps(3)*dy);
tez3=dt/(epsz*eps(3)*dz);
%Field arrays
ex=zeros(ib,jb,kb);
ey=zeros(ib,jb,kb);
ez=zeros(ib,jb,kb);
hx=zeros(ib,jb,kb);
hy=zeros(ib,jb,kb);
hz=zeros(ib,jb,kb);
%*************************************
%begin time-stepping loop
%*************************************
%for n=1:nmax;
for n=1:200
%main grid
%*********************************************
%substrate grid
exnl=ex(2:ie,2,2:ke);
exnr=ex(2:ie,je,2:ke);
eznl=ez(2:ie,2,2:ke);
eznr=ez(2:ie,je,2:ke);
exnt=ex(2:ie,2:je,ke);
eynt=ey(2:ie,2:je,ke);
eznt=ez(2:ie,2:je,ke);
eynf=ey(2,2:je,2:ke);
eynb=ey(ie,2:je,2:ke);
eznf=ez(2,2:je,2:ke);
eznb=ez(ie,2:je,2:ke);
ex(2:ie,2:je,2:3)=ex(2:ie,2:je,2:3)+tey1*(hz(2:ie,3:jb,2:3)-hz(2:ie,2:je,2:3
))-tez1*(hy(2:ie,2:je,3:4)-hy(2:ie,2:je,2:3));
ey(2:ie,2:je,2:3)=ey(2:ie,2:je,2:3)+tez1*(hx(2:ie,2:je,3:4)-hx(2:ie,2:je,2:3
))-tex1*(hz(3:ib,2:je,2:3)-hz(2:ie,2:je,2:3));
ez(2:ie,2:je,2:3)=ez(2:ie,2:je,2:3)+tex1*(hy(3:ib,2:je,2:3)-hy(2:ie,2:je,2:3
))-tey1*(hx(2:ie,3:jb,2:3)-hx(2:ie,2:je,2:3));
hx(2:ie,2:je,2:3)=hx(2:ie,2:je,2:3)+tuz*(ey(2:ie,2:je,2:3)-ey(2:ie,2:je,1:2)
)-tuy*(ez(2:ie,2:je,2:3)-ez(2:ie,1:je-1,2:3));
hy(2:ie,2:je,2:3)=hy(2:ie,2:je,2:3)+tux*(ez(2:ie,2:je,2:3)-ez(1:ie-1,2:je,2:
3))-tuz*(ex(2:ie,2:je,2:3)-ex(2:ie,2:je,1:2));
hz(2:ie,2:je,2:3)=hz(2:ie,2:je,2:3)+tuy*(ex(2:ie,2:je,2:3)-ex(2:ie,1:je-1,2:
3))-tux*(ey(2:ie,2:je,2:3)-ey(1:ie-1,2:je,2:3));
%dielectric-air interface
ex(2:ie,2:je,4)=ex(2:ie,2:je,4)+tey2*(hz(2:ie,3:jb,4)-hz(2:ie,2:je,4))-tez2*
(hy(2:ie,2:je,5)-hy(2:ie,2:je,4));
ey(2:ie,2:je,4)=ey(2:ie,2:je,4)+tez2*(hx(2:ie,2:je,5)-hx(2:ie,2:je,4))-tex2*
(hz(3:ib,2:je,4)-hz(2:ie,2:je,4));
ez(2:ie,2:je,4)=ez(2:ie,2:je,4)+tex2*(hy(3:ib,2:je,4)-hy(2:ie,2:je,4))-tey2*
(hx(2:ie,3:jb,4)-hx(2:ie,2:je,4));
hx(2:ie,2:je,4)=hx(2:ie,2:je,4)+tuz*(ey(2:ie,2:je,4)-ey(2:ie,2:je,3))-tuy*(e
z(2:ie,2:je,4)-ez(2:ie,1:je-1,4));
hy(2:ie,2:je,4)=hy(2:ie,2:je,4)+tux*(ez(2:ie,2:je,4)-ez(1:ie-1,2:je,4))-tuz*
(ex(2:ie,2:je,4)-ex(2:ie,2:je,3));
hz(2:ie,2:je,4)=hz(2:ie,2:je,4)+tuy*(ex(2:ie,2:je,4)-ex(2:ie,1:je-1,4))-tux*
(ey(2:ie,2:je,4)-ey(1:ie-1,2:je,4));
%conductor treatment
ex(16:46,12:50,4)=0;
ey(16:46,12:50,4)=0;
hz(16:46,12:50,4)=0;
ex(21:25,51:99,4)=0;
ey(21:25,51:99,4)=0;
hz(21:25,51:99,4)=0;
%air grid
ex(2:ie,2:je,5:ke)=ex(2:ie,2:je,5:ke)+tey3*(hz(2:ie,3:jb,5:ke)-hz(2:ie,2:je,
5:ke))-tez3*(hy(2:ie,2:je,6:kb)-hy(2:ie,2:je,5:ke));
ey(2:ie,2:je,5:ke)=ey(2:ie,2:je,5:ke)+tez3*(hx(2:ie,2:je,6:kb)-hx(2:ie,2:je,
5:ke))-tex3*(hz(3:ib,2:je,5:ke)-hz(2:ie,2:je,5:ke));
ez(2:ie,2:je,5:ke)=ez(2:ie,2:je,5:ke)+tex3*(hy(3:ib,2:je,5:ke)-hy(2:ie,2:je,
5:ke))-tey3*(hx(2:ie,3:jb,5:ke)-hx(2:ie,2:je,5:ke));
hx(2:ie,2:je,5:ke)=hx(2:ie,2:je,5:ke)+tuz*(ey(2:ie,2:je,5:ke)-ey(2:ie,2:je,4
:ke-1))-tuy*(ez(2:ie,2:je,5:ke)-ez(2:ie,1:je-1,5:ke));
hy(2:ie,2:je,5:ke)=hy(2:ie,2:je,5:ke)+tux*(ez(2:ie,2:je,5:ke)-ez(1:ie-1,2:je
,5:ke))-tuz*(ex(2:ie,2:je,5:ke)-ex(2:ie,2:je,4:ke-1));
hz(2:ie,2:je,5:ke)=hz(2:ie,2:je,5:ke)+tuy*(ex(2:ie,2:je,5:ke)-ex(2:ie,1:je-1
,5:ke))-tux*(ey(2:ie,2:je,5:ke)-ey(1:ie-1,2:je,5:ke));
%*************************************
%ground k=1
ex(2:ie,2:je,1)=0;
ey(2:ie,2:je,1)=0;
hz(2:ie,2:je,1)=0;
%left j=1
ca=(cc*dt-dy)/(cc*dt+dy);
ex(2:ie,1,2:ke)=exnl+ca*(ex(2:ie,2,2:ke)-ex(2:ie,1,2:ke));
ez(2:ie,1,2:ke)=eznl+ca*(ez(2:ie,2,2:ke)-ez(2:ie,1,2:ke));
%right j=jb
hx(2:ie,jb,2:kb)=-hx(2:ie,je,2:kb);
hz(2:ie,jb,2:kb)=-hz(2:ie,je,2:kb);
ez(20:26,jb,1:4)=exp(-((n-ndelay)^2/tau^2));
%top k=kb
cz=(cc*dt-dz)/(cc*dt+dz);
ex(2:ie,2:je,kb)=exnt+cz*(ex(2:ie,2:je,ke)-ex(2:ie,2:je,kb));
ey(2:ie,2:je,kb)=eynt+cz*(ey(2:ie,2:je,ke)-ey(2:ie,2:je,kb));
ez(2:ie,2:je,kb)=ez(2:ie,2:je,kb)+tex1*(hy(3:ib,2:je,kb)-hy(2:ie,2:je,kb))+t
ey1*(hx(2:ie,3:jb,kb)-hx(2:ie,2:je,kb));
hx(2:ie,2:je,kb)=hx(2:ie,2:je,kb)+tuz*(ey(2:ie,2:je,kb)-ey(2:ie,2:je,ke))-tu
y*(ez(2:ie,2:je,kb)-ez(2:ie,1:je-1,kb));
hy(2:ie,2:je,kb)=hy(2:ie,2:je,kb)+tux*(ez(2:ie,2:je,kb)-ez(1:ie-1,2:je,kb))-
tuz*(ex(2:ie,2:je,kb)-ex(2:ie,2:je,ke));
%hz(2:ib,2:jb,kb)=hz(2:ib,2:jb,kb)+tuy*(ex(2:ib,2:jb,kb)-ex(2:ib,1:je,kb))-t
ux*(ey(2:ib,2:jb,kb)-ey(1:ie,2:jb,kb));
%front i=1
cb=(cc*dt-dx)/(cc*dt+dx);
ey(1,2:je,2:ke)=eynf+cb*(ey(2,2:je,2:ke)-ex(1,2:je,2:ke));
ez(1,2:je,2:ke)=eznf+cb*(ez(2,2:je,2:ke)-ez(1,2:je,2:ke));
%back i=ib
ey(ib,2:je,2:ke)=eynb+cb*(ey(ie,2:je,2:ke)-ey(ib,2:je,2:ke));
ez(ib,2:je,2:ke)=eznb+cb*(ez(ie,2:je,2:ke)-ez(ib,2:je,2:ke));
ex(ib,2:je,2:3)=ex(ib,2:je,2:3)+tey1*(hz(ib,3:jb,2:3)-hz(ib,2:je,2:3))-tez1*
(hy(ib,2:je,3:4)-hy(ib,2:je,2:3));
ex(ib,2:je,4)=ex(ib,2:je,4)+tey2*(hz(ib,3:jb,4)-hz(ib,2:je,4))-tez2*(hy(ib,2
:je,5)-hy(ib,2:je,4));
ex(ib,2:je,5:ke)=ex(ib,2:je,5:ke)+tey3*(hz(ib,3:jb,5:ke)-hz(ib,2:je,5:ke))-t
ez3*(hy(ib,2:je,6:kb)-hy(ib,2:je,5:ke));
hy(ib,2:je,2:ke)=hy(ib,2:je,2:ke)+tux*(ez(ib,2:je,2:ke)-ez(ie,2:je,2:ke))-tu
z*(ex(ib,2:je,2:ke)-ez(ib,2:je,1:ke-1));
hz(ib,2:je,2:ke)=hz(ib,2:je,2:ke)+tuy*(ex(ib,2:je,2:ke)-ey(ib,1:je-1,2:ke))-
tux*(ey(ib,2:je,2:ke)-ey(ie,2:je,2:ke));
%end
end;
%x=1:ib;
%y=1:jb;
%Z=abs(ez(x,y,4));
%mesh(Z);

申明：网友回复良莠不齐，仅供参考。如需专业解答，请学习本站推出的微波射频专业培训课程。