环形光波导两边射出程序.doc

clear;%tic%Initial(原始的) parameters and other things.W=0.36; %Normalized（规范的） frequency%The following parameters are control parameters.WaveGuide=1; %If this program is for Wave Guide? If so, please specify 1.IsMovie=0; %If you want to play movie, please use 1.IsFigure=1; %If it will plot the figures? If so, please specify 1.WantToSeeEp=1; %Do you want to see the distrubution of Ep? If so, please specify 1.%End of defining control parametersMLatx=11; %11How many Lattice（格子） cell in x direction.MLaty=11; %11How many Lattice cell in y direction.NMlat=21; %21The grid number in each Lattice Cell. %SHOULD BE ODD INTEGER（奇整数）!if mod(NMlat,2)=0 NMlat=NMlat+1;end %Force it to be a odd(奇数) integer!NTx=MLatx*NMlat+1; %It is the number of the Grid along x axis.NTy=MLaty*NMlat+1; %It is the number of the Grid along y axis.if WaveGuide=1 Nrow=5; %4位置问题The row number of columns between the PML boundary and the waveguide.endNPML=12; %12大小就是时间长短问题How many PML layers will be used in our computation.NTimeSteps=2500; %2500Total number of Time StepsMeach=20; %20Define the interval(间隔) for plot figures if IsFigure=1. %This also works for saving intervals. R=0.2; %刚好全反射The radius of dielectric columnsea=11.4; %11.4The dielectric constant of these columns.Zmax=1.2; %0.6The maximum value for z axis when plotting figures.Colormax=0.8; %0.6越大越暗The maximum value for colormap when plotting figures.%Some constantsmu0=4*pi*1.0e-7; %Epsilon Zero, if using Gauss Unit, it equals to 1.e0=8.85*1e-12; %Mu Zero, if using Gauss Unit, it equals to 1.c=1/sqrt(mu0*e0); %The light speed.factor=mu0/e0; %The factor between conductivity(导电率) and permeability（浸透性）. %Permeability=Conductivity*factor, in PML. a=1;%e-6; % 1 The lattice constant.W=W*(2*pi*c/a); %frequency Dx=a/NMlat; %Delta x.Dy=Dx; %Delta y.Dt=1/sqrt(1/(Dx*Dx)+1/(Dy*Dy)/c; %Time interval%tic%In the following partm we define the dielectric constants:Ep=ones(NTx-1,NTy-1)*e0; Ep_cell=ones(NMlat,NMlat)*e0;x=-(NMlat-1)/2*Dx:Dx:(NMlat-1)/2*Dx;X,Y=meshgrid(x);X=X;Y=Y;flag=find(sqrt(X.2+Y.2)R);Ep_cell(flag)=e0*ea;Ep=repmat(Ep_cell,MLatx,MLaty);if WaveGuide=1 eb=1; %The dielectric constant in the waveguide Ep( 1:(MLatx-Nrow)*NMlat, Nrow*NMlat+1:(Nrow+1)*NMlat)=e0*eb; Ep( (MLatx-Nrow-1)*NMlat+1:(MLatx-Nrow)*NMlat, Nrow+1:MLaty*NMlat)=e0*eb; Ep( (MLatx-Nrow-1)*NMlat+1:(MLatx-Nrow)*NMlat,Nrow*NMlat+1:(Nrow+1)*NMlat)=e0*eb;%Ep_cell Ep( (MLatx-Nrow-2)*NMlat+1:(MLatx-Nrow-1)*NMlat,(Nrow+1)*NMlat+1:(Nrow+2)*NMlat)=Ep_cell;%e0*ebend%tocif WantToSeeEp=1 x=0:Dx:(NMlat*MLatx-1)*Dx; x=x-(NMlat*MLatx-1)*Dx/2; y=0:Dy:(NMlat*MLaty-1)*Dy; y=y-(NMlat*MLaty-1)*Dy/2; X,Y=meshgrid(x,y); X=X; Y=Y; surf(X,Y,Ep/e0); shading interp; view(0,90) axis(min(x), max(x),min(y), max(y) axis off; disp(Press any key to continue.); pauseend%以上是没光线进去前的光子晶体图%End of defining the Ep.if IsFigure=1 %Define the X Y coordinate for figures x=0:Dx:(NMlat*MLatx-1+2*NPML)*Dx; x=x-(NMlat*MLatx-1+2*NPML)*Dx/2; y=0:Dy:(NMlat*MLaty-1+2*NPML)*Dy; y=y-(NMlat*MLaty-1+2*NPML)*Dy/2; X,Y=meshgrid(x,y); X=X; Y=Y;end%Define the Ez, Hx, Hy, which are in the inside region.Ez=zeros(NTx-1,NTy-1); Hx=zeros(NTx-1,NTy);Hy=zeros(NTx,NTy-1);%Parameters about PML:n=4; %The order of the polynomial that decribes the conductivity profile.R=1e-10;Delta=NPML*Dx;SigmaMax=-(n+1)*e0*c*log(R)/(Delta*2); NUM=NPML*2:-1:1;Sigmax=SigmaMax*(NUM*Dx/2+Dx/2).(n+1)-(NUM*Dx/2-Dx/2).(n+1)/(Deltan*Dx*(n+1);Sigmay=Sigmax;SigmaBound=SigmaMax*(Dx/2).(n+1)/(Deltan*Dx*(n+1);%Sigmax=SigmaMax*(NUM/6).(n+1); %Another way, the definition.EzxPML1=zeros(NPML,NPML);EzyPML1=zeros(NPML,NPML);HxPML1=zeros(NPML,NPML);HyPML1=zeros(NPML,NPML); %Zone 1Sigmax_z1=repmat(Sigmax(2:2:NPML*2),1,NPML);Sigmax_x1=repmat(Sigmax(2:2:NPML*2),1,NPML);Sigmax_y1=repmat(Sigmax(1:2:NPML*2-1),1,NPML); Sigmay_z1=fliplr(repmat(Sigmax(2:2:NPML*2),NPML,1);Sigmay_x1=fliplr(repmat(Sigmax(1:2:NPML*2-1),NPML,1);Sigmay_y1=fliplr(repmat(Sigmax(2:2:NPML*2),NPML,1); %Zone 1EzxPML2=zeros(NPML,NPML);EzyPML2=zeros(NPML,NPML);HxPML2=zeros(NPML,NPML);HyPML2=zeros(NPML,NPML); %Zone 2Sigmax_z2=flipud(Sigmax_z1);Sigmax_x2=flipud(Sigmax_x1);Sigmax_y2=flipud(Sigmax_y1); Sigmay_z2=Sigmay_z1;Sigmay_x2=Sigmay_x1;Sigmay_y2=Sigmay_y1; %Zone 2EzxPML3=zeros(NPML,NPML);EzyPML3=zeros(NPML,NPML);HxPML3=zeros(NPML,NPML);HyPML3=zeros(NPML,NPML); %Zone 3Sigmax_z3=Sigmax_z1;Sigmax_x3=Sigmax_x1;Sigmax_y3=Sigmax_y1;Sigmay_z3=fliplr(Sigmay_z1);Sigmay_x3=fliplr(Sigmay_x1);Sigmay_y3=fliplr(Sigmay_y1); % Zone 3EzxPML4=zeros(NPML,NPML);EzyPML4=zeros(NPML,NPML);HxPML4=zeros(NPML,NPML);HyPML4=zeros(NPML,NPML); %Zone 4Sigmax_z4=flipud(Sigmax_z1);Sigmax_x4=flipud(Sigmax_x1);Sigmax_y4=flipud(Sigmax_y1); Sigmay_z4=fliplr(Sigmay_z1);Sigmay_x4=fliplr(Sigmay_x1);Sigmay_y4=fliplr(Sigmay_y1); %Zone 4EzxPMLA=zeros(NTx-1,NPML);EzyPMLA=zeros(NTx-1,NPML);HxPMLA=zeros(NTx-1,NPML);HyPMLA=zeros(NTx,NPML); %Zone ASigmay_zA=repmat(Sigmay_z1(1,:),NTx-1,1);Sigmay_xA=repmat(Sigmay_x1(1,:),NTx-1,1);Sigmay_yA=repmat(Sigmay_y1(1,:),NTx,1); %Zone AEzxPMLB=zeros(NTx-1,NPML);EzyPMLB=zeros(NTx-1,NPML);HxPMLB=zeros(NTx-1,NPML);HyPMLB=zeros(NTx,NPML); %Zone BSigmay_zB=fliplr(Sigmay_zA);Sigmay_xB=fliplr(Sigmay_xA);Sigmay_yB=fliplr(Sigmay_yA); %Zone BEzxPMLC=zeros(NPML,NTy-1);EzyPMLC=zeros(NPML,NTy-1);HxPMLC=zeros(NPML,NTy);HyPMLC=zeros(NPML,NTy-1);%Zone CSigmax_zC=repmat(flipud(Sigmay_z1(1,:),1,NTy-1);Sigmax_xC=repmat(flipud(Sigmay_x1(1,:),1,NTy);Sigmax_yC=repmat(flipud(Sigmay_y1(1,:),1,NTy-1); %Zone CEzxPMLD=zeros(NPML,NTy-1);EzyPMLD=zeros(NPML,NTy-1);HxPMLD=zeros(NPML,NTy);HyPMLD=zeros(NPML,NTy-1);%Zone DSigmax_zD=flipud(Sigmax_zC);Sigmax_xD=flipud(Sigmax_xC);Sigmax_yD=flipud(Sigmax_yC); %Zone Dif IsMovie=1 Movie=moviein(NTimeSteps/Meach+1); Mnum=1; EzAll=EzxPML3+EzyPML3,EzxPMLC+EzyPMLC,EzxPML1+EzyPML1; EzxPMLB+EzyPMLB,Ez,EzxPMLA+EzyPMLA; EzxPML4+EzyPML4,EzxPMLD+EzyPMLD,EzxPML2+EzyPML2; HxAll=HxPML3,HxPMLC,HxPML1;HxPMLB,Hx,HxPMLA;HxPML4,HxPMLD,HxPML2; HyAll=HyPML3,HyPMLC,HyPML1;HyPMLB,Hy,HyPMLA;HyPML4,HyPMLD,HyPML2; figure(1) surf(X,Y,real(EzAll) shading interp; axis(-0.5*a*MLatx-Dx*NPML 0.5*a*MLatx+Dx*NPML -0.5*a*MLaty-Dx*NPML 0.5*a*MLaty+Dx*NPML -Zmax Zmax) zlabel(Ez); Movie(:,1)=getframe;endif WaveGuide=1 Npy=NMlat*Nrow+round(NMlat+1)/2); Npx=NMlat*Nrow+round(NMlat+1)/2);else Npx=round(NTx/2); Npy=round(NTy/2);endexpboundary=exp(-SigmaBound*Dt/e0);ticfor m=1:NTimeSteps %tic %The following part is for the inside region. %H components. Hx(:,2:NTy-1)=Hx(:,2:NTy-1)-Dt*(Ez(:,2:NTy-1)-Ez(:,1:NTy-2)/(Dy*mu0); Hy(2:NTx-1,:)=Hy(2:NTx-1,:)+Dt*(Ez(2:NTx-1,:)-Ez(1:NTx-2,:)/(Dx*mu0); Hx(:,NTy)=expboundary*Hx(:,NTy)-(1-expboundary)*. (EzxPMLA(:,1)+EzyPMLA(:,1)-Ez(:,NTy-1)/(SigmaBound*factor*Dy); %Boundary A Hx(:,1)=expboundary*Hx(:,1)-(1-expboundary)*. (Ez(:,1)-EzxPMLB(:,NPML)-EzyPMLB(:,NPML)/(SigmaBound*factor*Dy); %Boundary B Hy(1,:)=expboundary*Hy(1,:)+(1-expboundary)*. (Ez(1,:)-EzxPMLC(NPML,:)-EzyPMLC(NPML,:)/(SigmaBound*factor*Dx); %Boundary C Hy(NTx,:)=expboundary*Hy(NTx,:)+(1-expboundary)*. (EzxPMLD(1,:)+EzyPMLD(1,:)-Ez(NTx-1,:)/(SigmaBound*factor*Dx); %Boundary D %The following part is for ZONE A. H components. HxPMLA(:,1:NPML-1)=exp(-Sigmay_xA(:,1:NPML-1)*Dt/e0).*HxPMLA(:,1:NPML-1)-. (1-exp(-Sigmay_xA(:,1:NPML-1)*Dt/e0)./(Sigmay_xA(:,1:NPML-1)*factor*Dy).*. (EzxPMLA(:,2:NPML)+EzyPMLA(:,2:NPML)-EzxPMLA(:,1:NPML-1)-EzyPMLA(:,1:NPML-1); HyPMLA(2:NTx-1,:)=HyPMLA(2:NTx-1,:)+. Dt*(EzxPMLA(2:NTx-1,:)+EzyPMLA(2:NTx-1,:)-EzxPMLA(1:NTx-2,:)-EzyPMLA(1:NTx-2,:)/(mu0*Dx); HyPMLA(1,:)=expboundary*HyPMLA(1,:)+(1-expboundary)*. (EzxPMLA(1,:)+EzyPMLA(1,:)-EzxPML1(NPML,:)-EzyPML1(NPML,:)/(SigmaBound*factor*Dx); %Boundary Left HyPMLA(NTx,:)=expboundary*HyPMLA(NTx,:)+(1-expboundary)*. (EzxPML2(1,:)+EzyPML2(1,:)-EzxPMLA(NTx-1,:)-EzyPMLA(NTx-1,:)/(SigmaBound*factor*Dx); %Boundary Right HxPMLA(:,NPML)=HxPMLA(:,NPML-1); %Boundary Upper; %The following part is for ZONE B. %H components. HxPMLB(:,2:NPML)=exp(-Sigmay_xB(:,2:NPML)*Dt/e0).*HxPMLB(:,2:NPML)-. (1-exp(-Sigmay_xB(:,2:NPML)*Dt/e0)./(Sigmay_xB(:,2:NPML)*factor*Dy).*. (EzxPMLB(:,2:NPML)+EzyPMLB(:,2:NPML)-EzxPMLB(:,1:NPML-1)-EzyPMLB(:,1:NPML-1); HyPMLB(2:NTx-1,:)=HyPMLB(2:NTx-1,:)+. Dt*(EzxPMLB(2:NTx-1,:)+EzyPMLB(2:NTx-1,:)-EzxPMLB(1:NTx-2,:)-EzyPMLB(1:NTx-2,:)/(mu0*Dx); HyPMLB(1,:)=expboundary*HyPMLB(1,:)+(1-expboundary)*. (EzxPMLB(1,:)+EzyPMLB(1,:)-EzxPML3(NPML,:)-EzyPML3(NPML,:)/(SigmaBound*factor*Dx); %Boundary Left HyPMLB(NTx,:)=expboundary*HyPMLB(NTx,:)+(1-expboundary)*. (EzxPML4(1,:)+EzyPML4(1,:)-EzxPMLB(NTx-1,:)-EzyPMLB(NTx-1,:)/(SigmaBound*factor*Dx); %Boundary Right HxPMLB(:,1)=HxPMLB(:,2); %Boundary Bottom; %The following part is for ZONE C. %H components. HxPMLC(:,2:NTy-1)=HxPMLC(:,2:NTy-1)-. Dt*(EzxPMLC(:,2:NTy-1)+EzyPMLC(:,2:NTy-1)-EzxPMLC(:,1:NTy-2)-EzyPMLC(:,1:NTy-2)/(mu0*Dy); HyPMLC(2:NPML,:)=exp(-Sigmax_yC(2:NPML,:)*Dt/e0).*HyPMLC(2:NPML,:)+. (1-exp(-Sigmax_yC(2:NPML,:)*Dt/e0)./(Sigmax_yC(2:NPML,:)*factor*Dx).*. (EzxPMLC(2:NPML,:)+EzyPMLC(2:NPML,:)-EzxPMLC(1:NPML-1,:)-EzyPMLC(1:NPML-1,:); HxPMLC(:,1)=expboundary*HxPMLC(:,1)-(1-expboundary)*. (EzxPMLC(:,1)+EzyPMLC(:,1)-EzxPML3(:,NPML)-EzyPML3(:,NPML)/(SigmaBound*factor*Dy); %Boundary Down HxPMLC(:,NTy)=expboundary*HxPMLC(:,NTy)-(1-expboundary)*. (EzxPML1(:,1)+EzyPML1(:,1)-EzxPMLC(:,NTy-1)-EzyPMLC(:,NTy-1)/(SigmaBound*factor*Dy); %Boundary Upper HyPMLC(1,:)=HyPMLC(2,:); %Boundary Left; %The following part is for ZONE D. %H components. HxPMLD(:,2:NTy-1)=HxPMLD(:,2:NTy-1)-. Dt*(EzxPMLD(:,2:NTy-1)+EzyPMLD(:,2:NTy-1)-EzxPMLD(:,1:NTy-2)-EzyPMLD(:,1:NTy-2)/(mu0*Dy); HyPMLD(1:NPML-1,:)=exp(-Sigmax_yD(1:NPML-1,:)*Dt/e0).*HyPMLD(1:NPML-1,:)+. (1-exp(-Sigmax_yD(1:NPML-1,:)*Dt/e0)./(Sigmax_yD(1:NPML-1,:)*factor*Dx).*. (EzxPMLD(2:NPML,:)+EzyPMLD(2:NPML,:)-EzxPMLD(1:NPML-1,:)-EzyPMLD(1:NPML-1,:); HxPMLD(:,1)=expboundary*HxPMLD(:,1)-(1-expboundary)*. (EzxPMLD(:,1)+EzyPMLD(:,1)-EzxPML4(:,NPML)-EzyPML4(:,NPML)/(SigmaBound*factor*Dy); %Boundary Down HxPMLD(:,NTy)=expboundary*HxPMLD(:,NTy)-(1-expboundary)*. (EzxPML2(:,1)+EzyPML2(:,1)-EzxPMLD(:,NTy-1)-EzyPMLD(:,NTy-1)/(SigmaBound*factor*Dy); %Boundary Upper HyPMLD(NPML,:)=HyPMLD(NPML-1,:); %Boundary Right; %The following part is for ZONE 1. %H components. HxPML1(:,1:NPML-1)=exp(-Sigmay_x1(:,1:NPML-1)*Dt/e0).*HxPML1(:,1:NPML-1)-. (1-exp(-Sigmay_x1(:,1:NPML-1)*Dt/e0)./(Sigmay_x1(:,1:NPML-1)*factor*Dy).*. (EzxPML1(:,2:NPML)+EzyPML1(:,2:NPML)-EzxPML1(:,1:NPML-1)-EzyPML1(:,1:NPML-1); HxPML1(:,NPML)=HxPML1(:,NPML-1); HyPML1(2:NPML,:)=exp(-Sigmax_y1(2:NPML,:)*Dt/e0).*HyPML1(2:NPML,:)+. (1-exp(-Sigmax_y1(2:NPML,:)*Dt/e0)./(Sigmax_y1(2:NPML,:)*factor*Dx).*. (EzxPML1(2:NPML,:)+EzyPML1(2:NPML,:)-EzxPML1(1:NPML-1,:)-EzyPML1(1:NPML-1,:); HyPML1(1,:)=HyPML1(2,:); %The following part is for ZONE 2. %H components. HxPML2(:,1:NPML-1)=exp(-Sigmay_x2(:,1:NPML-1)*Dt/e0).*HxPML2(:,1:NPML-1)-. (1-exp(-Sigmay_x2(:,1:NPML-1)*Dt/e0)./(Sigmay_x2(:,1:NPML-1)*factor*Dy).*. (EzxPML2(:,2:NPML)+EzyPML2(:,2:NPML)-EzxPML2(:,1:NPML-1)-EzyPML2(:,1:NPML-1); HxPML2(:,NPML)=HxPML2(:,NPML-1); HyPML2(1:NPML-1,:)=exp(-Sigmax_y2(1:NPML-1,:)*Dt/e0).*HyPML2(1:NPML-1,:)+. (1-exp(-Sigmax_y2(1:NPML-1,:)*Dt/e0)./(Sigmax_y2(1:NPML-1,:)*factor*Dx).*. (EzxPML2(2:NPML,:)+EzyPML2(2:NPML,:)-EzxPML2(1:NPML-1,:)-EzyPML2(1:NPML-1,:); HyPML2(NPML,:)=HyPML2(NPML-1,:); %Boundary Right; %The following part is for ZONE 3. %H components. HyPML3(2:NPML,:)=exp(-Sigmax_y3(2:NPML,:)*Dt/e0).*HyPML3(2:NPML,:)+. (1-exp(-Sigmax_y3(2:NPML,:)*Dt/e0)./(Sigmax_y3(2:NPML,:)*factor*Dx).*. (EzxPML3(2:NPML,:)+EzyPML3(2:NPML,:)-EzxPML3(1:NPML-1,:)-EzyPML3(1:NPML-1,:); HyPML3(1,:)=HyPML3(2,:); %Boundary Left; HxPML3(:,2:NPML)=exp(-Sigmay_x3(:,2:NPML)*Dt/e0).*HxPML3(:,2:NPML)-. (1-exp(-Sigmay_x3(:,2:NPML)*Dt/e0)./(Sigmay_x3(:,2:NPML)*factor*Dy).*. (EzxPML3(:,2:NPML)+EzyPML3(:,2:NPML)-EzxPML3(:,1:NPML-1)-EzyPML3(:,1:NPML-1); HxPML3(:,1)=HxPML3(:,2); %Boundary Bottom; %The following part is for ZONE 4. %H components. HxPML4(:,2:NPML)=exp(-Sigmay_x4(:,2:NPML)*Dt/e0).*HxPML4(:,2:NPML)-. (1-exp(-Sigmay_x4(:,2:NPML)*Dt/e0)./(Sigmay_x4(:,2:NPML)*factor*Dy).*. (EzxPML4(:,2:NPML)+EzyPML4(:,2:NPML)-EzxPML4(:,1:NPML-1)-EzyPML4(:,1:NPML-1); HxPML4(:,1)=HxPML4(:,2); %Boundary Bottom; HyPML4(1:NPML-1,:)=exp(-Sigmax_y4(1:NPML-1,:)*Dt/e0).*HyPML4(1:NPML-1,:)+. (1-exp(-Sigmax_y4(1:NPML-1,:)*Dt/e0)./(Sigmax_y4(1:NPML-1,:)*factor*Dx).*. (EzxPML4(2:NPML,:)+EzyPML4(2:NPML,:)-EzxPML4(1:NPML-1,:)-EzyPML4(1:NPML-1,:); HyPML4(NPML,:)=HyPML4(NPML-1,:); %Boundary Right; %The following part is for inside region. %Ez component. Ez=Ez+Dt*(Hy(2:NTx,1:NTy-1)-Hy(1:NTx-1,1:NTy-1)/Dx-. (Hx(1:NTx-1,2:NTy)-Hx(1:NTx-1,1:NTy-1)/Dy)./Ep; %Inside region. %The following part is for source Ez(1,Npy-(NMlat-1)/2:Npy+(NMlat-1)/2)=Ez(1,Npy-(NMlat-1)/2:Npy+(NMlat-1)/2)+sin(W*m*Dt);%*exp(-(m*W*Dt-3)2); %End of Source %The following part is for ZONE A. %Ez component. EzxPMLA=EzxPMLA+Dt*(HyPMLA(2:NTx,:)-HyPMLA(1:NTx-1,:)/(e0*Dx); EzyPMLA(:,2:NPML)=exp(-Sigmay_zA(:,2:NPML)*Dt/e0).*EzyPMLA(:,2:NPML)-. (1-exp(-Sigmay_zA(:,2:NPML)*Dt/e0)./(Sigmay_zA(:,2:NPML)*Dy).*. (HxPMLA(:,2:NPML)-HxPMLA(:,1:NPML-1); EzyPMLA(:,1)=exp(-Sigmay_zA(:,1)*Dt/e0).*EzyPMLA(:,1)-. (1-exp(-Sigmay_zA(:,1)*Dt/e0)./(Sigmay_zA(:,1)*Dy).*. (HxPMLA(:,1)-Hx(:,NTy); %The following part is for ZONE B. %Ez component. EzxPMLB=EzxPMLB+Dt*(HyPMLB(2:NTx,:)-HyPMLB(1:NTx-1,:)/(e0*Dx); EzyPMLB(:,1:NPML-1)=exp(-Sigmay_zB(:,1:NPML-1)*Dt/e0).*EzyPMLB(:,1:NPML-1)-. (1-exp(-Sigmay_zB(:,1:NPML-1)*Dt/e0)./(Sigmay_zB(:,1:NPML-1)*Dy).*. (HxPMLB(:,2:NPML)-HxPMLB(:,1:NPML-1); EzyPMLB(:,NPML)=exp(-Sigmay_zB(:,NPML)*Dt/e0).*EzyPMLB(:,NPML)-. (1-exp(-Sigmay_zB(:,NPML)*Dt/e0)./(Sigmay_zB(:,NPML)*Dy).*. (Hx(:,1)-HxPMLB(:,NPML); %The following part is for ZONE C. %Ez component. EzxPMLC(1:NPML-1,:)=exp(-Sigmax_zC(1:NPML-1,:)*Dt/e0).*EzxPMLC(1:NPML-1,:)+. (1-exp(-Sigmax_zC(1:NPML-1,:)*Dt/e0)./(Sigmax_zC(1:NPML-1,:)*Dx).*. (HyPMLC(2:NPML,:)-HyPMLC(1:NPML-1,:); EzxPMLC(NPML,:)=exp(-Sigmax_zC(NPML,:)*Dt/e0).*EzxPMLC(NPML,:)+. (1-exp(-Sigmax_zC(NPML,:)*Dt/e0)./(Sigmax_zC(NPML,:)*Dx).*. (Hy(1,:)-HyPMLC(NPML,:); EzyPMLC=EzyPMLC-Dt*(HxPMLC(:,2:NTy)-HxPMLC(:,1:NTy-1)/(e0*Dy); %The following part is for ZONE D. %Ez component. EzxPMLD(2:NPML,:)=exp(-Sigmax_zD(2:NPML,:)*Dt/e0).*EzxPMLD(2:NPML,:)+. (1-exp(-Sigmax_zD(2:NPML,:)*Dt/e0)./(Sigmax_zD(2:NPML,:)*Dx).*. (HyPMLD(2:NPML,:)-HyPMLD(1:NPML-1,:); EzxPMLD(1,:)=exp(-Sigmax_zD(1,:)*Dt/e0).*EzxPMLD(1,:)+. (1-exp(-Sigmax_zD(1,:)*Dt/e0)./(Sigmax_zD(1,:)*Dx).*. (HyPMLD(1,:)-Hy(NTx,:); EzyPMLD=EzyPMLD-Dt*(HxPMLD(:,2:NTy)-HxPMLD(:,1:NTy-1)/(e0*Dy); %The following part is for ZONE 1. %Ez component. EzxPML1(1:NPML-1,:)=exp(-Sigmax_z1(1:NPML-1,:)*Dt/e0).*EzxPML1(1:NPML-1,:)+. (1-exp(-Sigmax_z1(1:NPML-1,:)*Dt/e0)./(Sigmax_z1(1:NPML-1,:)*Dx).*. (HyPML1(2:NPML,:)-HyPML1(1:NPML-1,:); EzxPML1(NPML,:)=exp(-Sigmax_z1(NPML,:)*Dt/e0).*EzxPML1(NPML,:)+. (1-exp(-Sigmax_z1(NPML,:)*Dt/e0)./(Sigmax_z1(NPML,:)*Dx).*. (HyPMLA(1,:)-HyPML1(NPML,:); EzyPML1(:,2:NPML)=exp(-Sigmay_z1(:,2:NPML)*Dt/e0).*EzyPML1(:,2:NPML)-. (1-exp(-Sigmay_z1(:,2:NPML)*Dt/e0)./(Sigmay_z1(:,2:NPML)*Dy).*. (HxPML1(:,2:NPML)-HxPML1(:,1:NPML-1); EzyPML1(:,1)=exp(-Sigmay_z1(:,1)*Dt/e0).*EzyPML1(:,1)-. (1-exp(-Sigmay_z1(:,1)*Dt/e0)./(Sigmay_z1(:,1)*Dy).*. (HxPML1(:,1)-HxPMLC(:,NTy); %The following part is for ZONE 2. %Ez component. EzxPML2(2:NPML,:)=exp(-Sigmax_z2(2:NPML,:)*Dt/e0).*EzxPML2(2:NPML,:)+. (1-exp(-Sigmax_z2(2:NPML,:)*Dt/e0)./(Sigmax_z2(2:NPML,:)*Dx).*. (HyPML2(2:NPML,:)-HyPML2(1:NPML-1,:); EzxPML2(1,:)=exp(-Sigmax_z2(1,:)*Dt/e0).*EzxPML2(1,:)+. (1-exp(-Sigmax_z2(1,:)*Dt/e0)./(Sigmax_z2(1,:)*Dx).*. (HyPML2(1,:)-HyPMLA(NTx,:); EzyPML2(:,2:NPML)=exp(-Sigmay_z2(:,2:NPML)*Dt/e0).*EzyPML2(:,2:NPML)-. (1-exp(-Sigmay_z2(:,2:NPML)*Dt/e0)./(Sigmay_z2(:,2:NPML)*Dy).*. (HxPML2(:,2:NPML)-HxPML2(:,1:NPML-1); EzyPML2(:,1)=exp(-Sigmay_z2(:,1)*Dt/e0).*EzyPML2(:,1)-. (1-exp(-Sigmay_z2(:,1)*Dt/e0)./(Sigmay_z2(:,1)*Dy).*. (HxPML2(:,1)-HxPMLD(:,NTy);