www.gusucode.com > FDTD计算二维光子带隙的程序 > FDTD计算二维光子带隙的程序/247787/InitializePML.m
function InitializePML %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Electromagnetic Finite-Difference Time-Domain % % Version 1.20, Release 1 % % % % (C) Copyright 2005 % % Sharif University of Technology % % School of Electrical Engineering % % All Rights Reserved % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% global xCnt yCnt a X Y xPMLCnt yPMLCnt LX LY global BoundaryType dT global Mu Epsilon global f1 f2 g1 g2 h1 h2 k1 k2 a=1; LX=1; if BoundaryType<=3 X=xCnt; Y=yCnt; LY=1*Y/X; xPMLCnt=0; yPMLCnt=0; eSigmax=zeros(X,Y); eSigmay=zeros(X,Y); mSigmax=zeros(X,Y); mSigmay=zeros(X,Y); f=dT/2*mSigmay./Mu; f1=(1-f)./(1+f); f2=dT./Mu./(1+f); g=dT/2*mSigmax./Mu; g1=(1-g)./(1+g); g2=dT./Mu./(1+g); f=dT/2*eSigmax./Epsilon; h1=(1-f)./(1+f); h2=dT./Epsilon./(1+f); g=dT/2*eSigmay./Epsilon; k1=(1-g)./(1+g); k2=dT./Epsilon./(1+g); return end PMLCnt=20; PMLExp=4; PMLeSigma=20; X=PMLCnt+xCnt+PMLCnt; Y=PMLCnt+yCnt+PMLCnt; LY=1*Y/X; EpTemp=zeros(X,Y); MuTemp=zeros(X,Y); EpTemp(PMLCnt+1:X-PMLCnt,PMLCnt+1:Y-PMLCnt)=Epsilon; MuTemp(PMLCnt+1:X-PMLCnt,PMLCnt+1:Y-PMLCnt)=Mu; eSigmax=zeros(X,Y); eSigmay=zeros(X,Y); mSigmax=zeros(X,Y); mSigmay=zeros(X,Y); EpCL=Epsilon(:,1); MuCL=Mu(:,1); EpCR=Epsilon(:,yCnt); MuCR=Mu(:,yCnt); EpRU=[Row(Epsilon(1,1),PMLCnt) Epsilon(1,:) Row(Epsilon(1,yCnt),PMLCnt)]; EpRD=[Row(Epsilon(xCnt,1),PMLCnt) Epsilon(xCnt,:) Row(Epsilon(xCnt,yCnt),PMLCnt)]; MuRU=[Row(Mu(1,1),PMLCnt) Mu(1,:) Row(Mu(1,yCnt),PMLCnt)]; MuRD=[Row(Mu(xCnt,1),PMLCnt) Mu(xCnt,:) Row(Mu(xCnt,yCnt),PMLCnt)]; R=Row(PMLeSigma,Y); C=Column(PMLeSigma,X); for m=1:PMLCnt EpTemp(PMLCnt+1:X-PMLCnt,m)=EpCL; MuTemp(PMLCnt+1:X-PMLCnt,m)=MuCL; EpTemp(PMLCnt+1:X-PMLCnt,Y+1-m)=EpCR; MuTemp(PMLCnt+1:X-PMLCnt,Y+1-m)=MuCR; EpTemp(m,:)=EpRU; MuTemp(m,:)=MuRU; EpTemp(X+1-m,:)=EpRD; MuTemp(X+1-m,:)=MuRD; eSigmax(PMLCnt-m+1,:)=(m/PMLCnt)^PMLExp*R; eSigmax(X-PMLCnt+m,:)=(m/PMLCnt)^PMLExp*R; eSigmay(:,m)=((PMLCnt+1-m)/PMLCnt)^PMLExp*C; eSigmay(:,Y-m+1)=((PMLCnt+1-m)/PMLCnt)^PMLExp*C; end Epsilon=EpTemp; Mu=MuTemp; mSigmax=eSigmax.*Mu./Epsilon; mSigmay=eSigmay.*Mu./Epsilon; f=dT/2*mSigmay./Mu; f1=(1-f)./(1+f); f2=dT./Mu./(1+f); g=dT/2*mSigmax./Mu; g1=(1-g)./(1+g); g2=dT./Mu./(1+g); f=dT/2*eSigmax./Epsilon; h1=(1-f)./(1+f); h2=dT./Epsilon./(1+f); g=dT/2*eSigmay./Epsilon; k1=(1-g)./(1+g); k2=dT./Epsilon./(1+g); if BoundaryType==4 Epsilon=Epsilon(PMLCnt+1:X-PMLCnt,:); Mu=Mu(PMLCnt+1:X-PMLCnt,:); f1=f1(PMLCnt+1:X-PMLCnt,:); f2=f2(PMLCnt+1:X-PMLCnt,:); g1=g1(PMLCnt+1:X-PMLCnt,:); g2=g2(PMLCnt+1:X-PMLCnt,:); h1=h1(PMLCnt+1:X-PMLCnt,:); h2=h2(PMLCnt+1:X-PMLCnt,:); k1=k1(PMLCnt+1:X-PMLCnt,:); k2=k2(PMLCnt+1:X-PMLCnt,:); X=xCnt; xPMLCnt=0; yPMLCnt=PMLCnt; elseif BoundaryType==5 Epsilon=Epsilon(:,PMLCnt+1:Y-PMLCnt); Mu=Mu(:,PMLCnt+1:Y-PMLCnt); f1=f1(:,PMLCnt+1:Y-PMLCnt); f2=f2(:,PMLCnt+1:Y-PMLCnt); g1=g1(:,PMLCnt+1:Y-PMLCnt); g2=g2(:,PMLCnt+1:Y-PMLCnt); h1=h1(:,PMLCnt+1:Y-PMLCnt); h2=h2(:,PMLCnt+1:Y-PMLCnt); k1=k1(:,PMLCnt+1:Y-PMLCnt); k2=k2(:,PMLCnt+1:Y-PMLCnt); Y=yCnt; xPMLCnt=PMLCnt; yPMLCnt=0; else xPMLCnt=PMLCnt; yPMLCnt=PMLCnt; end