1411阅读
0回复

[求助]求牛人帮我改一下2D-TM的MATLAB程序

UID ：55468

注册：2010-03-23
登录：2010-10-18
发帖：9
等级：旁观者

0楼发表于: 2010-04-20 10:11:17

此悬赏帖已过期
最佳答案：5 rf币，热心助人剩余点数: 1 rf币。

想问一下哪些参数要改，我的递推有没有问题，PML的设置有没有问题，还有画的图总是不对，如果输出的场为数值的话不是0就是NAN，谁能帮我改一下呢，谢啦O(∩_∩)O~

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%      fundamental constants (in the dimesion of m)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
lambda=1.55e-6;             %center wavelength of source excitation
freq=cc/lambda;                %center frequency of source excitation
omega=2.0*pi*freq;
%***********************************************************************
%     Grid parameters (in the dimesion of m)
%***********************************************************************
ie=320;           %number of grid cells in x-direction
je=44;            %number of grid cells in y-direction
ib=ie+1;
jb=je+1;
is=5;            %location of x-directed hard source
js=je/2;          %location of x-directed hard source
dx=5.0e-9;        %space increment of square lattice
dt=lambda/(pi*cc);   %time step
nmax=3300;         %total number of time steps
iebc=8;           %thickness of left and right PML region
jebc=8;           %thickness of front and back PML region
ibbc=iebc+1;
jbbc=jebc+1;
iefbc=ie+2*iebc;
jefbc=je+2*jebc;
ibfbc=iefbc+1;
jbfbc=jefbc+1;
%***********************************************************************
%     Material parameters
%***********************************************************************
media=2;
omegap=1.34639e16;
epsm1=1.999+omegap^2/((9.61712e-13i)*omega-omega^2);
epsm=real(epsm1);
eps=[1.0 epsm ];
sig=[0.0 8.85e-18];
mur=[1.0 4.0*pi*1.0e-7];
sim=[0.0 0.0];
%***********************************************************************
%     Wave excitation
%***********************************************************************
rtau=10*lambda/pi/cc;
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
subplot(2,1,1)
plot(source)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                fields
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ez=zeros(ib,jb);
hx=zeros(ib,je);            %fields in main grid
hy=zeros(ie,jb);
ezbcf=zeros(ibfbc,jebc);
ezxbcf=zeros(ibfbc,jebc);
ezybcf=zeros(ibfbc,jebc);
hxbcf=zeros(ibfbc,jebc);       %fields in front region
hybcf=zeros(iefbc,jebc);
ezbcb=zeros(ibfbc,jebc);
ezxbcb=zeros(ibfbc,jebc);
ezybcb=zeros(ibfbc,jebc);
hxbcb=zeros(ibfbc,jebc);       %fields in back region
hybcb=zeros(iefbc,jebc);
ezbcl=zeros(iebc,jb);
ezxbcl=zeros(iebc,jb);
ezybcl=zeros(iebc,jb);
hxbcl=zeros(iebc,je);        %fields in left region
hybcl=zeros(iebc,jb);
ezbcr=zeros(iebc,jb);
ezxbcr=zeros(iebc,jb);
ezybcr=zeros(iebc,jb);
hxbcr=zeros(iebc,je);        %fields in right region
hybcr=zeros(iebc,jb);
%************************************************************************
%ezxbcfn(1:ibfbc,1:jebc)=ezxbcf(1:ibfbc,1:jebc);
%ezybcfn(1:ibfbc,1:jebc)=ezybcf(1:ibfbc,1:jebc);
%hxbcfn(1:ibfbc,1:jebc)=hxbcf(1:ibfbc,1:jebc);       %fields in front region******************************************************************************************
%hybcfn(1:iefbc,1:jebc)=hybcf(1:iefbc,1:jebc);

%ezxbcbn(1:ibfbc,1:jebc)=ezxbcb(1:ibfbc,1:jebc);
%ezybcbn(1:ibfbc,1:jebc)=ezybcb(1:ibfbc,1:jebc);
%hxbcbn(1:ibfbc,1:jebc)=hxbcb(1:ibfbc,1:jebc);       %fields in back region********************************************************************************************
%hybcbn(1:iefbc,1:jebc)=hybcb(1:iefbc,1:jebc);

%ezxbcln(1:iebc,1:jb)=ezxbcl(1:iebc,1:jb);
%ezybcln(1:iebc,1:jb)=ezybcl(1:iebc,1:jb);
%hxbcln(1:iebc,1:je)=hxbcl(1:iebc,1:je);        %fields in left region
%hybcln(1:iebc,1:jb)=hybcl(1:iebc,1:jb);

%ezxbcrn(1:iebc,1:jb)=ezxbcr(1:iebc,1:jb);
%ezybcrn(1:iebc,1:jb)=ezybcr(1:iebc,1:jb);
%hxbcrn(1:iebc,1:je)=hxbcr(1:iebc,1:je);        %fields in right region
%hybcrn(1:iebc,1:jb)=hybcr(1:iebc,1:jb);
%***********************************************************************
%     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));
cp(i)=(1.0-haf)/(1.0+haf);
cq(i)=dt/muz/mur(i)./dx./(1.0+haf);
end
%for i=2
%eaf =dt*sig(i)./(2.0*epsz*eps(i));
%ca2=(1.0-eaf)/(1.0+eaf);
%cb2=dt/epsz/eps(i)./dx/(1.0+eaf);
%haf =dt*sim(i)./(2.0*muz*mur(i));%
% cp2=(1.0-haf)/(1.0+haf);
%cq2=dt/muz/mur(i)./dx./(1.0+haf);
%end

%***********************************************************************
%     Geometry specification (main grid)
%***********************************************************************
%     Initialize entire main grid to free space
caez(1:ib,1:jb)=ca1;
cbez(1:ib,1:jb)=cb1;
cphx(1:ib,1:je)=cp1;
cqhx(1:ib,1:je)=cq1;
cphy(1:ie,1:jb)=cp1;
cqhy(1:ie,1:jb)=cq1;
% boundary of main zone
sigmamax=1e-17;               %????????????????????????????????
%for i=140:180
    for j=3
eaf =dt*(sig(2)+sigmamax/64)/(2.0*epsz*eps(2));
caez(140:180,j)=(1.0-eaf)/(1.0+eaf);
cbez(140:180,j)=dt/epsz/eps(2)/dx/(1.0+eaf);
    end
%end
%for i=140:180
    for j=42
eaf =dt*(sig(2)+sigmamax/64)/(2.0*epsz*eps(2));
caez(140:180,j)=(1.0-eaf)/(1.0+eaf);
cbez(140:180,j)=dt/epsz/eps(2)/dx/(1.0+eaf);
    end
%end
for i=140
    %for j=3:42
eaf =dt*(sig(2)+sigmamax/64)/(2.0*epsz*eps(2));
caez(i,3:42)=(1.0-eaf)/(1.0+eaf);
cbez(i,3:42)=dt/epsz/eps(2)/dx/(1.0+eaf);
    %end
end
for i=180
   % for j=3:42
eaf =dt*(sig(2)+sigmamax/64)/(2.0*epsz*eps(2));
caez(i,3:42)=(1.0-eaf)/(1.0+eaf);
cbez(i,3:42)=dt/epsz/eps(2)/dx/(1.0+eaf);
    %end
end
%     add SMG
for i=141:179
    for j=4:41
        caez(i,j)=ca2;
        cbez(i,j)=cb2;
        cphx(i,j)=cp2;
        cqhx(i,j)=cq2;
        cphy(i,j)=cp2;
        cqhy(i,j)=cq2;
    end
end

%     Fill the PML regions
caezxbcf=zeros(ibfbc,jebc);
cbezxbcf=zeros(ibfbc,jebc);
caezybcf=zeros(ibfbc,jebc);
cbezybcf=zeros(ibfbc,jebc);
cphxbcf=zeros(ibfbc,jebc);       %fields in front region
cqhxbcf=zeros(ibfbc,jebc);
cphybcf=zeros(iefbc,jebc);
cqhybcf=zeros(iefbc,jebc);
caezxbcb=zeros(ibfbc,jebc);
cbezxbcb=zeros(ibfbc,jebc);
caezybcb=zeros(ibfbc,jebc);
cbezybcb=zeros(ibfbc,jebc);
cphxbcb=zeros(ibfbc,jebc);       %fields in back region
cqhxbcb=zeros(ibfbc,jebc);
cphybcb=zeros(iefbc,jebc);
cqhybcb=zeros(iefbc,jebc);
caezxbcl=zeros(iebc,jb);
cbezxbcl=zeros(iebc,jb);
caezybcl=zeros(iebc,jb);
cbezybcl=zeros(iebc,jb);
cphxbcl=zeros(iebc,je);        %fields in left region
cqhxbcl=zeros(iebc,je);
cphybcl=zeros(iebc,jb);
cqhybcl=zeros(iebc,jb);

caezxbcr=zeros(iebc,jb);
cbezxbcr=zeros(iebc,jb);
caezybcr=zeros(iebc,jb);
cbezybcr=zeros(iebc,jb);
cphxbcr=zeros(iebc,je);        %fields in right region
cqhxbcr=zeros(iebc,je);
cphybcr=zeros(iebc,jb);
cqhybcr=zeros(iebc,jb);
%a=ones(iebc,je);
%***********************************************************************
%     Fill the PML regions
%***********************************************************************
rmax=0.00001;
orderbc=2;
delbc=iebc*dx;
sigmam=-log(rmax/100.0)*epsz*cc*(orderbc+1)/(2*delbc);
bcfactor=eps(1)*sigmam/(dx*(delbc^orderbc)*(orderbc+1));
rmax=0.00001;
orderbc=2;
%     FRONT region
caexbcf(1:iefbc,1)=1.0;
cbexbcf(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));
ca1=exp(-sigmay*dt/(epsz*eps(1)));
cb1=(1.0-ca1)/(sigmay*dx);
caexbcf(1:iefbc,j)=ca1;
cbexbcf(1:iefbc,j)=cb1;
end
sigmay = bcfactor*(0.5*dx)^(orderbc+1);
ca1=exp(-sigmay*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmay*dx);
caex(1:ie,1)=ca1;
cbex(1:ie,1)=cb1;
caexbcl(1:iebc,1)=ca1;
cbexbcl(1:iebc,1)=cb1;
caexbcr(1:iebc,1)=ca1;
cbexbcr(1:iebc,1)=cb1;
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(1)));
cp1=exp(-sigmays*dt/muz);
cq1=(1-cp1)/(sigmays*dx);
cphzybcf(1:iefbc,j)=cp1;
cqhzybcf(1:iefbc,j)=cq1;
caeybcf(1:ibfbc,j)=ca(1);
cbeybcf(1:ibfbc,j)=cb(1);
cphzxbcf(1:iefbc,j)=cp(1);
cqhzxbcf(1:iefbc,j)=cq(1);
end
%     BACK region
caexbcb(1:iefbc,jbbc)=1.0;
cbexbcb(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));
ca1=exp(-sigmay*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmay*dx);
caexbcb(1:iefbc,j)=ca1;
cbexbcb(1:iefbc,j)=cb1;
end
sigmay = bcfactor*(0.5*dx)^(orderbc+1);
ca1=exp(-sigmay*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmay*dx);
caex(1:ie,jb)=ca1;
cbex(1:ie,jb)=cb1;
caexbcl(1:iebc,jb)=ca1;
cbexbcl(1:iebc,jb)=cb1;
caexbcr(1:iebc,jb)=ca1;
cbexbcr(1:iebc,jb)=cb1;
for j=1:jebc
y1=j*dx;
y2=(j-1)*dx;
sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
sigmays=sigmay*(muz/(epsz*eps(1)));
cp1=exp(-sigmays*dt/muz);
ccq1=(1-cp1)/(sigmays*dx);
cphzybcb(1:iefbc,j)=cp1;
cqhzybcb(1:iefbc,j)=cq1;
caeybcb(1:ibfbc,j)=ca(1);
cbeybcb(1:ibfbc,j)=cb(1);
cphzxbcb(1:iefbc,j)=cp(1);
cqhzxbcb(1:iefbc,j)=cq(1);
end
%     LEFT region
caeybcl(1,1:je)=1.0;
cbeybcl(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));
ca1=exp(-sigmax*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmax*dx);
caeybcl(i,1:je)=ca1;
cbeybcl(i,1:je)=cb1;
caeybcf(i,1:jebc)=ca1;
cbeybcf(i,1:jebc)=cb1;
caeybcb(i,1:jebc)=ca1;
cbeybcb(i,1:jebc)=cb1;
end
sigmax=bcfactor*(0.5*dx)^(orderbc+1);
ca1=exp(-sigmax*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmax*dx);
caey(1,1:je)=ca1;
cbey(1,1:je)=cb1;
caeybcf(iebc+1,1:jebc)=ca1;
cbeybcf(iebc+1,1:jebc)=cb1;
caeybcb(iebc+1,1:jebc)=ca1;
cbeybcb(iebc+1,1:jebc)=cb1;
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(1)));
cp1=exp(-sigmaxs*dt/muz);
cq1=(1-da1)/(sigmaxs*dx);
cphzxbcl(i,1:je)=cp1;
cqhzxbcl(i,1:je)=cq1;
cphzxbcf(i,1:jebc)=cp1;
cqhzxbcf(i,1:jebc)=cq1;
cphzxbcb(i,1:jebc)=cp1;
cqhzxbcb(i,1:jebc)=cq1;
caexbcl(i,2:je)=ca(1);
cbexbcl(i,2:je)=cb(1);
cphzybcl(i,1:je)=cp(1);
cqhzybcl(i,1:je)=cq(1);
end


%***********************************************************************
%     Movie initialization
%***********************************************************************
subplot(2,1,2),pcolor(ez);
shading flat;
caxis([-1000000000 1000000000]);
axis([1 ib 1 jb]);
%axis([1 jb 1 ib]);
colorbar;
axis image;
axis off;
title(('ez at time step = 0'));

%***********************************************************************
%     BEGIN TIME-STEPPING LOOP
%***********************************************************************
ezbcl(1:iebc,1:jb)=ezxbcl(1:iebc,1:jb)+ezybcl(1:iebc,1:jb);
    ezbcr(1:iebc,1:jb)=ezxbcr(1:iebc,1:jb)+ezybcr(1:iebc,1:jb);
    ezbcf(1:ibfbc,1:jebc)=ezxbcf(1:ibfbc,1:jebc)+ezybcf(1:ibfbc,1:jebc);
    ezbcb(1:ibfbc,1:jebc)=ezxbcb(1:ibfbc,1:jebc)+ezybcb(1:ibfbc,1:jebc);

for n=1:nmax
%***********************************************************************
%     Update electric fields (ez,hx,hy) in main grid
%***********************************************************************

    x=source(n).*10e20;

     ez(is,js)=ez(is,js)+x;%
    ez(2:ie,2:je)=caez(2:ie,2:je).*ez(2:ie,2:je)+...
        cbez(2:ie,2:je).*((hy(2:ie,2:je)-hy(1:ie-1,2:je))-...
        (hx(2:ie,2:je)-hx(2:ie,1:je-1)))

    hx(1:ie,1:je)=cphx(1:ie,1:je).*hx(1:ie,1:je)+...
        cqhx(1:ie,1:je).*(ez(1:ie,2:jb)-ez(1:ie,1:je));
    hy(1:ie,1:jb)=cphy(1:ie,1:jb).*hy(1:ie,1:jb)+...
        cqhy(1:ie,1:jb).*(ez(2:ib,1:jb)-ez(1:ie,1:jb));

ezxbcfn(1:ibfbc,1:jebc)=ezxbcf(1:ibfbc,1:jebc);
ezybcfn(1:ibfbc,1:jebc)=ezybcf(1:ibfbc,1:jebc);
%ezybcf(2:iefbc,2)=ezybcf(2:iefbc,2);       %fields in front region
ezbcfn(1:ibfbc,1:jebc)=ezbcf(1:ibfbc,1:jebc);

ezxbcbn(1:ibfbc,1:jebc)=ezxbcb(1:ibfbc,1:jebc);
ezybcbn(1:ibfbc,1:jebc)=ezybcb(1:ibfbc,1:jebc);
%ezybcbn(2:iefbc,jebc-1)=ezybcb(2:iefbc,jebc-1);       %fields in back region
ezbcbn(1:ibfbc,1:jebc)=ezbcb(1:ibfbc,1:jebc);
%ezbcln(iebc-1,1)=ezbcl(iebc-1,1);
%ezbcln(2,jb)=ezbcln(2,jb);
%ezbcf(ibbc,1)=ezbcbn(ibbc-1,2)
ezxbcln(1:iebc,1:jb)=ezxbcl(1:iebc,1:jb);
ezybcln(1:iebc,1:jb)=ezybcl(1:iebc,1:jb);        %fields in left region
ezbcln(1:iebc,1:jb)=ezbcl(1:iebc,1:jb);
%ezbcf(1,1)=ezbcfn(2,2)
ezxbcrn(1:iebc,1:jb)=ezxbcr(1:iebc,1:jb);
ezybcrn(1:iebc,1:jb)=ezybcr(1:iebc,1:jb);       %fields in right region
ezbcrn(1:iebc,1:jb)=ezbcr(1:iebc,1:jb);

    % left
    hxbcl(1:iebc,1:je)=cphxbcl(1:iebc,1:je).*hxbcl(1:iebc,1:je)+...
        cqhxbcl(1:iebc,1:je).*(ezbcl(1:iebc,1:je)-ezbcl(1:iebc,2:jb));
    hybcl(1:iebc-1,1:jb)=cphybcl(1:iebc-1,1:jb).*hybcl(1:iebc-1,1:jb)+...
        cqhybcl(1:iebc-1,1:jb).*(ezbcl(2:iebc,1:jb)-ezbcl(1:iebc-1,1:jb));
    ezxbcl(2:iebc,1:jb)=caezxbcl(2:iebc,1:jb).*ezxbcl(2:iebc,1:jb)+...
        cbezxbcl(2:iebc,1:jb).*(hybcl(2:iebc,1:jb)-hybcl(1:iebc-1,1:jb));
    ezybcl(1:iebc,2:je)=caezybcl(1:iebc,2:je).*ezybcl(1:iebc,2:je)+...
        cbezybcl(1:iebc,2:je).*(hxbcl(1:iebc,2:je)-hxbcl(1:iebc,1:je-1));
    % right
    hxbcr(1:iebc,1:je)=cphxbcr(1:iebc,1:je).*hxbcr(1:iebc,1:je)+...
        cqhxbcr(1:iebc,1:je).*(ezbcr(1:iebc,1:je)-ezbcr(1:iebc,2:jb));
    hybcr(2:iebc,1:jb)=cphybcr(2:iebc,1:jb).*hybcr(2:iebc,1:jb)+...
        cqhybcr(2:iebc,1:jb).*(ezbcr(2:iebc,1:jb)-ezbcr(1:iebc-1,1:jb));
    ezxbcr(1:iebc-1,1:jb)=caezxbcr(1:iebc-1,1:jb).*ezxbcr(1:iebc-1,1:jb)+...
        cbezxbcl(1:iebc-1,1:jb).*(hybcr(2:iebc,1:jb)-hybcr(1:iebc-1,1:jb));
    ezybcr(1:iebc,2:je)=caezybcr(1:iebc,2:je).*ezybcr(1:iebc,2:je)+...
        cbezybcl(1:iebc,2:je).*(hxbcr(1:iebc,2:je)-hxbcr(1:iebc,1:je-1));
    % down
    hxbcf(1:ibfbc,1:jebc-1)=cphxbcf(1:ibfbc,1:jebc-1).*hxbcf(1:ibfbc,1:jebc-1)+...
        cqhxbcf(1:ibfbc,1:jebc-1).*(ezbcf(1:ibfbc,2:jebc)-ezbcf(1:ibfbc,1:jebc-1));
    hybcf(1:iefbc,1:jebc)=cphybcf(1:iefbc,1:jebc).*hybcf(1:iefbc,1:jebc)+...
        cqhybcf(1:iefbc,1:jebc).*(ezbcf(2:ibfbc,1:jebc)-ezbcf(1:iefbc,1:jebc));
    ezxbcf(2:iefbc,1:jebc)=caezxbcf(2:iefbc,1:jebc).*ezxbcf(2:iefbc,1:jebc)+...
        cbezxbcf(2:iefbc,1:jebc).*(hybcf(2:iefbc,1:jebc)-hybcf(1:iefbc-1,1:jebc));
    ezybcf(1:ibbc,2:jebc)=caezybcf(1:ibbc,2:jebc).*ezybcf(1:ibbc,2:jebc)+...
        cbezybcf(1:ibbc,2:jebc).*(hxbcf(1:ibbc,2:jebc)-hxbcf(1:ibbc,1:jebc-1));
    % up
    hxbcb(1:ibfbc,2:jebc)=cphxbcb(1:ibfbc,2:jebc).*hxbcb(1:ibfbc,2:jebc)+...
        cqhxbcb(1:ibfbc,2:jebc).*(ezbcb(1:ibfbc,2:jebc)-ezbcb(1:ibfbc,1:jebc-1));
    hybcb(1:iefbc,1:jebc)=cphxbcb(1:iefbc,1:jebc).*hxbcb(1:iefbc,1:jebc)+...
        cqhybcb(1:iefbc,1:jebc).*(ezbcb(2:ibfbc,1:jebc)-ezbcb(1:iefbc,1:jebc));
    ezxbcb(2:iefbc,1:jebc)=caezxbcb(2:iefbc,1:jebc).*ezxbcb(2:iefbc,1:jebc)+...
        cbezxbcb(2:iefbc,1:jebc).*(hybcb(2:iefbc,1:jebc)-hybcb(1:iefbc-1,1:jebc));
    ezybcb(1:ibbc,1:jebc-1)=caezybcb(1:ibbc,1:jebc-1).*ezybcb(1:ibbc,1:jebc-1)+...
        cbezybcb(1:ibbc,1:jebc-1).*(hxbcb(1:ibbc,2:jebc)-hxbcb(1:ibbc,1:jebc-1));


% ezx
% left line
ezbcl(1,2:je)=ezbcln(2,2:je)+(cc*dt-dx)/(cc*dt+dx)*(ezbcl(2,2:je)-...
    ezbcl(1,2:je))+cc^2*muz*dt/(2*(cc*dt+dx))*(hxbcl(1,2:je)-...
    hxbcl(1,1:je-1)+hxbcl(2,2:je)-hxbcl(2,1:je-1));
ezbcf(1,2:jebc)=ezbcfn(2,2:jebc)+(cc*dt-dx)/(cc*dt+dx)*(ezbcf(2,2:jebc)-...
    ezbcf(1,2:jebc))+cc^2*muz*dt/(2*(cc*dt+dx))*(hxbcf(1,2:jebc)-...
    hxbcf(1,1:jebc-1)+hxbcf(2,2:jebc)-hxbcf(2,1:jebc-1));
ezbcb(1,1:jebc-1)=ezbcbn(2,1:jebc-1)+(cc*dt-dx)/(cc*dt+dx)*(ezbcb(2,1:jebc-1)-...
    ezbcb(1,1:jebc-1))+cc^2*muz*dt/(2*(cc*dt+dx))*(hxbcb(1,2:jebc)-...
    hxbcb(1,1:jebc-1)+hxbcb(2,2:jebc)-hxbcb(2,1:jebc-1));
ezbcl(1,1)=ezbcln(2,1)+(cc*dt-dx)/(cc*dt+dx)*(ezbcl(2,1)-...
    ezbcl(1,1))+cc^2*muz*dt/(2*(cc*dt+dx))*(hxbcl(1,1)-...
    hxbcf(1,jebc)+hxbcl(2,1)-hxbcf(2,jebc));
ezbcl(1,jb)=ezbcln(2,jb)+(cc*dt-dx)/(cc*dt+dx)*(ezbcl(2,jb)-...
    ezbcl(1,jb))+cc^2*muz*dt/(2*(cc*dt+dx))*(hxbcb(1,1)-...
    hxbcl(1,je)+hxbcb(2,1)-hxbcl(2,je));
% right line
ezbcr(iebc,2:je)=ezbcrn(iebc-1,2:je)+(cc*dt-dx)/(cc*dt+dx)*(ezbcr(iebc-1,2:je)-...
    ezbcr(iebc,2:je))+cc^2*muz*dt/(2*(cc*dt+dx))*(hxbcr(iebc,2:je)-...
    hxbcr(iebc,1:je-1)+hxbcr(iebc-1,2:je)-hxbcr(iebc-1,1:je-1));
ezbcf(ibfbc,2:jebc)=ezbcfn(ibfbc-1,2:jebc)+...
    (cc*dt-dx)/(cc*dt+dx)*(ezbcf(ibfbc-1,2:jebc)-...
    ezbcf(ibfbc,2:jebc))+cc^2*muz*dt/(2*(cc*dt+dx))*(hxbcf(ibfbc,2:jebc)-...
    hxbcf(ibfbc,1:jebc-1)+hxbcf(ibfbc-1,2:jebc)-hxbcf(ibfbc-1,1:jebc-1));
ezbcb(ibfbc,1:jebc-1)=ezbcbn(ibfbc-1,1:jebc-1)+...
    (cc*dt-dx)/(cc*dt+dx)*(ezbcb(ibfbc-1,1:jebc-1)-...
    ezbcb(ibfbc,1:jebc-1))+cc^2*muz*dt/(2*(cc*dt+dx))*(hxbcb(ibfbc,2:jebc)-...
    hxbcb(ibfbc,1:jebc-1)+hxbcb(ibfbc-1,2:jebc)-hxbcb(ibfbc-1,1:jebc-1));
ezbcr(iebc,1)=ezbcln(iebc-1,1)+(cc*dt-dx)/(cc*dt+dx)*(ezbcl(iebc-1,1)-...
    ezbcl(iebc,1))+cc^2*muz*dt/(2*(cc*dt+dx))*(hxbcr(iebc,1)-...
    hxbcf(ibbc,jebc)+hxbcr(2,1)-hxbcf(ibbc-1,jebc));
ezbcr(iebc,jb)=ezbcln(iebc,jb)+(cc*dt-dx)/(cc*dt+dx)*(ezbcl(iebc-1,jb)-...
    ezbcl(iebc,jb))+cc^2*muz*dt/(2*(cc*dt+dx))*(hxbcb(iebc,1)-...
    hxbcl(iebc,je)+hxbcb(iebc-1,1)-hxbcl(iebc-1,je));
% left down corner
ezbcf(1,1)=ezbcfn(2,2)+...
    (cc*dt-2^(1/2)*dx)/(cc*dt+2^(1/2)*dx)*(ezbcf(2,2)-...
    ezbcf(1,1));
% right down corner
ezbcf(ibbc,1)=ezbcbn(ibbc-1,2)+...
    (cc*dt-2^(1/2)*dx)/(cc*dt+2^(1/2)*dx)*(ezbcb(ibbc-1,2)-...
    ezbcf(ibbc,1));
% left up corner
ezbcb(1,jebc)=ezbcbn(2,jebc-1)+...
    (cc*dt-2^(1/2)*dx)/(cc*dt+2^(1/2)*dx)*(ezbcb(2,jebc-1)-...
    ezbcb(1,jebc));
% right up corner
ezbcb(ibbc,jebc)=ezbcbn(ibbc-1,jebc-1)+...
    (cc*dt-2^(1/2)*dx)/(cc*dt+2^(1/2)*dx)*(ezbcb(ibbc-1,jebc-1)-...
    ezbcb(ibbc,jebc));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ezy
% bcl:up down
    ezybcl(1:iebc,1)=caezybcl(1:iebc,1).*ezybcl(1:iebc,1)+...
        cbezybcl(1:iebc,1).*(hxbcl(1:iebc,1)-hxbcf(1:iebc,jebc));
    ezybcl(1:iebc,jb)=caezybcl(1:iebc,jb).*ezybcl(1:iebc,jb)+...
        cbezybcl(1:iebc,jb).*(hxbcb(1:iebc,1)-hxbcl(1:iebc,je));
    %bcr: up down
    ezybcr(1:iebc,1)=caezybcr(1:iebc,1).*ezybcr(1:iebc,1)+...
        cbezybcl(1:iebc,1).*(hxbcr(1:iebc,1)-hxbcf(1:iebc,jebc));
    ezybcr(1:iebc,jb)=caezybcr(1:iebc,jb).*ezybcr(1:iebc,jb)+...
        cbezybcl(1:iebc,jb).*(hxbcb(1:iebc,1)-hxbcr(1:iebc,je));
    % bcf: down
    ezybcf(2:iefbc,1)=ezybcfn(2:iefbc,2)+...
        (cc*dt-dx)/(cc*dt+dx)*(ezybcf(2:iefbc,2)-ezybcf(2:iefbc,1))-...              %注意公式中正负
        cc*cc*muz*dt/2*(cc*dt+dx)*(hybcf(1:iefbc-1,1)-...
        hybcf(2:iefbc,1)+hybcf(1:iefbc-1,2)-hybcf(2:iefbc,2));
    % bcb: up
    ezybcb(2:iefbc,jebc)=ezybcbn(2:iefbc,jebc-1)+...
        (cc*dt-dx)/(cc*dt+dx)*(ezbcb(2:iefbc,jebc-1)-ezbcb(2:iefbc,jebc))-...
        cc*cc*muz*dt/2*(cc*dt+dx)*(hybcb(1:iefbc-1,jebc)-...
   &n ..

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

发帖回复