4217阅读
15回复

[资料共享]程序__方柱散射

离线weixiao2050

UID ：11674

注册：2008-04-26
登录：2024-09-14
发帖：85
等级：仿真一级

0楼发表于: 2008-07-11 11:59:36

— 本帖被 casey 执行合并操作(2008-07-11) —

%%%%%%%%%%PML%%%%%%%%
for i=1:N
gi2(i)=1;
gi3(i)=1;
fi1(i)=0;
fi2(i)=1;
fi3(i)=1;
end
for j=1:N
gj2(j)=1;
gj3(j)=1;
fj1(j)=0;
fj2(j)=1;
fj3(j)=1;
end
for i=1:npml+1
xnum=npml-i+1;
xxn=xnum/npml;
xn=0.33*(xxn^3);
gi2(i)=1/(1+xn);
gi2(N-i+1)=1/(1+xn);
gi3(i)=(1-xn)/(1+xn);
gi3(N-i+1)=(1-xn)/(1+xn);
xxn=(xnum-0.5)/npml;
xn=0.25*(xxn^3);
fi1(i)=xn;
fi1(N-i)=xn;
fi2(i)=1/(1+xn);
fi2(N-i)=1/(1+xn);
fi3(i)=(1-xn)/(1+xn);
fi3(N-i)=(1-xn)/(1+xn);
end
for j=1:npml+1
xnum=npml-j+1;
xxn=xnum/npml;
xn=0.33*(xxn^3);
gj2(j)=1/(1+xn);
gj2(N-j+1)=1/(1+xn);
gj3(j)=(1-xn)/(1+xn);
gj3(N-j+1)=(1-xn)/(1+xn);
xxn=(xnum-0.5)/npml;
xn=0.25*(xxn^3);
fj1(j)=xn;
fj1(N-j)=xn;
fj2(j) ..

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

离线weixiao2050

UID ：11674

注册：2008-04-26
登录：2024-09-14
发帖：85
等级：仿真一级

1楼发表于: 2008-07-11 12:00:26

程序

以下是以一个方柱为例来仿真其散射的波形：（只限制于近场）
程序如下：
clear;clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c=3*10^8; % 波速
f=3*10^9; % 频率
lamda=c/f; % 波长
k=2*pi/lamda; % 波数
epsz=1/(4*pi*9*10^9); % 真空介电常数
mu=4*pi*10.^(-7); % 真空磁导率
Z=120*pi; % 真空波阻抗
epsilon=1; % 相对介电常数
sigma=0; % 电导率
N=100; % 网格数量
L=800; % 迭代次数
ddx=lamda/20; % 网格尺寸
dt=ddx/(2*c); % 时间间隔
ia=N/4; % 总场区域x左
ib=3*N/4; % 总场区域x右
ja=ia; % 总场区域x下
jb=ib; % 总场区域x上
npml=N/8; % PML点数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r=2*lamda;
M=r/(2*ddx);
for i=1:N
for j=1:N
ga(i,j)=1/(epsilon+sigma*dt/epsz); % 求和参量
gb(i,j)=sigma*dt/epsz; % 求和参量
end
end
spread=8; % 脉冲宽度
t0=25; % 脉冲高度
is=N/2; % 源的X位置
js=N/2; % 源的Y位置
ez_inc=zeros(1,N);
hx_inc=zeros(1,N);
ez_low1=0;
ez_low2=0;
ez_low3=0;
ez_low4=0;
dz=zeros(N,N); % z方向电荷密度
ez=zeros(N,N); % z方向电场
iz=zeros(N,N); % z方向电场求和参量
hx=zeros(N,N); % x方向磁场
hy=zeros(N,N); % y方向磁场
ihx=zeros(N,N); % x方向磁场参量
ihy=zeros(N,N); % y方向磁场参量
%%%%%%%%%%%%%%%%%%%%%%%%PML%%%%%%%%%%%%%%%%%%%
for i=1:N
gi2(i)=1;
gi3(i)=1;
fi1(i)=0;
fi2(i)=1;
fi3(i)=1;
end
for j=1:N
gj2(j)=1;
gj3(j)=1;
fj1(j)=0;
fj2(j)=1;
fj3(j)=1;
end
for i=1:npml+1
xnum=npml-i+1;
xxn=xnum/npml;
xn=0.33*(xxn^3);
gi2(i)=1/(1+xn);
gi2(N-i+1)=1/(1+xn);
gi3(i)=(1-xn)/(1+xn);
gi3(N-i+1)=(1-xn)/(1+xn);
xxn=(xnum-0.5)/npml;
xn=0.25*(xxn^3);
fi1(i)=xn;
fi1(N-i)=xn;
fi2(i)=1/(1+xn);
fi2(N-i)=1/(1+xn);
fi3(i)=(1-xn)/(1+xn);
fi3(N-i)=(1-xn)/(1+xn);
end
for j=1:npml+1
xnum=npml-j+1;
xxn=xnum/npml;
xn=0.33*(xxn^3);
gj2(j)=1/(1+xn);
gj2(N-j+1)=1/(1+xn);
gj3(j)=(1-xn)/(1+xn);
gj3(N-j+1)=(1-xn)/(1+xn);
xxn=(xnum-0.5)/npml;
xn=0.25*(xxn^3);
fj1(j)=xn;
fj1(N-j)=xn;
fj2(j)=1/(1+xn);
fj2(N-j)=1/(1+xn);
fj3(j)=(1-xn)/(1+xn);
fj3(N-j)=(1-xn)/(1+xn);
end
for T=1:L
%%%%%%%%%%%%%%TM波Y方向传播%%%%%%%%%%%%%%%%%%%%%
for j=2:N-1
ez_inc(j)=ez_inc(j)+0.5*(hx_inc(j-1)-hx_inc(j));
end
ez_inc(1)=ez_low2;
ez_low2=ez_low1;
ez_low1=ez_inc(2);
ez_inc(N)=ez_low3;
ez_low3=ez_low4;
ez_low4=ez_inc(N-1);
%% %%%%%%%%%%%%%% 电荷密度dz%%%%%%%%%%%%%%%%%%%%%%%%%
for i=2:N-1
for j=2:N
dz(i,j)=gi3(i)*gj3(j)*dz(i,j)+gi2(i)*gj2(j)*0.5*( hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1) );
end
end
%%%%%%%%%%%%%%%%%%脉冲的加入%%%%%%%%%%%%%%%%%%%
source=exp((-0.5)*( (t0-T)/spread ).^2);
ez_inc(5)=source;
for i=ia:ib
dz(i,ja)=dz(i,ja)+0.5*hx_inc(ja-1);
dz(i,jb)=dz(i,jb)-0.5*hx_inc(jb);
end
%%%%%%%%%%%%%% 电场ez%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:N
for j=1:N
ez(i,j)=ga(i,j)*( dz(i,j)-iz(i,j) );
iz(i,j)=iz(i,j)+gb(i,j)*ez(i,j) ;
end
end
%%%%%%%%%%%%%%%%%%%%边界电场%%%%%%%%%%%%%%%%%%%%%%%
for j=1:N
ez(1,j)=0;
ez(N,j)=0;
end
for i=1:N
ez(i,1)=0;
ez(i,N)=0;
end
%%%%%%%%%%%%%%%%%%%%%%边界条件%%%%%%%%%%%%%%%%%%%%%%%%
for i=N/2-M:N/2+M-1
for j=N/2-M:N/2+M-1
ez(i,j)=0;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%磁场%%%%%%%%%%%%%%%%%%%%%%%%
for j=1:N-1
hx_inc(j)=hx_inc(j)+0.5*(ez_inc(j)-ez_inc(j+1));
end
%%%%%%%%%%%%%%%%%%%%%%%%磁场Hx%%%%%%%%%%%%%%%%%%%%%%
for i=1:N
for j=1:N-1
curl_e=ez(i,j)-ez(i,j+1);
ihx(i,j)=ihx(i,j)+fi1(i)*curl_e;
hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*0.5*(curl_e+ihx(i,j));
end
end;
for i=ia:ib
hx(i,ja-1)=hx(i,ja-1)+0.5*ez_inc(ja);
hx(i,jb)=hx(i,jb)-0.5*ez_inc(jb);
end
%%%%%%%%磁场Hy%%%%%%%%%%%%%%%%%%%
for i=1:N-1
for j=1:N
curl_e=ez(i+1,j)-ez(i,j);
ihy(i,j)=ihy(i,j)+fj1(j)*curl_e;
hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5*(curl_e+ihy(i,j));
end
end
for j=ja:jb;
hy(ia-1,j)=hy(ia-1,j)-0.5*ez_inc(j);
hy(ib,j)=hy(ib,j)+0.5*ez_inc(j);
end
mesh(ez)
drawnow;
end

离线weixiao2050

UID ：11674

注册：2008-04-26
登录：2024-09-14
发帖：85
等级：仿真一级

2楼发表于: 2008-07-11 12:08:34

程序

clear
clc
IE=40;
JE=40;
KE=40;
ic=IE/2;
jc=JE/2;
ia=7;
ja=7;
ka=7;
ib=IE-ia;
jb=JE-ja;
kb=KE-ka;
epsilon=ones(1,10);
sigma=zeros(1,10);
syms n k kc ke nsteps T t0 spread pulse kstart epsilon;
epsz=8.85419e-12;%空气介电常数
ddx=0.1; %仿真步长
dt=ddx/(2*3e8); %计算时间步长
ez_inc=zeros(1,JE);
hx_inc=zeros(1,JE);
ex_low_m1=0; %完全吸收边界
ex_low_m2=0;
ex_high_m2=0;
ex_high_m1=0;
hx=zeros(IE,JE,KE);
hy=zeros(IE,JE,KE);
hz=zeros(IE,JE,KE);
ex=zeros(IE,JE,KE);
ey=zeros(IE,JE,KE);
ez=zeros(IE,JE,KE);
dx=zeros(IE,JE,KE);
dy=zeros(IE,JE,KE);
dz=zeros(IE,JE,KE);
gax=ones(IE,JE,KE); %求解电场参数，包含介电常数，电导率等
gay=ones(IE,JE,KE);
gaz=ones(IE,JE,KE);
gbx=zeros(IE,JE,KE); %求解电场参数，包含介电常数，电导率等
gby=zeros(IE,JE,KE);
gbz=zeros(IE,JE,KE);
idxl=zeros(IE,JE,KE);
idxh=zeros(IE,JE,KE);
ihxl=zeros(IE,JE,KE);
ihxh=zeros(IE,JE,KE);
idyl=zeros(IE,JE,KE);
idyh=zeros(IE,JE,KE);
ihyl=zeros(IE,JE,KE);
ihyh=zeros(IE,JE,KE);
idzl=zeros(IE,JE,KE);
idzh=zeros(IE,JE,KE);
ihzl=zeros(IE,JE,KE);
ihzh=zeros(IE,JE,KE);
gi1=zeros(1,IE);
fi1=zeros(1,IE);
gi2=ones(1,IE);
fi2=ones(1,IE);
gi3=ones(1,IE);
fi3=ones(1,IE);
gj1=zeros(1,JE);
fj1=zeros(1,JE);
gj2=ones(1,JE);
fj2=ones(1,JE);
gj3=ones(1,JE);
fj3=ones(1,JE);
gk1=zeros(1,KE);
fk1=zeros(1,KE);
gk2=ones(1,KE);
fk2=ones(1,KE);
gk3=ones(1,KE);
fk3=ones(1,KE);
npml=input('npml-->');%匹配层宽度
n_pml=npml;
%*********************************************************************
% PML参数设置
%*********************************************************************
for i=1:n_pml
xnum=npml-i+1;
xd=npml;
xxn=xnum/xd;
xn=0.33*power(xxn,3);
fi1(i)=xn;
fi1(IE-i+1)=xn;
gi2(i)=1/(1+xn);
gi2(IE-1+i)=1/(1+xn);
gi3(i)=(1-xn)/(1+xn);
gi3(IE-i+1)=(1-xn)/(1+xn);
xxn=(xnum-0.5)/xd;
xn=0.33*power(xxn,3);
gi1(i)=xn;
gil(IE-i)=xn;
fi2(i)=1/(1+xn);
fi2(IE-i)=1/(1+xn);
fi3(i)=(1-xn)/(1+xn);
fi3(IE-i)=(1-xn)/(1+xn);
end
for j=1:n_pml
xnum=npml-j+1;
xd=npml;
xxn=xnum/xd;
xn=0.33*power(xxn,3);
fj1(j)=xn;
fj1(JE-j+1)=xn;
gj2(j)=1/(1+xn);
gj2(JE-j+1)=1/(1+xn);
gj3(j)=(1-xn)/(1+xn);
gj3(JE-j+1)=(1-xn)/(1+xn);
xxn=(xnum-0.5)/xd;
xn=0.33*power(xxn,3);
gj1(j)=xn;
gjl(JE-j)=xn;
fj2(j)=1/(1+xn);
fj2(JE-j)=1/(1+xn);
fj3(j)=(1-xn)/(1+xn);
fj3(JE-j)=(1-xn)/(1+xn);
end
for k=1:n_pml
xnum=npml-k+1;
xd=npml;
xxn=xnum/xd;
xn=0.33*power(xxn,3);
fk1(k)=xn;
fk1(KE-k+1)=xn;
gk2(k)=1/(1+xn);
gk2(KE-k+1)=1/(1+xn);
gk3(k)=(1-xn)/(1+xn);
gk3(KE-k+1)=(1-xn)/(1+xn);
xxn=(xnum-0.5)/xd;
xn=0.33*power(xxn,3);
gk1(k)=xn;
gkl(KE-k)=xn;
fk2(k)=1/(1+xn);
fk2(KE-k)=1/(1+xn);
fk3(k)=(1-xn)/(1+xn);
fk3(KE-k)=(1-xn)/(1+xn);
end
%*******************************************************
%
%*******************************************************

t0=40; %脉冲中心
spread=10; %电磁脉冲宽度
T=0;
nsteps=1;
while nsteps>0
nsteps=input('nsteps-->');
%fprintf('nsteps is%4.0f\n',nsteps);

for n=1:nsteps
T=T+1;
%************************ main loop ***************************
%************************ DX *****************************
for i=1:ia
for j=2:JE
for k=2:KE
curl_h=(hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1));
idxl(i,j,k)=idxl(i,j,k)+curl_h;
dx(i,j,k)=gj3(j)*gk3(k)*dx(i,j,k)+gj2(j)*gk2(k)*0.5*(gi1(i)*idxl(i,j,k)+curl_h);
%dx(i,j,k)=gj3(j)*gk3(k)*dx(i,j,k)+gj2(j)*gk2(k)*0.5*(curl_h);
end
end
end
for i=(ia+1):ib
for j=2:JE
for k=2:KE
curl_h=(hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1));
dx(i,j,k)=gj3(j)*gk3(k)*dx(i,j,k)+gj2(j)*gk2(k)*0.5*(curl_h);
end
end
end
for i=(ib+1):IE
ixh=i-ib;
for j=2:JE
for k=2:KE
curl_h=(hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1));
idxh(ixh,j,k)=idxh(ixh,j,k)+curl_h;
dx(i,j,k)=gj3(j)*gk3(k)*dx(i,j,k)+gj2(j)*gk2(k)*0.5*(gi1(i)*idxh(ixh,j,k)+curl_h);
%dx(i,j,k)=gj3(j)*gk3(k)*dx(i,j,k)+gj2(j)*gk2(k)*0.5*(curl_h);
end
end
end
%******************************* DY ***********************************
for i=2:IE
for j=1:ja
for k=2:KE
curl_h=(hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k));
idyl(i,j,k)=idyl(i,j,k)+curl_h;
dy(i,j,k)=gi3(i)*gk3(k)*dy(i,j,k)+0.5*gi2(i)*gk2(k)*(gj1(j)*idyl(i,j,k)+curl_h);
%dy(i,j,k)=gi3(i)*gk3(k)*dy(i,j,k)+0.5*gi2(i)*gk2(k)*(curl_h);
end
end
end
for i=2:IE
for j=(ja+1):jb
for k=2:KE
curl_h=(hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k));
dy(i,j,k)=gi3(i)*gk3(k)*dy(i,j,k)+0.5*gi2(i)*gk2(k)*(curl_h);
end
end
end
for i=2:IE
for j=(jb+1):JE
jyh=j-jb;
for k=2:KE
curl_h=(hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k));
idyh(i,jyh,k)=idyh(i,jyh,k)+curl_h;
dy(i,j,k)=gi3(i)*gk3(k)*dy(i,j,k)+0.5*gi2(i)*gk2(k)*(gj1(j)*idyh(i,jyh,k)+curl_h);
end
end
end

%******************************* DZ ***********************************

for i=2:IE
for j=2:JE
for k=1:ka
curl_h=(hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k));
idzl(i,j,k)=idzl(i,j,k)+curl_h;
dz(i,j,k)=gi3(i)*gj3(j)*dz(i,j,k)+0.5*gi2(i)*gj2(j)*(gk1(k)*idzl(i,j,k)+curl_h);
end
end
end
for i=2:IE
for j=2:JE
for k=(ka+1):kb
curl_h=(hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k));
dz(i,j,k)=gi3(i)*gj3(j)*dz(i,j,k)+0.5*gi2(i)*gj2(j)*(curl_h);
end
end
end
for i=2:IE
for j=2:JE
for k=(kb+1):KE
kzh=k-kb;
curl_h=(hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k));
idzh(i,j,kzh)=idzh(i,j,kzh)+curl_h;
dz(i,j,k)=gi3(i)*gj3(j)*dz(i,j,k)+0.5*gi2(i)*gj2(j)*(gk1(k)*idzh(i,j,kzh)+curl_h);
end
end
end
%**************************************************************************
% source
%**************************************************************************
for k=(KE/2-2):(KE/2+2)
dz(ic,jc,k)=0;
end
pulse=exp(-0.5*(power((t0-T)/spread,2))); %入射干扰
dz(ic,jc,KE/2)=pulse;

for i=1:IE
for j=1:JE
for k=1:KE
ex(i,j,k)=gax(i,j,k)*dx(i,j,k);
ey(i,j,k)=gay(i,j,k)*dy(i,j,k);
ez(i,j,k)=gaz(i,j,k)*dz(i,j,k);
end
end
end

%**************************************************************************
% H
%**************************************************************************
%******************************* HX ***********************************
for i=1:ia
for j=1:JE-1
for k=1:KE-1
curl_e=(ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k));
ihxl(i,j,k)=ihxl(i,j,k)+curl_e;
hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*(fi1(i)*ihxl(i,j,k)+curl_e);
% hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*(curl_e);
end
end
end
for i=(ia+1):ib
for j=1:JE-1
for k=1:KE-1
curl_e=(ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k));
%ihxl(i,j,k)=ihxl(i,j,k)+curl_e;
%hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*(fi1(i)*ihxl(i,j,k)+curl_e);
hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*(curl_e);
end
end
end
for i=ib+1:IE
ixh=i-ib;
for j=1:JE-1
for k=1:KE-1
curl_e=(ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k));
ihxh(ixh,j,k)=ihxh(ixh,j,k)+curl_e;
hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*(fi1(i)*ihxh(ixh,j,k)+curl_e);
% hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*(curl_e);
end
end
end
%******************************* Hy ***********************************
for i=1:IE-1
for j=1:ja
for k=1:KE-1
curl_e=(ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k));
ihyl(i,j,k)=ihyl(i,j,k)+curl_e;
hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+0.5*fi2(i)*fk2(k)*(fj1(j)*ihyl(i,j,k)+curl_e);
% hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+0.5*fi2(i)*fk2(k)*(curl_e);
end
end
end
for i=1:IE-1
for j=(ja+1):jb
for k=1:KE-1
curl_e=(ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k));
%ihyl(i,j,k)=ihyl(i,j,k)+curl_e;
%hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+0.5*fi2(i)*fk2(k)*(fj1(j)*ihyl(i,j,k)+curl_e);
hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+0.5*fi2(i)*fk2(k)*(curl_e);
end
end
end
for i=1:IE-1
for j=jb+1:JE
jyh=j-jb;
for k=1:KE-1
curl_e=(ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k));
ihyh(i,jyh,k)=ihyh(i,jyh,k)+curl_e;
hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+0.5*fi2(i)*fk2(k)*(fj1(j)*ihyh(i,jyh,k)+curl_e);
%hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+0.5*fi2(i)*fk2(k)*(curl_e);
end
end
end

%******************************* Hz ***********************************
for i=1:IE-1
for j=1:JE-1
for k=1:ka
curl_e=(ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k));
ihzl(i,j,k)=ihzl(i,j,k)+curl_e;
hz(i,j,k)=fi3(i)*fj3(j)*hz(i,j,k)+0.5*fi2(i)*fj2(j)*(fk1(k)*ihzl(i,j,k)+curl_e);
%hz(i,j,k)=fi3(i)*fj3(j)*hz(i,j,k)+0.5*fi2(i)*fj2(j)*(curl_e);
end
end
end
for i=1:IE-1
for j=1:JE-1
for k=(ka+1):kb
curl_e=(ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k));
%ihzl(i,j,k)=ihzl(i,j,k)+curl_e;
%hz(i,j,k)=fi3(i)*fj3(j)*hz(i,j,k)+0.5*fi2(i)*fj2(j)*(fk1(k)*ihzl(i,j,k)+curl_e);
hz(i,j,k)=fi3(i)*fj3(j)*hz(i,j,k)+0.5*fi2(i)*fj2(j)*(curl_e);
end
end
end
for i=1:IE-1
for j=1:JE-1
for k=(kb+1):KE
kzh=k-kb;
curl_e=(ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k));
ihzh(i,j,kzh)=ihzh(i,j,kzh)+curl_e;
hz(i,j,k)=fi3(i)*fj3(j)*hz(i,j,k)+0.5*fi2(i)*fj2(j)*(fk1(k)*ihzh(i,j,kzh)+curl_e);
%hz(i,j,k)=fi3(i)*fj3(j)*hz(i,j,k)+0.5*fi2(i)*fj2(j)*(curl_e);
end
end
end
end

x=1:IE;
y=1:JE;
subplot(1,2,1);
plot(x,ez(x,20,20));
subplot(1,2,2);
meshc(x,y,dz(x,y,20));
%pcolor(x,y,ez(x,y,10))
hidden off;
title('Ez');
xlabel('X');
ylabel('Y');
%zlabel('Z');
grid on
fprintf('T=%5.0f\n',T);
end