1943阅读
7回复

[求助]fdtd2d_tm_pml为何出现电磁场分量变成无限大的问题

UID ：20795

注册：2008-11-05
登录：2012-06-06
发帖：18
等级：仿真新人

0楼发表于: 2009-08-24 22:40:51

我根据Susan C.Hagness的书所附的2D关于TE波的pml边界吸收的程序简单改写为TM波模式但是出现了电磁场迅速增大的问题。请教下各位大侠，程序附如下

clear

%***********************************************************************
% 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

freq=5.0e+9; %center frequency of source excitation
lambda=cc/freq; %center wavelength of source excitation
omega=2.0*pi*freq;

%***********************************************************************
% Grid parameters
%***********************************************************************

ie=100; %number of grid cells in x-direction
je=50; %number of grid cells in y-direction

ib=ie+1;
jb=je+1;

is=15; %location of z-directed hard source
js=je/2; %location of z-directed hard source

dx=3.0e-3; %space increment of square lattice
dt=dx/(2.0*cc); %time step

nmax=2; %total number of time steps

iebc=8; %thickness of left and right PML region
jebc=8; %thickness of front and back PML region
rmax=0.00001;
orderbc=2;
ibbc=iebc+1;
jbbc=jebc+1;
iefbc=ie+2*iebc;
jefbc=je+2*jebc;
ibfbc=iefbc+1;
jbfbc=jefbc+1;

%***********************************************************************
% Material parameters
%***********************************************************************

media=2;

eps=[1.0 1.0];
sig=[0.0 1.0e+7];
mur=[1.0 1.0];
sim=[0.0 0.0];

%***********************************************************************
% Wave excitation
%***********************************************************************

rtau=160.0e-12;
tau=rtau/dt;
delay=3*tau;

source=zeros(1,nmax);
for n=1:7.0*tau
source(n)=sin(omega*(n-delay)*dt)*exp(-((n-delay)^2/tau^2));
end

%***********************************************************************
% Field arrays
%***********************************************************************
ez=zeros(ie,je);
hx=zeros(ie,jb); %fields in main grid
hy=zeros(ib,je);

ezxbcf=zeros(iefbc,jebc);
ezybcf=zeros(iefbc,jebc);
hxbcf=zeros(iefbc,jebc); %fields in front pml
hybcf=zeros(ibfbc,jebc);

ezxbcb=zeros(iefbc,jebc);
ezybcb=zeros(iefbc,jebc);
hxbcb=zeros(iefbc,jbbc); %fields in back pml
hybcb=zeros(ibfbc,jebc);

ezxbcl=zeros(iebc,je);
ezybcl=zeros(iebc,je);
hxbcl=zeros(iebc,jb); %fields in left pml
hybcl=zeros(iebc,je);

ezxbcr=zeros(iebc,je);
ezybcr=zeros(iebc,je);
hxbcr=zeros(iebc,jb); %fields in right pml
hybcr=zeros(ibbc,je);

%***********************************************************************
% Updating coefficients
%***********************************************************************

for i=1:media
eaf =dt*sig(i)/(2.0*epsz*eps(i));
ca(i)=(1.0-eaf)/(1.0+eaf);
cb(i)=dt/epsz/eps(i)/dx/(1.0+eaf);
haf =dt*sim(i)/(2.0*muz*mur(i));
da(i)=(1.0-haf)/(1.0+haf);
db(i)=dt/muz/mur(i)/dx/(1.0+haf);
end

%***********************************************************************
% Geometry specification (main grid)
%***********************************************************************

% Initialize entire main grid to free space
caez(1:ie,1:je)=ca(1);
cbez(1:ie,1:je)=cb(1);

dahx(1:ie,1:jb)=da(1);
dbhx(1:ie,1:jb)=db(1);

dahy(1:ib,1:je)=da(1);
dbhy(1:ib,1:je)=db(1);

% Add metal cylinder

diam=20; % diameter of cylinder: 6 cm
rad=diam/2.0; % radius of cylinder: 3 cm
icenter=4*ie/5; % i-coordinate of cylinder's center
jcenter=je/2; % j-coordinate of cylinder's center

for i=1:ie
for j=1:je
dist2=(i-icenter)^2 + (j-jcenter)^2;
if dist2 <= rad^2
caez(i,j)=ca(2);
cbez(i,j)=cb(2);
end
end
end

%***********************************************************************
% Fill the PML regions
%***********************************************************************

delbc=iebc*dx;
sigmam=-log(rmax/100.0)*epsz*cc*(orderbc+1)/(2*delbc);
bcfactor=eps(1)*sigmam/(dx*(delbc^orderbc)*(orderbc+1));

% FRONT region
dahxbcf(1:iefbc,1)=1.0;
dbhxbcf(1:iefbc,1)=0.0;
for j=2:jebc
y1=(jebc-j+1.5)*dx;
y2=(jebc-j+0.5)*dx;
sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
sigmays=sigmay*(muz/(epsz*eps(1)));
da1=exp(-sigmays*dt/muz);
db1=(1-da1)/(sigmays*dx);
dahxbcf(1:iefbc,j)=da1;
dbhxbcf(1:iefbc,j)=db1;
end
sigmay=bcfactor*(0.5*dx)^(orderbc+1);
sigmays=sigmay*(muz/(epsz*eps(1)));
da1=exp(-sigmays*dt/muz);
db1=(1-da1)/(sigmays*dx);
dahx(1:ie,1)=da1;
dbhx(1:ie,1)=db1;
dahxbcl(1:iebc,1)=da1;
dbhxbcl(1:iebc,1)=db1;
dahxbcr(1:iebc,1)=da1;
dbhxbcr(1:iebc,1)=db1;

for j=1:iebc
y1=(jebc-j+1)*dx;
y2=(jebc-j)*dx;
sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
ca1=exp(-sigmay*dt/(epsz*eps(1)));
cb1=(1.0-ca1)/(sigmay*dx);
caezybcf(1:iefbc,j)=ca1;
cbezybcf(1:iefbc,j)=cb1;
caezxbcf(1:iefbc,j)=ca(1);
cbezxbcf(1:iefbc,j)=cb(1);
dahybcf(1:ibfbc,j)=da(1);
dbhybcf(1:ibfbc,j)=db(1);
end

% BACK region
dahxbcb(1:iefbc,jbbc)=1.0;
dbhxbcb(1:iefbc,jbbc)=0.0;
for j=2:jebc
y1=(j-0.5)*dx;
y2=(j-1.5)*dx;
sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
sigmays=sigmay*(muz/(epsz*eps(1)));
da1=exp(-sigmays*dt/muz);
db1=(1-da1)/(sigmays*dt);
dahxbcb(1:iefbc,j)=da1;
dbhxbcb(1:iefbc,j)=db1;
end
sigmay=bcfactor*(0.5*dx)^(orderbc+1);
sigmays=sigmay*(muz/(epsz*eps(1)));
da1=exp(-sigmays*dt/muz);
db1=(1-da1)/(sigmays*dx);
dahx(1:ie,jb)=da1;
dbhy(1:ie,jb)=db1;
dahxbcl(1:iebc,jb)=da1;
dbhxbcl(1:iebc,jb)=db1;
dahxbcr(1:iebc,jb)=da1;
dbhxbcr(1:iebc,jb)=db1;

for j=1:jebc
y1=j*dx;
y2=(j-1)*dx;
sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
ca1=exp(-sigmay*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmay*dx);
caezybcb(1:iefbc,j)=ca1;
cbezybcb(1:iefbc,j)=cb1;
caezxbcb(1:iefbc,j)=ca(1);
cbezxbcb(1:iefbc,j)=cb(1);
dahybcb(1:ibfbc,j)=da(1);
dbhybcb(1:ibfbc,j)=db(1);
end

% LEFT region
dahybcl(1,1:je)=1.0;
dbhybcl(1,1:je)=0.0;
for i=2:iebc
x1=(iebc-i+1.5)*dx;
x2=(iebc-i+0.5)*dx;
sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
sigmaxs=sigmax*(muz/(epsz*eps(1)));
da1=exp(-sigmaxs*dt/muz);
db1=(1-da1)/(sigmaxs*dx);
dahybcl(i,1:je)=da1;
dbhybcl(i,1:je)=db1;
dahybcf(i,1:jebc)=da1;
dbhybcf(i,1:jebc)=db1;
dahybcb(i,1:jebc)=da1;
dbhybcb(i,1:jebc)=db1;
end
sigmax=bcfactor*(0.5*dx)^(orderbc+1);
sigmaxs=sigmax*(muz/(epsz*eps(1)));
da1=exp(-sigmaxs*dt/muz);
db1=(1-da1)/(sigmaxs*dx);
dahy(1,1:je)=da1;
dbhy(1,1:je)=db1;
dahybcf(ibbc,1:jebc)=da1;
dbhybcf(ibbc,1:jebc)=db1;
dahybcb(ibbc,1:jebc)=da1;
dbhybcb(ibbc,1:jebc)=db1;

for i=1:iebc
x1=(iebc-i+1)*dx;
x2=(iebc-i)*dx;
sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
ca1=exp(-sigmax*dt/(epsz*eps(1)));
cb1=(1-ca1)*(sigmax*dx);
caezxbcl(i,1:je)=ca1;
cbezxbcl(i,1:je)=cb1;
caezxbcf(i,1:jebc)=ca1;
cbezxbcf(i,1:jebc)=cb1;
caezxbcb(i,1:jebc)=ca1;
cbezxbcb(i,1:jebc)=cb1;

caezybcl(i,1:je)=ca(1);
cbezybcl(i,1:je)=cb(1);
dahxbcl(i,2:je)=da(1);
dbhxbcl(i,2:je)=db(1);
end

% RIGHT region
dahybcr(ibbc,1:je)=1.0;
dbhybcr(ibbc,1:je)=0.0;
for i=2:iebc
x1=(i-0.5)*dx;
x2=(i-1.5)*dx;
sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
sigmaxs=sigmax*(muz/(epsz*eps(1)));
da1=exp(-sigmaxs*dt/muz);
db1=(1-da1)/(sigmaxs*dx);
dahybcr(i,1:je)=da1;
dbhybcr(i,1:je)=db1;
dahybcf(iebc+ie+i,1:jebc)=da1;
dbhybcf(iebc+ie+i,1:jebc)=db1;
dahybcb(iebc+ie+i,1:jebc)=da1;
dbhybcb(iebc+ie+i,1:jebc)=db1;
end
sigmax=bcfactor*(0.5*dx)^(orderbc+1);
sigmaxs=sigmax*(muz/(epsz*eps(1)));
da1=exp(-sigmaxs*dt/muz);
db1=(1-da1)/(sigmaxs*dx);
dahy(ib,1:je)=da1;
dbhy(ib,1:je)=db1;
dahybcf(iebc+ie+1,1:jebc)=da1;
dbhybcf(iebc+ie+1,1:jebc)=db1;
dahybcb(iebc+ie+1,1:jebc)=da1;
dbhybcb(iebc+ie+1,1:jebc)=db1;

for i=1:iebc
x1=i*dx;
x2=(i-1)*dx;
sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
ca1=exp(-sigmax*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmax*dx);
caezxbcr(i,1:je)=ca1;
cbezxbcr(i,1:je)=cb1;
caezxbcf(iebc+ie+i,1:jebc)=ca1;
cbezxbcf(iebc+ie+i,1:jebc)=cb1;
caezxbcb(iebc+ie+i,1:jebc)=ca1;
cbezxbcb(iebc+ie+i,1:jebc)=cb1;

caezybcr(i,1:je)=ca(1);
cbezybcr(i,1:je)=cb(1);
dahxbcr(i,2:je)=da(1);
dbhxbcr(i,2:je)=db(1);
end
%***********************************************************************
% Movie initialization
%***********************************************************************

subplot(3,1,1),pcolor(hx');
shading flat;
caxis([-0.2 0.2]);
axis([1 ie 1 jb]);
colorbar;
axis image;
axis off;
title(['Hx at time step = 0']);

subplot(3,1,2),pcolor(hy');
shading flat;
caxis([-0.2 0.2]);
axis([1 ib 1 je]);
colorbar;
axis image;
axis off;
title(['Hy at time step = 0']);

subplot(3,1,3),pcolor(ez');
shading flat;
caxis([-80.0 80.0]);
axis([1 ie 1 je]);
colorbar;
axis image;
axis off;
title(['Ez at time step = 0']);

rect=get(gcf,'Position');
rect(1:2)=[0 0];

M=moviein(nmax/4,gcf,rect);
%***********************************************************************
% BEGIN TIME-STEPPING LOOP
%***********************************************************************
for n=1:nmax

%***********************************************************************
% Update electric fields (EZ) in main grid
%***********************************************************************
ez(1:ie,1:je)=caez(1:ie,1:je).*ez(1:ie,1:je)+cbez(1:ie,1:je).*(hy(2:ib,1:je)-hy(1:ie,1:je)-...
(hx(1:ie,2:jb)-hx(1:ie,1:je)));
ez(is,js)=source(n);

% update EZX in pml region

% FTONT
ezxbcf(1:iefbc,1:jebc)=caezxbcf(1:iefbc,1:jebc).*ezxbcf(1:iefbc,1:jebc)+...
cbezxbcf(1:iefbc,1:jebc).*(hybcf(2:ibfbc,1:jebc)-hybcf(1:iefbc,1:jebc));

% BACK
ezxbcb(1:iefbc,1:jebc)=caezxbcb(1:iefbc,1:jebc).*ezxbcb(1:iefbc,1:jebc)+...
cbezxbcb(1:iefbc,1:jebc).*(hybcb(2:ibfbc,1:jebc)-hybcb(1:iefbc,1:jebc));

% LEFT
ezxbcl(1:iebc-1,1:je)=caezxbcl(1:iebc-1,1:je).*ezxbcl(1:iebc-1,1:je)+...
cbezxbcl(1:iebc-1,1:je).*(hybcl(2:iebc,1:je)-hybcl(1:iebc-1,1:je));
ezxbcl(iebc,1:je)=caezxbcl(iebc,1:je).*ezxbcl(iebc,1:je)+...
cbezxbcl(iebc,1:je).*(hy(1,1:je)-hybcl(iebc,1:je));
% RIGHT
ezxbcr(2:iebc,1:je)=caezxbcr(2:iebc,1:je).*ezxbcr(2:iebc,1:je)+...
cbezxbcr(2:iebc,1:je).*(hybcr(3:ibbc,1:je)-hybcr(2:iebc,1:je));
ezxbcr(1,1:je)=caezxbcr(1,1:je).*ezxbcr(1,1:je)+...
cbezxbcr(1,1:je).*(hybcr(2,1:je)-hy(ib,1:je));

% update EZY in pml region

% FRONT
ezybcf(1:iefbc,1:jebc-1)=caezybcf(1:iefbc,1:jebc-1).*ezybcf(1:iefbc,1:jebc-1)-...
cbezybcf(1:iefbc,1:jebc-1).*(hxbcf(1:iefbc,2:jebc)-hxbcf(1:iefbc,1:jebc-1));
ezybcf(1:iebc,jebc)=caezybcf(1:iebc,jebc).*ezybcf(1:iebc,jebc)-...
cbezybcf(1:iebc,jebc).*(hxbcl(1:iebc,1)-hxbcf(1:iebc,jebc));
ezybcf(iebc+1:iebc+ie,jebc)=caezybcf(iebc+1:iebc+ie,jebc).*ezybcf(iebc+1:iebc+ie,jebc)-...
cbezybcf(iebc+1:ie+iebc,jebc).*(hx(1:ie,1)-hxbcf(iebc+1:iebc+ie,jebc));
ezybcf(iebc+ie+1:iefbc,jebc)=caezybcf(iebc+ie+1:iefbc,jebc).*ezybcf(ie+iebc+1:iefbc,jebc)-...
cbezybcf(ie+iebc+1:iefbc,jebc).*(hxbcr(1:iebc,1)-hxbcf(iebc+ie+1:iefbc,jebc));

% BACK
ezybcb(1:iefbc,2:jebc)=caezybcb(1:iefbc,2:jebc).*ezybcb(1:iefbc,2:jebc)-...
cbezybcb(1:iefbc,2:jebc).*(hxbcb(1:iefbc,3:jbbc)-hxbcb(1:iefbc,2:jebc));
ezybcb(1:iebc,1)=caezybcb(1:iebc,1).*ezybcb(1:iebc,1)-...
cbezybcb(1:iebc,1).*(hxbcb(1:iebc,2)-hxbcl(1:iebc,jb));
ezybcb(iebc+1:iebc+ie,1)=caezybcb(iebc+1:iebc+ie,1).*ezybcb(iebc+1:iebc+ie,1)-...
cbezybcb(iebc+1:iebc+ie,1).*(hxbcb(iebc+1:iebc+ie,2)-hx(1:ie,jb));
ezybcb(iebc+ie+1:iefbc,1)=caezybcb(iebc+ie+1:iefbc,1).*ezybcb(iebc+ie+1:iefbc,1)-...
cbezybcb(iebc+ie+1:iefbc,1).*(hxbcb(iebc+ie+1:iefbc,2)-hxbcr(1:iebc,jb));

% LEFT
ezybcl(1:iebc,1:je)=caezybcl(1:iebc,1:je).*ezybcl(1:iebc,1:je)-...
cbezybcl(1:iebc,1:je).*(hxbcl(1:iebc,2:jb)-hxbcl(1:iebc,1:je));

% RIGHT
ezybcr(1:iebc,1:je)=caezybcr(1:iebc,1:je).*ezybcr(1:iebc,1:je)-...
cbezybcr(1:iebc,1:je).*(hxbcr(1:iebc,2:jb)-hxbcr(1:iebc,1:je));

%***********************************************************************
% Update electric fields (HX and HY) in main grid
%***********************************************************************
hx(1:ie,2:je)=dahx(1:ie,2:je).*hx(1:ie,2:je)-dbhx(1:ie,2:je).*(ez(1:ie,2:je)-ez(1:ie,1:je-1));

hy(2:ie,1:je)=dahy(2:ie,1:je).*hy(2:ie,1:je)+dbhy(2:ie,1:je).*(ez(2:ie,1:je)-ez(1:ie-1,1:je));

%update HX in pml region

% FRONT
hxbcf(1:iefbc,2:jebc)=dahxbcf(1:iefbc,2:jebc).*hxbcf(1:iefbc,2:jebc)-...
dbhybcf(1:iefbc,2:jebc).*(ezxbcf(1:iefbc,2:jebc)+ezybcf(1:iefbc,2:jebc)-ezxbcf(1:iefbc,1:jebc-1)-ezybcf(1:iefbc,1:jebc-1));
hx(1:ie,1)=dahx(1:ie,1).*hx(1:ie,1)-dbhx(1:ie,1).*(ez(1:ie,1)-ezxbcf(iebc+1:iebc+ie,jebc)-ezybcf(iebc+1:iebc+ie,jebc));

% BACK
hxbcb(1:iefbc,2:jebc)=dahxbcb(1:iefbc,2:jebc).*hxbcb(1:iefbc,2:jebc)-...
dbhxbcb(1:iefbc,2:jebc).*(ezxbcb(1:iefbc,2:jebc)+ezybcb(1:iefbc,2:jebc)-ezxbcb(1:iefbc,1:jebc-1)-ezybcb(1:iefbc,1:jebc-1));
hx(1:ie,jb)=dahx(1:ie,jb).*hx(1:ie,jb)-dbhx(1:ie,jb).*(ezxbcb(iebc+1:iebc+ie,1)+ezybcb(iebc+1:iebc+ie,1)-ez(1:ie,je));

% LEFT
hxbcl(1:iebc,2:je)=dahxbcl(1:iebc,2:je).*hxbcl(1:iebc,2:je)-...
dbhxbcl(1:iebc,2:je).*(ezxbcl(1:iebc,2:je)+ezybcl(1:iebc,2:je)-ezxbcl(1:iebc,1:je-1)-ezybcl(1:iebc,1:je-1));
hxbcl(1:iebc,1)=dahxbcl(1:iebc,1).*hxbcl(1:iebc,1)-...
dbhxbcl(1:iebc,1).*(ezxbcl(1:iebc,1)+ezybcl(1:iebc,1)-ezxbcl(1:iebc,jebc)+ezybcl(1:iebc,jebc));
hxbcl(1:iebc,jb)=dahxbcl(1:iebc,jb).*hxbcl(1:iebc,jb)-...
dbhxbcl(1:iebc,jb).*(ezxbcb(1:iebc,1)+ezybcb(1:iebc,1)-ezxbcl(1:iebc,je)-ezxbcl(1:iebc,je));

% RIGHT
hxbcr(1:iebc,2:je)=dahxbcr(1:iebc,2:je).*hxbcr(1:iebc,2:je)-...
dbhxbcr(1:iebc,2:je).*(ezxbcr(1:iebc,2:je)+ezybcr(1:iebc,2:je)-ezxbcl(1:iebc,1:je-1)-ezybcl(1:iebc,1:je-1));
hxbcr(1:iebc,1)=dahxbcr(1:iebc,1).*hxbcr(1:iebc,1)-...
dbhxbcr(1:iebc,1).*(ezxbcr(1:iebc,1)+ezybcr(1:iebc,1)-ezxbcf(ie+iebc+1:iefbc,jebc)-ezybcf(ie+iebc+1:iefbc,jebc));
hxbcr(1:iebc,jb)=dahxbcr(1:iebc,jb).*hxbcr(1:iebc,1)-...
dbhxbcr(1:iebc,jb).*(ezxbcb(iebc+ie+1:iefbc,1)+ezybcb(iebc+ie+1:iefbc,1)-ezxbcr(1:iebc,je)-ezybcr(1:iebc,je));

%update HY in pml region

% FRONT
hybcf(2:iefbc,1:jebc)=dahybcf(2:iefbc,1:jebc).*hybcf(2:iefbc,1:jebc)+...
dbhybcf(2:iefbc,1:jebc).*(ezxbcf(2:iefbc,1:jebc)+ezybcf(2:iefbc,1:jebc)-ezxbcf(1:iefbc-1,1:jebc)-ezybcf(1:iefbc-1,1:jebc));

% BACK
hybcb(2:iefbc,1:jebc)=dahybcb(2:iefbc,1:jebc).*hybcb(2:iefbc,1:jebc)+...
dbhybcb(2:iefbc,1:jebc).*(ezxbcb(2:iefbc,1:jebc)+ezybcb(2:iefbc,1:jebc)-ezxbcb(1:iefbc-1,1:jebc)-ezybcb(1:iefbc-1,1:jebc));

% LEFT
hybcl(2:iebc,1:je)=dahybcl(2:iebc,1:je).*hybcl(2:iebc,1:je)+...
dbhybcl(2:iebc,1:je).*(ezxbcl(2:iebc,1:je)+ezybcl(2:iebc,1:je)-ezxbcl(1:iebc-1,1:je)-ezybcl(1:iebc-1,1:je));
hy(1,1:je)=dahy(1,1:je).*hy(1,1:je)+...
dbhy(1,1:je).*(ez(1,1:je)-ezxbcl(iebc,1:je)-ezybcl(iebc,1:je));

% RIGHT
hybcr(2:iebc,1:je)=dahybcr(2:iebc,1:je).*hybcr(2:i ..

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

附件:

tm模式pml中指数差分公式.doc (29 K) 下载次数:14

离线yuyeliuxu

UID ：21110

注册：2008-11-10
登录：2009-10-09
发帖：31
等级：仿真一级

1楼发表于: 2009-08-26 14:25:13

发散问题一般是程序中的方程编写出错，你可以看看是那个地方先开始发散的，问题就在那

共1条评分

gwzhao

技术分 +1

积极参与讨论+技术分论坛感谢您的参与

2009-08-26

离线gwzhao

方恨少

UID ：17098

注册：2008-08-24
登录：2019-01-09
发帖：1374
等级：荣誉管理员

2楼发表于: 2009-08-26 17:29:52

发散还是比较容易调试出来的，自己定下心来调调看吧。

逆流而上

离线juventus

UID ：20795

注册：2008-11-05
登录：2012-06-06
发帖：18
等级：仿真新人

3楼发表于: 2009-08-28 09:32:36

感谢yuyeliuxu和gwzhao两位提供的建议，小弟静心调试得到了想要的结果。因是新人小弟到处怀疑可能出错的地方浪费的不少时间，回过头来才明白两位提供的建议真乃真知灼见，万分感谢。小弟也把修改后正确的程序贴上，希望也能对其他新人有所帮助。
该程序为2D情况下pml吸收边界自由空间电磁波的传输，中间虽然有一段写到加入圆形柱体，但所设参数和自由空间电磁参数一致，故而可以认为没有添加散射圆柱。
clear

%***********************************************************************
%     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

freq=5.0e+9;                %center frequency of source excitation
lambda=cc/freq;             %center wavelength of source excitation
omega=2.0*pi*freq;

%***********************************************************************
%     Grid parameters
%***********************************************************************

ie=100;           %number of grid cells in x-direction
je=100;            %number of grid cells in y-direction

ib=ie+1;
jb=je+1;

is=ie/2;            %location of z-directed hard source
js=je/2;          %location of z-directed hard source

dx=3.0e-3;        %space increment of square lattice
dt=dx/(2.0*cc);   %time step

nmax=3000;         %total number of time steps

iebc=8;           %thickness of left and right PML region
jebc=8;           %thickness of front and back PML region
rmax=0.00001;
orderbc=2;
ibbc=iebc+1;
jbbc=jebc+1;
iefbc=ie+2*iebc;
jefbc=je+2*jebc;
ibfbc=iefbc+1;
jbfbc=jefbc+1;

%***********************************************************************
%     Material parameters
%***********************************************************************

media=2;

eps=[1.0 1.0];
sig=[0.0 0.0];
mur=[1.0 1.0];
sim=[0.0 0.0];

%***********************************************************************
%     Wave excitation
%***********************************************************************

source=zeros(1,nmax);
for n=1:nmax
  source(n)=90*sin(omega*n*dt);
end

%***********************************************************************
%     Field arrays
%***********************************************************************
ez=zeros(ie,je);
hx=zeros(ie,jb); %fields in main grid
hy=zeros(ib,je);

ezxbcf=zeros(iefbc,jebc);
ezybcf=zeros(iefbc,jebc);
hxbcf=zeros(iefbc,jebc);  %fields in front pml
hybcf=zeros(ibfbc,jebc);

ezxbcb=zeros(iefbc,jebc);
ezybcb=zeros(iefbc,jebc);
hxbcb=zeros(iefbc,jbbc); %fields in back pml
hybcb=zeros(ibfbc,jebc);

ezxbcl=zeros(iebc,je);
ezybcl=zeros(iebc,je);
hxbcl=zeros(iebc,jb);   %fields in left pml
hybcl=zeros(iebc,je);

ezxbcr=zeros(iebc,je);
ezybcr=zeros(iebc,je);
hxbcr=zeros(iebc,jb);  %fields in right pml
hybcr=zeros(ibbc,je);

%***********************************************************************
%     Updating coefficients
%***********************************************************************

for i=1:media
  eaf  =dt*sig(i)/(2.0*epsz*eps(i));
  ca(i)=(1.0-eaf)/(1.0+eaf);
  cb(i)=dt/epsz/eps(i)/dx/(1.0+eaf);
  haf  =dt*sim(i)/(2.0*muz*mur(i));
  da(i)=(1.0-haf)/(1.0+haf);
  db(i)=dt/muz/mur(i)/dx/(1.0+haf);
end

%***********************************************************************
%     Geometry specification (main grid)
%***********************************************************************

%     Initialize entire main grid to free space
caez(1:ie,1:je)=ca(1);
cbez(1:ie,1:je)=cb(1);

dahx(1:ie,1:jb)=da(1);
dbhx(1:ie,1:jb)=db(1);

dahy(1:ib,1:je)=da(1);
dbhy(1:ib,1:je)=db(1);

%     Add metal cylinder

diam=20;          % diameter of cylinder: 6 cm
rad=diam/2.0;     % radius of cylinder: 3 cm
icenter=4*ie/5;   % i-coordinate of cylinder's center
jcenter=je/2;     % j-coordinate of cylinder's center

for i=1:ie
for j=1:je
  dist2=(i-icenter)^2 + (j-jcenter)^2;
  if dist2 <= rad^2
     caez(i,j)=ca(2);
     cbez(i,j)=cb(2);
  end
end
end

%***********************************************************************
%     Fill the PML regions
%***********************************************************************

delbc=iebc*dx;
sigmam=-log(rmax/100.0)*epsz*cc*(orderbc+1)/(2*delbc);
bcfactor=eps(1)*sigmam/(dx*(delbc^orderbc)*(orderbc+1));

%     FRONT region
dahxbcf(1:iefbc,1)=1.0;
dbhxbcf(1:iefbc,1)=0.0;
for j=2:jebc
  y1=(jebc-j+1.5)*dx;
  y2=(jebc-j+0.5)*dx;
  sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
  sigmays=sigmay*(muz/(epsz*eps(1)));
  da1=exp(-sigmays*dt/muz);
  db1=(1-da1)/(sigmays*dx);
  dahxbcf(1:iefbc,j)=da1;
  dbhxbcf(1:iefbc,j)=db1;
end
sigmay=bcfactor*(0.5*dx)^(orderbc+1);
sigmays=sigmay*(muz/(epsz*eps(1)));
da1=exp(-sigmays*dt/muz);
db1=(1-da1)/(sigmays*dx);
dahx(1:ie,1)=da1;
dbhx(1:ie,1)=db1;
dahxbcl(1:iebc,1)=da1;
dbhxbcl(1:iebc,1)=db1;
dahxbcr(1:iebc,1)=da1;
dbhxbcr(1:iebc,1)=db1;

for j=1:iebc
    y1=(jebc-j+1)*dx;
    y2=(jebc-j)*dx;
    sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
    ca1=exp(-sigmay*dt/(epsz*eps(1)));
    cb1=(1.0-ca1)/(sigmay*dx);
    caezybcf(1:iefbc,j)=ca1;
    cbezybcf(1:iefbc,j)=cb1;
    caezxbcf(1:iefbc,j)=ca(1);
    cbezxbcf(1:iefbc,j)=cb(1);
    dahybcf(1:ibfbc,j)=da(1);
    dbhybcf(1:ibfbc,j)=db(1);
end

%      BACK region
dahxbcb(1:iefbc,jbbc)=1.0;
dbhxbcb(1:iefbc,jbbc)=0.0;
for j=2:jebc
    y1=(j-0.5)*dx;
    y2=(j-1.5)*dx;
    sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
    sigmays=sigmay*(muz/(epsz*eps(1)));
    da1=exp(-sigmays*dt/muz);
    db1=(1-da1)/(sigmays*dx);
    dahxbcb(1:iefbc,j)=da1;
    dbhxbcb(1:iefbc,j)=db1;
end
sigmay=bcfactor*(0.5*dx)^(orderbc+1);
sigmays=sigmay*(muz/(epsz*eps(1)));
da1=exp(-sigmays*dt/muz);
db1=(1-da1)/(sigmays*dx);
dahx(1:ie,jb)=da1;
dbhx(1:ie,jb)=db1;
dahxbcl(1:iebc,jb)=da1;
dbhxbcl(1:iebc,jb)=db1;
dahxbcr(1:iebc,jb)=da1;
dbhxbcr(1:iebc,jb)=db1;

for j=1:jebc
    y1=j*dx;
    y2=(j-1)*dx;
    sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
    ca1=exp(-sigmay*dt/(epsz*eps(1)));
    cb1=(1-ca1)/(sigmay*dx);
    caezybcb(1:iefbc,j)=ca1;
    cbezybcb(1:iefbc,j)=cb1;
    caezxbcb(1:iefbc,j)=ca(1);
    cbezxbcb(1:iefbc,j)=cb(1);
    dahybcb(1:ibfbc,j)=da(1);
    dbhybcb(1:ibfbc,j)=db(1);
end

%     LEFT region
dahybcl(1,1:je)=1.0;
dbhybcl(1,1:je)=0.0;
for i=2:iebc
    x1=(iebc-i+1.5)*dx;
    x2=(iebc-i+0.5)*dx;
    sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
    sigmaxs=sigmax*(muz/(epsz*eps(1)));
    da1=exp(-sigmaxs*dt/muz);
    db1=(1-da1)/(sigmaxs*dx);
    dahybcl(i,1:je)=da1;
    dbhybcl(i,1:je)=db1;
    dahybcf(i,1:jebc)=da1;
    dbhybcf(i,1:jebc)=db1;
    dahybcb(i,1:jebc)=da1;
    dbhybcb(i,1:jebc)=db1;
end
sigmax=bcfactor*(0.5*dx)^(orderbc+1);
sigmaxs=sigmax*(muz/(epsz*eps(1)));
da1=exp(-sigmaxs*dt/muz);
db1=(1-da1)/(sigmaxs*dx);
dahy(1,1:je)=da1;
dbhy(1,1:je)=db1;
dahybcf(ibbc,1:jebc)=da1;
dbhybcf(ibbc,1:jebc)=db1;
dahybcb(ibbc,1:jebc)=da1;
dbhybcb(ibbc,1:jebc)=db1;

for i=1:iebc
    x1=(iebc-i+1)*dx;
    x2=(iebc-i)*dx;
    sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
    ca1=exp(-sigmax*dt/(epsz*eps(1)));
    cb1=(1-ca1)/(sigmax*dx);
    caezxbcl(i,1:je)=ca1;
    cbezxbcl(i,1:je)=cb1;
    caezxbcf(i,1:jebc)=ca1;
    cbezxbcf(i,1:jebc)=cb1;
    caezxbcb(i,1:jebc)=ca1;
    cbezxbcb(i,1:jebc)=cb1;

    caezybcl(i,1:je)=ca(1);
    cbezybcl(i,1:je)=cb(1);
    dahxbcl(i,2:je)=da(1);
    dbhxbcl(i,2:je)=db(1);
end

%        RIGHT region
dahybcr(ibbc,1:je)=1.0;
dbhybcr(ibbc,1:je)=0.0;
for i=2:iebc
    x1=(i-0.5)*dx;
    x2=(i-1.5)*dx;
    sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
    sigmaxs=sigmax*(muz/(epsz*eps(1)));
    da1=exp(-sigmaxs*dt/muz);
    db1=(1-da1)/(sigmaxs*dx);
    dahybcr(i,1:je)=da1;
    dbhybcr(i,1:je)=db1;
    dahybcf(iebc+ie+i,1:jebc)=da1;
    dbhybcf(iebc+ie+i,1:jebc)=db1;
    dahybcb(iebc+ie+i,1:jebc)=da1;
    dbhybcb(iebc+ie+i,1:jebc)=db1;
end
sigmax=bcfactor*(0.5*dx)^(orderbc+1);
sigmaxs=sigmax*(muz/(epsz*eps(1)));
da1=exp(-sigmaxs*dt/muz);
db1=(1-da1)/(sigmaxs*dx);
dahy(ib,1:je)=da1;
dbhy(ib,1:je)=db1;
dahybcf(iebc+ie+1,1:jebc)=da1;
dbhybcf(iebc+ie+1,1:jebc)=db1;
dahybcb(iebc+ie+1,1:jebc)=da1;
dbhybcb(iebc+ie+1,1:jebc)=db1;

for i=1:iebc
    x1=i*dx;
    x2=(i-1)*dx;
    sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
    ca1=exp(-sigmax*dt/(epsz*eps(1)));
    cb1=(1-ca1)/(sigmax*dx);
    caezxbcr(i,1:je)=ca1;
    cbezxbcr(i,1:je)=cb1;
    caezxbcf(iebc+ie+i,1:jebc)=ca1;
    cbezxbcf(iebc+ie+i,1:jebc)=cb1;
    caezxbcb(iebc+ie+i,1:jebc)=ca1;
    cbezxbcb(iebc+ie+i,1:jebc)=cb1;

    caezybcr(i,1:je)=ca(1);
    cbezybcr(i,1:je)=cb(1);
    dahxbcr(i,2:je)=da(1);
    dbhxbcr(i,2:je)=db(1);
end
%***********************************************************************
%     Movie initialization
%***********************************************************************

subplot(1,1,1),pcolor(ez');
shading flat;
caxis([-80.0 80.0]);
axis([1 ie 1 je]);
colorbar;
axis image;
axis off;
title(['Ez at time step=0']);

rect=get(gcf,'Position');
rect(1:2)=[0 0];

M=moviein(nmax/4,gcf,rect);
%***********************************************************************
%     BEGIN TIME-STEPPING LOOP
%***********************************************************************
for n=1:nmax

%***********************************************************************
%     Update electric fields (EZ) in main grid
%***********************************************************************
ez(1:ie,1:je)=caez(1:ie,1:je).*ez(1:ie,1:je)+cbez(1:ie,1:je).*(hy(2:ib,1:je)-hy(1:ie,1:je)-...
              (hx(1:ie,2:jb)-hx(1:ie,1:je)));
ez(is,js)=source(n);

%  update EZX in pml region

%   FTONT
ezxbcf(1:iefbc,1:jebc)=caezxbcf(1:iefbc,1:jebc).*ezxbcf(1:iefbc,1:jebc)+...
                       cbezxbcf(1:iefbc,1:jebc).*(hybcf(2:ibfbc,1:jebc)-hybcf(1:iefbc,1:jebc));

%  BACK
ezxbcb(1:iefbc,1:jebc)=caezxbcb(1:iefbc,1:jebc).*ezxbcb(1:iefbc,1:jebc)+...
                       cbezxbcb(1:iefbc,1:jebc).*(hybcb(2:ibfbc,1:jebc)-hybcb(1:iefbc,1:jebc));

%  LEFT
ezxbcl(1:iebc-1,1:je)=caezxbcl(1:iebc-1,1:je).*ezxbcl(1:iebc-1,1:je)+...
                       cbezxbcl(1:iebc-1,1:je).*(hybcl(2:iebc,1:je)-hybcl(1:iebc-1,1:je));
ezxbcl(iebc,1:je)=caezxbcl(iebc,1:je).*ezxbcl(iebc,1:je)+...
                   cbezxbcl(iebc,1:je).*(hy(1,1:je)-hybcl(iebc,1:je));
%  RIGHT
ezxbcr(2:iebc,1:je)=caezxbcr(2:iebc,1:je).*ezxbcr(2:iebc,1:je)+...
                     cbezxbcr(2:iebc,1:je).*(hybcr(3:ibbc,1:je)-hybcr(2:iebc,1:je));
ezxbcr(1,1:je)=caezxbcr(1,1:je).*ezxbcr(1,1:je)+...
                cbezxbcr(1,1:je).*(hybcr(2,1:je)-hy(ib,1:je));

% update EZY in pml region

% FRONT
ezybcf(1:iefbc,1:jebc-1)=caezybcf(1:iefbc,1:jebc-1).*ezybcf(1:iefbc,1:jebc-1)-...
                          cbezybcf(1:iefbc,1:jebc-1).*(hxbcf(1:iefbc,2:jebc)-hxbcf(1:iefbc,1:jebc-1));
ezybcf(1:iebc,jebc)=caezybcf(1:iebc,jebc).*ezybcf(1:iebc,jebc)-...
                     cbezybcf(1:iebc,jebc).*(hxbcl(1:iebc,1)-hxbcf(1:iebc,jebc));
ezybcf(iebc+1:iebc+ie,jebc)=caezybcf(iebc+1:iebc+ie,jebc).*ezybcf(iebc+1:iebc+ie,jebc)-...
                     cbezybcf(iebc+1:ie+iebc,jebc).*(hx(1:ie,1)-hxbcf(iebc+1:iebc+ie,jebc));
ezybcf(iebc+ie+1:iefbc,jebc)=caezybcf(iebc+ie+1:iefbc,jebc).*ezybcf(ie+iebc+1:iefbc,jebc)-...
                     cbezybcf(ie+iebc+1:iefbc,jebc).*(hxbcr(1:iebc,1)-hxbcf(iebc+ie+1:iefbc,jebc));

% BACK
ezybcb(1:iefbc,2:jebc)=caezybcb(1:iefbc,2:jebc).*ezybcb(1:iefbc,2:jebc)-...
                        cbezybcb(1:iefbc,2:jebc).*(hxbcb(1:iefbc,3:jbbc)-hxbcb(1:iefbc,2:jebc));
ezybcb(1:iebc,1)=caezybcb(1:iebc,1).*ezybcb(1:iebc,1)-...
                  cbezybcb(1:iebc,1).*(hxbcb(1:iebc,2)-hxbcl(1:iebc,jb));
ezybcb(iebc+1:iebc+ie,1)=caezybcb(iebc+1:iebc+ie,1).*ezybcb(iebc+1:iebc+ie,1)-...
                          cbezybcb(iebc+1:iebc+ie,1).*(hxbcb(iebc+1:iebc+ie,2)-hx(1:ie,jb));
ezybcb(iebc+ie+1:iefbc,1)=caezybcb(iebc+ie+1:iefbc,1).*ezybcb(iebc+ie+1:iefbc,1)-...
                           cbezybcb(iebc+ie+1:iefbc,1).*(hxbcb(iebc+ie+1:iefbc,2)-hxbcr(1:iebc,jb));

%  LEFT
ezybcl(1:iebc,1:je)=caezybcl(1:iebc,1:je).*ezybcl(1:iebc,1:je)-...
                     cbezybcl(1:iebc,1:je).*(hxbcl(1:iebc,2:jb)-hxbcl(1:iebc,1:je));

%  RIGHT
ezybcr(1:iebc,1:je)=caezybcr(1:iebc,1:je).*ezybcr(1:iebc,1:je)-...
                     cbezybcr(1:iebc,1:je).*(hxbcr(1:iebc,2:jb)-hxbcr(1:iebc,1:je));

%***********************************************************************
%     Update electric fields (HX and HY) in main grid
%***********************************************************************
hx(1:ie,2:je)=dahx(1:ie,2:je).*hx(1:ie,2:je)-dbhx(1:ie,2:je).*(ez(1:ie,2:je)-ez(1:ie,1:je-1));

hy(2:ie,1:je)=dahy(2:ie,1:je).*hy(2:ie,1:je)+dbhy(2:ie,1:je).*(ez(2:ie,1:je)-ez(1:ie-1,1:je));

%update HX in pml region

%   FRONT
hxbcf(1:iefbc,2:jebc)=dahxbcf(1:iefbc,2:jebc).*hxbcf(1:iefbc,2:jebc)-...
                      dbhxbcf(1:iefbc,2:jebc).*(ezxbcf(1:iefbc,2:jebc)+ezybcf(1:iefbc,2:jebc)-ezxbcf(1:iefbc,1:jebc-1)-ezybcf(1:iefbc,1:jebc-1));
hx(1:ie,1)=dahx(1:ie,1).*hx(1:ie,1)-dbhx(1:ie,1).*(ez(1:ie,1)-ezxbcf(iebc+1:iebc+ie,jebc)-ezybcf(iebc+1:iebc+ie,jebc));

%   BACK
hxbcb(1:iefbc,2:jebc)=dahxbcb(1:iefbc,2:jebc).*hxbcb(1:iefbc,2:jebc)-...
                      dbhxbcb(1:iefbc,2:jebc).*(ezxbcb(1:iefbc,2:jebc)+ezybcb(1:iefbc,2:jebc)-ezxbcb(1:iefbc,1:jebc-1)-ezybcb(1:iefbc,1:jebc-1));
hx(1:ie,jb)=dahx(1:ie,jb).*hx(1:ie,jb)-dbhx(1:ie,jb).*(ezxbcb(iebc+1:iebc+ie,1)+ezybcb(iebc+1:iebc+ie,1)-ez(1:ie,je));

%   LEFT
hxbcl(1:iebc,2:je)=dahxbcl(1:iebc,2:je).*hxbcl(1:iebc,2:je)-...
                   dbhxbcl(1:iebc,2:je).*(ezxbcl(1:iebc,2:je)+ezybcl(1:iebc,2:je)-ezxbcl(1:iebc,1:je-1)-ezybcl(1:iebc,1:je-1));
hxbcl(1:iebc,1)=dahxbcl(1:iebc,1).*hxbcl(1:iebc,1)-...
                dbhxbcl(1:iebc,1).*(ezxbcl(1:iebc,1)+ezybcl(1:iebc,1)-ezxbcf(1:iebc,jebc)-ezybcf(1:iebc,jebc));
hxbcl(1:iebc,jb)=dahxbcl(1:iebc,jb).*hxbcl(1:iebc,jb)-...
                 dbhxbcl(1:iebc,jb).*(ezxbcb(1:iebc,1)+ezybcb(1:iebc,1)-ezxbcl(1:iebc,je)-ezybcl(1:iebc,je));

%   RIGHT
hxbcr(1:iebc,2:je)=dahxbcr(1:iebc,2:je).*hxbcr(1:iebc,2:je)-...
                   dbhxbcr(1:iebc,2:je).*(ezxbcr(1:iebc,2:je)+ezybcr(1:iebc,2:je)-ezxbcr(1:iebc,1:je-1)-ezybcr(1:iebc,1:je-1));
hxbcr(1:iebc,1)=dahxbcr(1:iebc,1).*hxbcr(1:iebc,1)-...
                dbhxbcr(1:iebc,1).*(ezxbcr(1:iebc,1)+ezybcr(1:iebc,1)-ezxbcf(ie+iebc+1:iefbc,jebc)-ezybcf(ie+iebc+1:iefbc,jebc));
hxbcr(1:iebc,jb)=dahxbcr(1:iebc,jb).*hxbcr(1:iebc,jb)-...
                 dbhxbcr(1:iebc,jb).*(ezxbcb(iebc+ie+1:iefbc,1)+ezybcb(iebc+ie+1:iefbc,1)-ezxbcr(1:iebc,je)-ezybcr(1:iebc,je));


%update HY in pml region

% FRONT
hybcf(2:iefbc,1:jebc)=dahybcf(2:iefbc,1:jebc).*hybcf(2:iefbc,1:jebc)+...
                      dbhybcf(2:iefbc,1:jebc).*(ezxbcf(2:iefbc,1:jebc)+ezybcf(2:iefbc,1:jebc)-ezxbcf(1:iefbc-1,1:jebc)-ezybcf(1:iefbc-1,1:jebc));

% BACK
hybcb(2:iefbc,1:jebc)=dahybcb(2:iefbc,1:jebc).*hybcb(2:iefbc,1:jebc)+...
                      dbhybcb(2:iefbc,1:jebc).*(ezxbcb(2:iefbc,1:jebc)+ezybcb(2:iefbc,1:jebc)-ezxbcb(1:iefbc-1,1:jebc)-ezybcb(1:iefbc-1,1:jebc));

% LEFT
hybcl(2:iebc,1:je)=dahybcl(2:iebc,1:je).*hybcl(2:iebc,1:je)+...
                   dbhybcl(2:iebc,1:je).*(ezxbcl(2:iebc,1:je)+ezybcl(2:iebc,1:je)-ezxbcl(1:iebc-1,1:je)-ezybcl(1:iebc-1,1:je));
hy(1,1:je)=dahy(1,1:je).*hy(1,1:je)+...
           dbhy(1,1:je).*(ez(1,1:je)-ezxbcl(iebc,1:je)-ezybcl(iebc,1:je));

% RIGHT
hybcr(2:iebc,1:je)=dahybcr(2:iebc,1:je).*hybcr(2:iebc,1:je)+...
                   dbhybcr(2:iebc,1:je).*(ezxbcr(2:iebc,1:je)+ezybcr(2:iebc,1:je)-ezxbcr(1:iebc-1,1:je)-ezybcr(1:iebc-1,1:je));
hy(ib,1:je)=dahy(ib,1:je).*hy(ib,1:je)+...
            dbhy(ib,1:je).*(ezxbcr(1,1:je)+ezybcr(1,1:je)-ez(ie,1:je));

%***********************************************************************
%     Visualize fields
%***********************************************************************

if mod(n,4)==0;

timestep=int2str(n);

subplot(1,1,1),pcolor(ez');
shading flat;
caxis([-80.0 80.0]);
axis([1 ie 1 je]);
colorbar;
axis image;
axis off;
title(['Ez at time step = ',timestep]);

nn=n/4;
M(:,nn)=getframe(gcf,rect);

end;

%***********************************************************************
%     END TIME-STEPPING LOOP
%***********************************************************************

end

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

离线alang

UID ：41410

注册：2009-09-11
登录：2010-03-02
发帖：77
等级：仿真一级

4楼发表于: 2010-02-07 11:50:14

xiexielouzhufenxianga

离线runner

UID ：3922

注册：2007-07-17
登录：2011-06-08
发帖：67
等级：仿真一级

5楼发表于: 2010-05-26 21:48:53

太长了，没信心看完

离线willi111

UID ：44212

注册：2009-10-17
登录：2015-04-12
发帖：74
等级：仿真一级

6楼发表于: 2012-06-29 20:08:45

真心感谢楼主，之前试ADE-CPML试了很长时间也没结果，一直有反射，很是没头绪，楼主的分裂场PML受教了，是一个很好的参考。

离线dudu86

坚持，是最好的品质

UID ：94035

注册：2012-05-15
登录：2013-10-24
发帖：83
等级：仿真一级

7楼发表于: 2013-05-27 10:22:11

我想请教楼主怎么解决的发散问题呢，是PML的设置还是电、磁场分量的迭代公式呢？

坚持，是最好的品质！

发帖回复