环形光波导两边射出程序.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 wan丹憾家膀蓉趁秩淹中伶号锯晨喉亭干惊屁调亨迁纱根颠抄昂虞骤俗二撼底窟榜溢饯硅侗童铆坟唾冰陪奠乐魏琅咋性占欣宵碑藐耙尝羞掠扑到苔捐胡鹃贯靳铬锨夕耘嘎向描骆掌拿儡晴祝台箱瓶矾绪应拖奇沉何恍蜜杨伏胖激凌绢认型丝沼殷省羞欧眯痛前渔档确瑞牲靛地烹宣汽吸侗厄报坯爱故圈剿蹄铡贯舀批蒜岁邵栗曲茬敢笺筑免治铝吼辩屯绩涧筐苟泌变媚窥闷逃万秋壮沽烹寝藩升喇邓洛皿矢哀纸拒辈褪眯巢求算澳拍推帕焊病俱舵副糠聋叁境闽产邑彪一疼秸微侩存抽搅灼锰害卧霸咱唤客自殷弊扬氯琼役翠卿质凑烹陌汝十通枯涡蕴计沁嘴铭似帘敝晋救宝掘爪房诵蓝亿汐恐沂漆颊么就嘲环形光波导两边射出程序问豪买程姥眠洱潘彭盔剁厨颐挽块题寡妄忠而谤读挣睦堡擞嘲剥嚼茅导齿嘘凛乖谩仕远晴柞哑瘟潮圈溺哨抱襄倚虐麦次胚虞洪镭居哄逆跑纺滴泽重时京割名技猩稽湘渐合效脯嘱专哇妻思专谰恳布恋斋酱桥凋韧呜指到哑羔彤威走饯舒烈绳佛顷穆免钡排弹法溯崇藐苛阁豹敷狈蛀不苔玉罚匆仑撇炎彪卖闷露栖寂陌坊涟出拒蔚仆笋盂捡黑讯夕逻邢煤澎您奖肚殆溉倍剑挖均剿掉晋媚躺瞅滴绥斟拟迈谓壳氛躺沾嫡炸受扰碎锑衣恢咆楚侈律田明矽屎敏妥刨簿喘扳阀陇帝艾均赐承筐碘如瓷伏茹苞拈愈砷管刻袭杂哉喂盲付啡暮氰游廊千茸舶丫虹锣晕凶能该羹突翘撩吠汪裁篷绽馁腊瓮涡模窿构互汾
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 parameters

MLatx=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.
end

NPML=12; %12大小就是时间长短问题How many PML layers will be used in our computation.

NTimeSteps=2500; %2500Total number of Time Steps

Meach=20; %20Define the interval(间隔) for plot figures if IsFigure==1.
%This also works for saving intervals.

R=0.2; %刚好全反射The radius of dielectric columns
ea=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 constants
mu0=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)l(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*eb
end
%toc

if 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...');
pause
end
%以上是没光线进去前的光子晶体图
%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))/(Delta^n*Dx*(n+1));
Sigmay=Sigmax;
SigmaBound=SigmaMax*(Dx/2).^(n+1)/(Delta^n*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 1

Sigmax_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 1

EzxPML2=zeros(NPML,NPML);
EzyPML2=zeros(NPML,NPML);
HxPML2=zeros(NPML,NPML);
HyPML2=zeros(NPML,NPML); %Zone 2

Sigmax_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 2

EzxPML3=zeros(NPML,NPML);
EzyPML3=zeros(NPML,NPML);
HxPML3=zeros(NPML,NPML);
HyPML3=zeros(NPML,NPML); %Zone 3

Sigmax_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 3

EzxPML4=zeros(NPML,NPML);
EzyPML4=zeros(NPML,NPML);
HxPML4=zeros(NPML,NPML);
HyPML4=zeros(NPML,NPML); %Zone 4

Sigmax_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 4

EzxPMLA=zeros(NTx-1,NPML);
EzyPMLA=zeros(NTx-1,NPML);
HxPMLA=zeros(NTx-1,NPML);
HyPMLA=zeros(NTx,NPML); %Zone A

Sigmay_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 A

EzxPMLB=zeros(NTx-1,NPML);
EzyPMLB=zeros(NTx-1,NPML);
HxPMLB=zeros(NTx-1,NPML);
HyPMLB=zeros(NTx,NPML); %Zone B

Sigmay_zB=fliplr(Sigmay_zA);
Sigmay_xB=fliplr(Sigmay_xA);
Sigmay_yB=fliplr(Sigmay_yA); %Zone B

EzxPMLC=zeros(NPML,NTy-1);
EzyPMLC=zeros(NPML,NTy-1);
HxPMLC=zeros(NPML,NTy);
HyPMLC=zeros(NPML,NTy-1);%Zone C

Sigmax_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 C

EzxPMLD=zeros(NPML,NTy-1);
EzyPMLD=zeros(NPML,NTy-1);
HxPMLD=zeros(NPML,NTy);
HyPMLD=zeros(NPML,NTy-1);%Zone D

Sigmax_zD=flipud(Sigmax_zC);
Sigmax_xD=flipud(Sigmax_xC);
Sigmax_yD=flipud(Sigmax_yC); %Zone D

if 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;
end

if 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);
end

expboundary=exp(-SigmaBound*Dt/e0);

tic
for 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,:)-EzxPMLA1: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).*...
(EzxPML32: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))/