initializepml.m

外国人开发的电磁时域有限差分方法工具包 Electromagnetic Finite-Difference Time-Domain (EmFDTD) is a basic two-dimensio

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