1302阅读
1回复

[求助]帮忙看看我的编程思路有没有问题(附程序)

UID ：21110

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

0楼发表于: 2009-07-10 10:02:08

得到的结果pml不吸收
tic;
clear; clc;
N_iteration=1;
%***********************************************************************
%     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=5e+9;                  %frequency of source excitation
Tc=0.2e-9                %Time shift, i.e. time for the center position of the pulse
Tau=1/freq;                 %impulse width parameter, e.g.: 0.2ns pulse width here
                            % corresponds frequency of 1/Tau Hz,  i.e.
                            % 0.2ns---5GHz,  0.5ns---2GHz
lambda=cc/freq;             %center wavelength of source excitation
omega=2.0*pi*freq;

M_factor=1;%1       % refine the mesh to get higher quality M_factor times meshed on 1e-3

nmax=500*M_factor;  %total number of calculating time steps   nmax=2200*M_factor;
nstep=10;           % =250 every nstep give a display

%***********************************************************************
%     Grid parameters
%***********************************************************************
ie=100*M_factor;   %100           %number of grid cells in x-direction
je=100*M_factor;  %100          %number of grid cells in y-direction

ib=ie+1;
jb=je+1;
imp_x=50*M_factor;   %imp_x=5*M_factor; %imp_x=10*M_factor;
imp_y=50*M_factor;         %(imp_x ,imp_y) is the location of the point source impulse

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%   fine meshed      %%%%
dx=1e-3/M_factor   ;%  dx=0.25e-3;    %dx=1e-3  %dx=0.25e-3;    %dx=1e-3/M_factor;     %space increment of square lattice
dt=dx/(2.0*cc);   %time step  dt=dx/(2.0*cc);

iebc=8*M_factor;           %thickness of left and right PML region
jebc=8*M_factor;           %thickness of front and back PML region
rmax=0.000001;   %reflection factor ?x?? concerning PML issue
orderbc=2;
ibbc=iebc+1;
jbbc=jebc+1;
iefbc=ie+2*iebc;
jefbc=je+2*jebc;
ibfbc=iefbc+1;
jbfbc=jefbc+1;

%***********************************************************************
%     Material parameters
%***********************************************************************
epss=10;
eps=7;
tau0=7e-12;
sigdisp=(epss-eps)*epsz/tau0;
sig=0.15+sigdisp;
mur=1.0;
sim=0;
%***********************************************************************
%     Wave excitation
%***********************************************************************
source=zeros(1,nmax);
%     for n=1:nmax
%         t(n)=dt*n;
%         source(n)=-sqrt(exp(1))*(2*pi/Tau)*(t(n)-Tc)*exp(-(1/2)*(2*pi*(t(n)-Tc)/Tau)^2);  % normalized
%     end
    for  n=1:nmax           % n=n1:n2
         t(n)=dt*n;
         source(n)=-sin(omega*(t(n)-Tc));   %  %Eq.(6) of Monocycle shapes paper
    end
%***********************************************************************
%     Field arrays
%***********************************************************************

ex=zeros(ie,jb);           %fields in main grid
px=zeros(ie,jb);
ey=zeros(ib,je);
py=zeros(ib,je);
hz=zeros(ie,je);

exbcf=zeros(iefbc,jebc);   %fields in front PML region
pxbcf=zeros(iefbc,jebc);
eybcf=zeros(ibfbc,jebc);
pybcf=zeros(ibfbc,jebc);
hzxbcf=zeros(iefbc,jebc);
hzybcf=zeros(iefbc,jebc);

exbcb=zeros(iefbc,jbbc);   %fields in back PML region
pxbcb=zeros(iefbc,jbbc);
eybcb=zeros(ibfbc,jebc);
pybcb=zeros(ibfbc,jebc);
hzxbcb=zeros(iefbc,jebc);
hzybcb=zeros(iefbc,jebc);

exbcl=zeros(iebc,jb);      %fields in left PML region
pxbcl=zeros(iebc,jb);
eybcl=zeros(iebc,je);
pybcl=zeros(iebc,je);
hzxbcl=zeros(iebc,je);
hzybcl=zeros(iebc,je);

exbcr=zeros(iebc,jb);      %fields in right PML region
pxbcr=zeros(iebc,jb);
eybcr=zeros(ibbc,je);
pybcr=zeros(ibbc,je);
hzxbcr=zeros(iebc,je);
hzybcr=zeros(iebc,je);

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


  haf=dt*sim/(2.0*muz*mur);
  da=(1.0-haf)/(1.0+haf);
  db=dt/muz/mur/dx/(1.0+haf);
dapm=(2*tau0-dt)/(2*tau0+dt);
dbpm=2*tau0*dt/(2*tau0+dt)*sigdisp;

%***********************************************************************
%     Geometry specification (main grid)
%***********************************************************************
  haf=dt*sim/(2.0*muz*mur);
dahz(1:ie,1:je)=(1.0-haf)/(1.0+haf);
dbhz(1:ie,1:je)=dt/muz/mur/dx/(1.0+haf);
for i=1:ie
    for j=1:jb

eaf=dt*sig/(2.0*epsz*eps);
caex(i,j)=(1.0-eaf)/(1.0+eaf);                % parameters refer to page 6 from
cpex(i,j)=dt/epsz/eps/tau0/(1.0+eaf);
cbex(i,j)=dt/epsz/eps/dx/(1.0+eaf);        % p.6 from the book of Yun Wenhua
dapx(i,j)=(2*tau0-dt)/(2*tau0+dt);
dbpx(i,j)=2*tau0*dt/(2*tau0+dt)*sigdisp;

    end
end

for i=1: ib
    for j=1: je
    eaf=dt*sig/(2.0*epsz*eps);
caey(i,j)=(1.0-eaf)/(1.0+eaf);                % parameters refer to page 6 from
cpey(i,j)=dt/epsz/eps/tau0/(1.0+eaf);
cbey(i,j)=dt/epsz/eps/dx/(1.0+eaf);        % p.6 from the book of Yun Wenhua
dapy(i,j)=(2*tau0-dt)/(2*tau0+dt);
dbpy(i,j)=2*tau0*dt/(2*tau0+dt)*sigdisp;

    end
end
%***********************************************************************
%     Fill the PML regions
%***********************************************************************
delbc=iebc*dx;
sigmam=-log(rmax/100.0)*epsz*cc*(orderbc+1)/(2*delbc);

  eaf=dt*sig/(epsz*eps);
  capml=exp(-eaf);
  cbpml=(1.0-capml)*dt/(eaf*epsz*eps*dx);
  cppml=(1.0-capml)*dt/(eaf*epsz*eps*tau0);      % p.6 from the book of Yun Wenhua
dap=exp(-dt/tau0);
dbp=(1-dap)*sigdisp*tau0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     FRONT region
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

caexbcf(1:iefbc,1)=1.0;        %0.0   %1.0;
cbexbcf(1:iefbc,1)=0.0;
cpexbcf(1:iefbc,1)=0.0;
dapxbcf(1:iefbc,1)=1.0;
dbpxbcf(1:iefbc,1)=0.0;
bcfactor=eps*sigmam/(dx*(delbc^orderbc)*(orderbc+1));
for j=2:jebc
  y1=(jebc-j+1.5)*dx;
  y2=(jebc-j+0.5)*dx;
  sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
   eaf1=dt*(sigmay+sigdisp)/(epsz*eps);
  ca1=exp(-eaf1);
  cb1=(1.0-ca1)/((sigmay+sigdisp)*dx);
  cp1=(1.0-ca1)/((sigmay+sigdisp)*tau0);
  caexbcf(1:iefbc,j)=ca1;
  cbexbcf(1:iefbc,j)=cb1;
  cpexbcf(1:iefbc,j)=cp1;
  dapxbcf(1:iefbc,j)=dap;
  dbpxbcf(1:iefbc,j)=dbp;
end
sigmay = bcfactor*(0.5*dx)^(orderbc+1);
  eaf1=dt*(sigmay+sigdisp)/(epsz*eps);
  ca1=exp(-eaf1);
  cb1=(1.0-ca1)/((sigmay+sigdisp)*dx);
  cp1=(1.0-ca1)/((sigmay+sigdisp)*tau0);

caex(1:ie,1)=ca1;
cbex(1:ie,1)=cb1;
cpex(1:ie,1)=cp1;
dapx(1:ie,1)=dap;
dbpx(1:ie,1)=dbp;
caexbcl(1:iebc,1)=ca1;
cbexbcl(1:iebc,1)=cb1;
cpexbcl(1:iebc,1)=cp1;
dapxbcl(1:iebc,1)=dap;
dbpxbcl(1:iebc,1)=dbp;
caexbcr(1:iebc,1)=ca1;
cbexbcr(1:iebc,1)=cb1;
cpexbcr(1:iebc,1)=cp1;
dapxbcr(1:iebc,1)=dap;
dbpxbcr(1:iebc,1)=dbp;
for j=1:jebc
  y1=(jebc-j+1)*dx;
  y2=(jebc-j)*dx;
  sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
  sigmays=sigmay*(muz/(epsz*eps));
  da1=exp(-sigmays*dt/muz);
  db1=(1-da1)/(sigmays*dx);
  dahzybcf(1:iefbc,j)=da1;
  dbhzybcf(1:iefbc,j)=db1;
  caeybcf(1:ibfbc,j)=capml;
  cbeybcf(1:ibfbc,j)=cbpml;
  cpeybcf(1:ibfbc,j)=cppml;

  dahzxbcf(1:iefbc,j)=da;
  dbhzxbcf(1:iefbc,j)=db;

   dapybcf(1:ibfbc,j)=dap;
  dbpybcf(1:ibfbc,j)=dbp;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     BACK region
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
caexbcb(1:iefbc,jbbc)=1.0;   % 0.0 %1.0;
cbexbcb(1:iefbc,jbbc)=0.0;
cpexbcb(1:iefbc,jbbc)=0.0;
   dapxbcb(1:iefbc,1)=1.0;
  dbpxbcb(1:iefbc,1)=0.0;
bcfactor=eps*sigmam/(dx*(delbc^orderbc)*(orderbc+1));
for j=2:jebc
  y1=(j-0.5)*dx;
  y2=(j-1.5)*dx;
  sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
eaf1=dt*(sigmay+sigdisp)/(epsz*eps);
  ca1=exp(-eaf1);
  cb1=(1.0-ca1)/((sigmay+sigdisp)*dx);
  cp1=(1.0-ca1)/((sigmay+sigdisp)*tau0);
  caexbcb(1:iefbc,j)=ca1;
  cbexbcb(1:iefbc,j)=cb1;
  cpexbcb(1:iefbc,j)=cp1;
  dapxbcb(1:iefbc,j)=dap;
  dbpxbcb(1:iefbc,j)=dbp;
end
sigmay = bcfactor*(0.5*dx)^(orderbc+1);

   eaf1=dt*(sigmay+sigdisp)/(epsz*eps);
  ca1=exp(-eaf1);
  cb1=(1.0-ca1)/((sigmay+sigdisp)*dx);
  cp1=(1.0-ca1)/((sigmay+sigdisp)*tau0);
caex(1:ie,jb)=ca1;
cbex(1:ie,jb)=cb1;
cpex(1:ie,jb)=cp1;
dapx(1:ie,jb)=dap;
  dbpx(1:ie,jb)=dbp;
caexbcl(1:iebc,jb)=ca1;
cbexbcl(1:iebc,jb)=cb1;
cpexbcl(1:iebc,jb)=cp1;
dapxbcl(1:iebc,jb)=dap;
  dbpxbcl(1:iebc,jb)=dbp;
caexbcr(1:iebc,jb)=ca1;
cbexbcr(1:iebc,jb)=cb1;
cpexbcr(1:iebc,jb)=cp1;
dapxbcr(1:iebc,jb)=dap;
dbpxbcr(1:iebc,jb)=dbp;
for j=1:jebc
  y1=j*dx;
  y2=(j-1)*dx;
  sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
  sigmays=sigmay*(muz/(epsz*eps));
  da1=exp(-sigmays*dt/muz);
  db1=(1-da1)/(sigmays*dx);
  dahzybcb(1:iefbc,j)=da1;
  dbhzybcb(1:iefbc,j)=db1;
  caeybcb(1:ibfbc,j)=capml;
  cbeybcb(1:ibfbc,j)=cbpml;
  cpeybcb(1:ibfbc,j)=cppml;
  dahzxbcb(1:iefbc,j)=da;
  dbhzxbcb(1:iefbc,j)=db;

  dapybcb(1:ibfbc,j)=dap;
  dbpybcb(1:ibfbc,j)=dbp;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     LEFT region
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

caeybcl(1,1:je)=1.0;   %0.0%1.0;
cbeybcl(1,1:je)=0.0;
cpeybcl(1,1:je)=0.0;
   dapybcl(1,1:je)=1.0;
  dbpybcl(1,1:je)=0.0;
bcfactor=eps*sigmam/(dx*(delbc^orderbc)*(orderbc+1));
for i=2:iebc
  x1=(iebc-i+1.5)*dx;
  x2=(iebc-i+0.5)*dx;
  sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
   eaf1=dt*(sigmax+sigdisp)/(epsz*eps);
  ca1=exp(-eaf1);
  cb1=(1.0-ca1)/((sigmax+sigdisp)*dx);
  cp1=(1.0-ca1)/((sigmax+sigdisp)*tau0);
  caeybcl(i,1:je)=ca1;
  cbeybcl(i,1:je)=cb1;
  cpeybcl(i,1:je)=cp1;
  dapybcl(i,1:je)=dap;
  dbpybcl(i,1:je)=dbp;
  caeybcf(i,1:jebc)=ca1;
  cbeybcf(i,1:jebc)=cb1;
  cpeybcf(i,1:jebc)=cp1;
  dapybcf(i,1:jebc)=dap;
  dbpybcf(i,1:jebc)=dbp;
  caeybcb(i,1:jebc)=ca1;
  cbeybcb(i,1:jebc)=cb1;
  cpeybcb(i,1:jebc)=cp1;
   dapybcb(i,1:jebc)=dap;
  dbpybcb(i,1:jebc)=dbp;
end
sigmax=bcfactor*(0.5*dx)^(orderbc+1);
   eaf1=dt*(sigmax+sigdisp)/(epsz*eps);
  ca1=exp(-eaf1);
  cb1=(1.0-ca1)/((sigmax+sigdisp)*dx);
  cp1=(1.0-ca1)/((sigmax+sigdisp)*tau0);
caey(1,1:je)=ca1;
cbey(1,1:je)=cb1;
cpey(1,1:je)=cp1;
dapy(1,1:je)=dap;
dbpy(1,1:je)=dbp;
caeybcf(iebc+1,1:jebc)=ca1;
cbeybcf(iebc+1,1:jebc)=cb1;
cpeybcf(iebc+1,1:jebc)=cp1;
dapybcf(iebc+1,1:jebc)=dap;
  dbpybcf(iebc+1,1:jebc)=dbp;
caeybcb(iebc+1,1:jebc)=ca1;
cbeybcb(iebc+1,1:jebc)=cb1;
cpeybcb(iebc+1,1:jebc)=cp1;
dapybcb(iebc+1,1:jebc)=dap;
  dbpybcb(iebc+1,1:jebc)=dbp;
for i=1:iebc
  x1=(iebc-i+1)*dx;
  x2=(iebc-i)*dx;
  sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
  sigmaxs=sigmax*(muz/(epsz*eps));
  da1=exp(-sigmaxs*dt/muz);
  db1=(1-da1)/(sigmaxs*dx);
  dahzxbcl(i,1:je)=da1;
  dbhzxbcl(i,1:je)=db1;

  dahzxbcf(i,1:jebc)=da1;
  dbhzxbcf(i,1:jebc)=db1;

  dahzxbcb(i,1:jebc)=da1;
  dbhzxbcb(i,1:jebc)=db1;

  caexbcl(i,2:je)=capml;
  cbexbcl(i,2:je)=cbpml;
  cpexbcl(i,2:je)=cppml;
   dapxbcl(i,2:je)=dap;
  dbpxbcl(i,2:je)=dbp;
  dahzybcl(i,1:je)=da;
  dbhzybcl(i,1:je)=db;

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     RIGHT region
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

caeybcr(ibbc,1:je)=1.0;     %0.0%1.0;
cbeybcr(ibbc,1:je)=0.0;
cpeybcr(ibbc,1:je)=0.0;
dapybcr(ibbc,1:je)=1.0;
dbpybcr(ibbc,1:je)=0.0;
bcfactor=eps*sigmam/(dx*(delbc^orderbc)*(orderbc+1));
for i=2:iebc
  x1=(i-0.5)*dx;
  x2=(i-1.5)*dx;
  sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
   eaf1=dt*(sigmax+sigdisp)/(epsz*eps);
  ca1=exp(-eaf1);
  cb1=(1.0-ca1)/((sigmax+sigdisp)*dx);
  cp1=(1.0-ca1)/((sigmax+sigdisp)*tau0);
  caeybcr(i,1:je)=ca1;
  cbeybcr(i,1:je)=cb1;
  cpeybcr(i,1:je)=cp1;

    dapybcr(i,1:je)=dap;
  dbpybcr(i,1:je)=dbp;

  caeybcf(i+iebc+ie,1:jebc)=ca1;
  cbeybcf(i+iebc+ie,1:jebc)=cb1;
  cpeybcf(i+iebc+ie,1:jebc)=cp1;

    dapybcf(i+iebc+ie,1:jebc)=dap;
  dbpybcf(i+iebc+ie,1:jebc)=dbp;

  caeybcb(i+iebc+ie,1:jebc)=ca1;
  cbeybcb(i+iebc+ie,1:jebc)=cb1;
  cpeybcb(i+iebc+ie,1:jebc)=cp1;

  dapybcb(i+iebc+ie,1:jebc)=dap;
  dbpybcb(i+iebc+ie,1:jebc)=dbp;
end
sigmax=bcfactor*(0.5*dx)^(orderbc+1);
eaf1=dt*(sigmay+sigdisp)/(epsz*eps);
  ca1=exp(-eaf1);
  cb1=(1.0-ca1)/((sigmax+sigdisp)*dx);
  cp1=(1.0-ca1)/((sigmax+sigdisp)*tau0);
caey(ib,1:je)=ca1;
cbey(ib,1:je)=cb1;
cpey(ib,1:je)=cp1;
dapy(ib,1:je)=dap;
dbpy(ib,1:je)=dbp;

caeybcf(iebc+ib,1:jebc)=ca1;
cbeybcf(iebc+ib,1:jebc)=cb1;
cpeybcf(iebc+ib,1:jebc)=cp1;
dapybcf(iebc+ib,1:jebc)=dap;
dbpybcf(iebc+ib,1:jebc)=dbp;

caeybcb(iebc+ib,1:jebc)=ca1;
cbeybcb(iebc+ib,1:jebc)=cb1;
cpeybcb(iebc+ib,1:jebc)=cp1;
dapybcb(iebc+ib,1:jebc)=dap;
dbpybcb(iebc+ib,1:jebc)=dbp;
for i=1:iebc
  x1=i*dx;
  x2=(i-1)*dx;
  sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
  sigmaxs=sigmax*(muz/(epsz*eps));
  da1=exp(-sigmaxs*dt/muz);
  db1=(1-da1)/(sigmaxs*dx);
  dahzxbcr(i,1:je) = da1;
  dbhzxbcr(i,1:je) = db1;

  dahzxbcf(i+ie+iebc,1:jebc)=da1;
  dbhzxbcf(i+ie+iebc,1:jebc)=db1;

  dahzxbcb(i+ie+iebc,1:jebc)=da1;
  dbhzxbcb(i+ie+iebc,1:jebc)=db1;

   dapxbcb(i,2:je)=dap;
  dbpxbcb(i,2:je)=dbp;
  caexbcr(i,2:je)=capml;
  cbexbcr(i,2:je)=cbpml;
  cpexbcr(i,2:je)=cppml;

  dahzybcr(i,1:je)=da;
  dbhzybcr(i,1:je)=db;
   end

%***********************************************************************
%     BEGIN TIME-STEPPING LOOP
%***********************************************************************
for n=1:nmax

%***********************************************************************
%     Update electric fields (EX and EY) in main grid
%***********************************************************************
ex(:,2:je)=caex(:,2:je).*ex(:,2:je)+...
           cbex(:,2:je).*(hz(:,2:je)-hz(:,1:je-1))+cpex(:,2:je).*px(:,2:je);

ey(2:ie,:)=caey(2:ie,:).*ey(2:ie,:)+...
           cbey(2:ie,:).*(hz(1:ie-1,:)-hz(2:ie,:))+cpey(2:ie,:).*py(2:ie,:);

%***********************************************************************
%     Update EX in PML regions
%***********************************************************************

%     FRONT

exbcf(:,2:jebc)=caexbcf(:,2:jebc).*exbcf(:,2:jebc)-...
  cbexbcf(:,2:jebc).*(hzxbcf(:,1:jebc-1)+hzybcf(:,1:jebc-1)-...
                      hzxbcf(:,2:jebc)-hzybcf(:,2:jebc))+...
    cpexbcf(:,2:jebc).*pxbcf(:,2:jebc);
ex(1:ie,1)=caex(1:ie,1).*ex(1:ie,1)-...
  cbex(1:ie,1).*(hzxbcf(ibbc:iebc+ie,jebc)+...
                hzybcf(ibbc:iebc+ie,jebc)-hz(1:ie,1))+...
      cpex(1:ie,1).*px(1:ie,1);

%     BACK

exbcb(:,2:jebc-1)=caexbcb(:,2:jebc-1).*exbcb(:,2:jebc-1)-...
  cbexbcb(:,2:jebc-1).*(hzxbcb(:,1:jebc-2)+hzybcb(:,1:jebc-2)-...
                        hzxbcb(:,2:jebc-1)-hzybcb(:,2:jebc-1))+...
   cpexbcb(:,2:jebc-1).*pxbcb(:,2:jebc-1);
ex(1:ie,jb)=caex(1:ie,jb).*ex(1:ie,jb)-...
  cbex(1:ie,jb).*(hz(1:ie,jb-1)-hzxbcb(ibbc:iebc+ie,1)-...
                 hzybcb(ibbc:iebc+ie,1))+...
                 cpex(1:ie,jb).*px(1:ie,jb);

%     LEFT

exbcl(:,2:je)=caexbcl(:,2:je).*exbcl(:,2:je)-...
  cbexbcl(:,2:je).*(hzxbcl(:,1:je-1)+hzybcl(:,1:je-1)-...
                    hzxbcl(:,2:je)-hzybcl(:,2:je))+...
                    cpexbcl(:,2:je).*pxbcl(:,2:je);
exbcl(:,1)=caexbcl(:,1).*exbcl(:,1)-...
  cbexbcl(:,1).*(hzxbcf(1:iebc,jebc)+hzybcf(1:iebc,jebc)-...
                 hzxbcl(:,1)-hzybcl(:,1))+cpexbcl(:,1).*pxbcl(:,1);
exbcl(:,jb)=caexbcl(:,jb).*exbcl(:,jb)-...
  cbexbcl(:,jb).*(hzxbcl(:,je)+hzybcl(:,je)-...
                  hzxbcb(1:iebc,1)-hzybcb(1:iebc,1))+...
          cpexbcl(:,jb).*pxbcl(:,je);

%     RIGHT

exbcr(:,2:je)=caexbcr(:,2:je).*exbcr(:,2:je)-...
  cbexbcr(:,2:je).*(hzxbcr(:,1:je-1)+hzybcr(:,1:je-1)-...
                    hzxbcr(:,2:je)-hzybcr(:,2:je))+...
       cpexbcr(:,2:je).*pxbcr(:,2:je);
exbcr(:,1)=caexbcr(:,1).*exbcr(:,1)-...
  cbexbcr(:,1).*(hzxbcf(1+iebc+ie:iefbc,jebc)+...
                 hzybcf(1+iebc+ie:iefbc,jebc)-...
                 hzxbcr(:,1)-hzybcr(:,1))+cpexbcr(:,1).*pxbcr(:,1);
exbcr(:,jb)=caexbcr(:,jb).*exbcr(:,jb)-...
  cbexbcr(:,jb).*(hzxbcr(:,je)+hzybcr(:,je)-...
                  hzxbcb(1+iebc+ie:iefbc,1)-...
                  hzybcb(1+iebc+ie:iefbc,1))+cpexbcr(:,jb).*pxbcr(:,jb);

%***********************************************************************
%     Update EY in PML regions
%***********************************************************************

%     FRONT

eybcf(2:iefbc,:)=caeybcf(2:iefbc,:).*eybcf(2:iefbc,:)-...
  cbeybcf(2:iefbc,:).*(hzxbcf(2:iefbc,:)+hzybcf(2:iefbc,:)-...
                       hzxbcf(1:iefbc-1,:)-hzybcf(1:iefbc-1,:))+...
       cpeybcf(2:iefbc,:).*pybcf(2:iefbc,:);

%     BACK

eybcb(2:iefbc,:)=caeybcb(2:iefbc,:).*eybcb(2:iefbc,:)-...
  cbeybcb(2:iefbc,:).*(hzxbcb(2:iefbc,:)+hzybcb(2:iefbc,:)-...
                       hzxbcb(1:iefbc-1,:)-hzybcb(1:iefbc-1,:))+...
       cpeybcb(2:iefbc,:).*pybcb(2:iefbc,:);

%     LEFT

eybcl(2:iebc,:)=caeybcl(2:iebc,:).*eybcl(2:iebc,:)-...
  cbeybcl(2:iebc,:).*(hzxbcl(2:iebc,:)+hzybcl(2:iebc,:)-...
                      hzxbcl(1:iebc-1,:)-hzybcl(1:iebc-1,:))+...
     cpeybcl(2:iebc,:).*pybcl(2:iebc,:);
ey(1,:)=caey(1,:).*ey(1,:)-...
  cbey(1,:).*(hz(1,:)-hzxbcl(iebc,:)-hzybcl(iebc,:))+cpey(1,:).*py(1,:);

%     RIGHT

eybcr(2:iebc,:)=caeybcr(2:iebc,:).*eybcr(2:iebc,:)-...
  cbeybcr(2:iebc,:).*(hzxbcr(2:iebc,:)+hzybcr(2:iebc,:)-...
                      hzxbcr(1:iebc-1,:)-hzybcr(1:iebc-1,:))+...
    cpeybcr(2:iebc,:).*pybcr(2:iebc,:);
ey(ib,:)=caey(ib,:).*ey(ib,:)-...
  cbey(ib,:).*(hzxbcr(1,:)+hzybcr(1,:)- hz(ie,:))+...
  cpey(ib,:).*py(ib,:);

%***********************************************************************
%     Update magnetic fields (HZ) in main grid
%***********************************************************************

% for dispersion consideration P equation (32) 1995 IEEE tran. antennas and
% Propagation Vol.43 No. 4, pp422 (Jeffrey L. Young)
px(:,2:je)=dapx(:,2:je).*px(:,2:je)+...
          dbpx(:,2:je).*ex(:,2:je);

py(2:ie,:)=dapy(2:ie,:).*py(2:ie,:)+...
        dbpy(2:ie,:).*ey(2:ie,:);

hz(1:ie,1:je)=dahz(1:ie,1:je).*hz(1:ie,1:je)+...
              dbhz(1:ie,1:je).*(ex(1:ie,2:jb)-ex(1:ie,1:je)+...
                                ey(1:ie,1:je)-ey(2:ib,1:je));

hz(imp_x,imp_y)=source(n);    %(imp_x ,imp_y) are the location of the point source

%hz(imp_x,imp_y+25)=source(n);

%***********************************************************************
%     Update HZX ,px in PML regions
%***********************************************************************

%     FRONT

pxbcf(:,2:jebc)=dapxbcf(:,2:jebc).*pxbcf(:,2:jebc)+...
    dbpxbcf(:,2:jebc).*exbcf(:,2:jebc);
px(1:ie,1)=dapx(1:ie,1).*px(1:ie,1)+...
     dbpx(1:ie,1).*ex(1:ie,1);

hzxbcf(1:iefbc,:)=dahzxbcf(1:iefbc,:).*hzxbcf(1:iefbc,:)-...
  dbhzxbcf(1:iefbc,:).*(eybcf(2:ibfbc,:)-eybcf(1:iefbc,:));

%     BACK
pxbcb(:,2:jebc-1)=dapxbcb(:,2:jebc-1).*pxbcb(:,2:jebc-1)+...
   dbpxbcb(:,2:jebc-1).*exbcb(:,2:jebc-1);
px(1:ie,jb)=dapx(1:ie,jb).*px(1:ie,jb)+...
                 dbpx(1:ie,jb).*ex(1:ie,jb);

hzxbcb(1:iefbc,:)=dahzxbcb(1:iefbc,:).*hzxbcb(1:iefbc,:)-...
  dbhzxbcb(1:iefbc,:).*(eybcb(2:ibfbc,:)-eybcb(1:iefbc,:));

%     LEFT

pxbcl(:,2:je)=dapxbcl(:,2:je).*pxbcl(:,2:je)+...
                    dbpxbcl(:,2:je).*exbcl(:,2:je);
pxbcl(:,1)=dapxbcl(:,1).*pxbcl(:,1)+dbpxbcl(:,1).*exbcl(:,1);
pxbcl(:,jb)=dapxbcl(:,jb).*pxbcl(:,jb)+dbpxbcl(:,jb).*exbcl(:,je);

hzxbcl(1:iebc-1,:)=dahzxbcl(1:iebc-1,:).*hzxbcl(1:iebc-1,:)-...
  dbhzxbcl(1:iebc-1,:).*(eybcl(2:iebc,:)-eybcl(1:iebc-1,:));
hzxbcl(iebc,:)=dahzxbcl(iebc,:).*hzxbcl(iebc,:)-...
  dbhzxbcl(iebc,:).*(ey(1,:)-eybcl(iebc,:));

%     RIGHT
pxbcr(:,2:je)=dapxbcr(:,2:je).*pxbcr(:,2:je)+dbpxbcr(:,2:je).*exbcr(:,2:je);
pxbcr(:,1)=dapxbcr(:,1).*pxbcr(:,1)+dbpxbcr(:,1).*exbcr(:,1);
pxbcr(:,jb)=dapxbcr(:,jb).*pxbcr(:,jb)+dbpxbcr(:,jb).*exbcr(:,jb);

hzxbcr(2:iebc,:)=dahzxbcr(2:iebc,:).*hzxbcr(2:iebc,:)-...
  dbhzxbcr(2:iebc,:).*(eybcr(3:ibbc,:)-eybcr(2:iebc,:));
hzxbcr(1,:)=dahzxbcr(1,:).*hzxbcr(1,:)-...
  dbhzxbcr(1,:).*(eybcr(2,:)-ey(ib,:));

%***********************************************************************
%     Update HZY ,py in PML regions
%***************************************************** ..

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

图片:未命名.JPG

附件:

FDTD.rar (5 K) 下载次数:13

离线yuyeliuxu

UID ：21110

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

1楼发表于: 2009-07-16 10:20:08

程序是按照论文编写的，现在才发现跟taflove的书有些不一样，论文理解也有些问题

附件:

propagation in linear dispersive media.pdf (509 K) 下载次数:13

发帖回复