mtlab pml.txt

摩尔边界的源码在计算电磁学中有很大应用

%% Simulation of a pulse hitting UPML
%% To test the UPML boundary conditions
%% 2D TM wave

% Diogram of the areas
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% A %%         3       %% D %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   %%                 %%   %%
%%   %%                 %%   %%
%%   %%                 %%   %%
%%   %%                 %%   %%
%% 1 %%       FDTD      %% 2 %%
%%   %%                 %%   %%
%%   %%                 %%   %%
%%   %%                 %%   %%
%%   %%                 %%   %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% B %%         4       %% C %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear,clc

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters to define the FDTD lattice
KX=100;    
KY=100;
kcx=KY/2;
kcy=KX/2;
Npml=10;
Lx=KX-2*Npml;
Ly=KY-2*Npml;
Lpml=Npml;
Rpml=KX-Npml+1;
Bpml=Npml;
Tpml=KY-Npml+1;
% Parameters for the Gaussian pulse and the FDTD cells
t0=40.0;
spread=10;
c0=3e8;
dx=2;
dt=0.5*dx/c0;
freq=c0/50/dx;
ep_r=5;  % Relative epsilon of the medium around UPML
mu_r=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters control the displayer
N_Steps=20;      % Increment for NSTEPS
NSTEPS=N_Steps;  % Temporary variable for time steps
T_Des=800;       % Destination value of time steps
n_start=1;        % Record the start value of n
n=n_start;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Want to see the contour plot of Ez,please specify it to 1
SeeContour=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Matrixs for the E and H fields
Ez=zeros(KY,KX);
Hx=zeros(KY-1,KX);
Hy=zeros(KY,KX-1);

% Matrics for the Edge areas
DzA=zeros(Npml,Npml);
DzB=zeros(Npml,Npml);
DzC=zeros(Npml,Npml);
DzD=zeros(Npml,Npml);

BxA=zeros(Npml-1,Npml);
BxB=zeros(Npml-1,Npml);
BxC=zeros(Npml-1,Npml);
BxD=zeros(Npml-1,Npml);

ByA=zeros(Npml,Npml-1);
ByB=zeros(Npml,Npml-1);
ByC=zeros(Npml,Npml-1);
ByD=zeros(Npml,Npml-1);

% Matrics for the Plane areas
BX1=zeros(Ly+1,Npml); % 由于电场算入了棱边区,导致在
BX2=zeros(Ly+1,Npml); % 平面区磁场比电场多一个抽样点
BY3=zeros(Npml,Lx+1);
BY4=zeros(Npml,Lx+1);

% Matrixs to record the previous value of the fields
DzA_prev=DzA;
DzB_prev=DzB;
DzC_prev=DzC;
DzD_prev=DzD;

BxA_prev=BxA;
BxB_prev=BxB;
BxC_prev=BxC;
BxD_prev=BxD;

ByA_prev=ByA;
ByB_prev=ByB;
ByC_prev=ByC;
ByD_prev=ByD;

BX1_prev=BX1;
BX2_prev=BX2;

BY3_prev=BY3;
BY4_prev=BY4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% UPML Medium parameters
epsz=8.85419e-12;
sigma_M=1/30/pi/dx/sqrt(ep_r);
% Electric Conductivity of the medium
% the parameters are for the Zone A area
sigma_xE=sigma_M*((Npml-1:-1:0)/Npml).^4;      % Npml elements
sigma_xH=sigma_M*((Npml-1.5:-1:0.5)/Npml).^4;  % Npml-1 elements
eaf_xE=sigma_xE*dt/2/epsz;
eaf_yE=rot90(eaf_xE);
eaf_xH=sigma_xH*dt/2/epsz;
eaf_yH=rot90(eaf_xH);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters for the Edge areas
% CDz
CDz_A=(1-eaf_xE)./(1+eaf_xE);
CDz_A=repmat(CDz_A,Npml,1);
CDz_B=flipud(CDz_A);
CDz_C=fliplr(CDz_B);
CDz_D=flipud(CDz_C);
% CHxDz
CHxDz_A=0.5./(1+eaf_xE);
CHxDz_A=repmat(CHxDz_A,Npml,1);
CHxDz_B=flipud(CHxDz_A);
CHxDz_C=fliplr(CHxDz_B);
CHxDz_D=flipud(CHxDz_C);
% CEz
CEz_A=(1-eaf_yE)./(1+eaf_yE);
CEz_A=repmat(CEz_A,1,Npml);
CEz_B=flipud(CEz_A);
CEz_C=fliplr(CEz_B);
CEz_D=flipud(CEz_C);
% CDzEz
CDzEz_A=1/ep_r./(1+eaf_yE);
CDzEz_A=repmat(CDzEz_A,1,Npml);
CDzEz_B=flipud(CDzEz_A);
CDzEz_C=fliplr(CDzEz_B);
CDzEz_D=flipud(CDzEz_C);
% CBx
CBx_A=(1-eaf_yH)./(1+eaf_yH);
CBx_A=repmat(CBx_A,1,Npml);
CBx_B=flipud(CBx_A);
CBx_C=fliplr(CBx_B);
CBx_D=flipud(CBx_C);
% CEzBx
CEzBx_A=0.5./(1+eaf_yH);
CEzBx_A=repmat(CEzBx_A,1,Npml);
CEzBx_B=flipud(CEzBx_A);
CEzBx_C=fliplr(CEzBx_B);
CEzBx_D=flipud(CEzBx_C);
% CBxHx_New
CBxHx_New_A=(1+eaf_xE)/mu_r;
CBxHx_New_A=repmat(CBxHx_New_A,Npml-1,1);
CBxHx_New_B=flipud(CBxHx_New_A);
CBxHx_New_C=fliplr(CBxHx_New_B);
CBxHx_New_D=flipud(CBxHx_New_C);
% CBxHx_Old
CBxHx_Old_A=(1-eaf_xE)/mu_r;
CBxHx_Old_A=repmat(CBxHx_Old_A,Npml-1,1);
CBxHx_Old_B=flipud(CBxHx_Old_A);
CBxHx_Old_C=fliplr(CBxHx_Old_B);
CBxHx_Old_D=flipud(CBxHx_Old_C);
% CBy
CBy_A=(1-eaf_xH)./(1+eaf_xH);
CBy_A=repmat(CBy_A,Npml,1);
CBy_B=flipud(CBy_A);
CBy_C=fliplr(CBy_B);
CBy_D=flipud(CBy_C);
% CEzBy
CEzBy_A=0.5./(1+eaf_xH);
CEzBy_A=repmat(CEzBy_A,Npml,1);
CEzBy_B=flipud(CEzBy_A);
CEzBy_C=fliplr(CEzBy_B);
CEzBy_D=flipud(CEzBy_C);
% CByHy_New
CByHy_New_A=(1+eaf_yE)/mu_r;
CByHy_New_A=repmat(CByHy_New_A,1,Npml-1);
CByHy_New_B=flipud(CByHy_New_A);
CByHy_New_C=fliplr(CByHy_New_B);
CByHy_New_D=flipud(CByHy_New_C);
% CByHy_Old
CByHy_Old_A=(1-eaf_yE)/mu_r;
CByHy_Old_A=repmat(CByHy_Old_A,1,Npml-1);
CByHy_Old_B=flipud(CByHy_Old_A);
CByHy_Old_C=fliplr(CByHy_Old_B);
CByHy_Old_D=flipud(CByHy_Old_C);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters for the Plane areas
%%%%%%%%%%%%%%%   X direction   %%%%%%%%%%%%%%%%%%%%
% CEz_Px
CEz_Px_I=(1-eaf_xE)./(1+eaf_xE);
CEz_Px_I=repmat(CEz_Px_I,Ly,1);
CEz_Px_II=fliplr(CEz_Px_I);
% CHxEz_Px
CHxEz_Px_I=0.5/ep_r./(1+eaf_xE);
CHxEz_Px_I=repmat(CHxEz_Px_I,Ly,1);
CHxEz_Px_II=fliplr(CHxEz_Px_I);
% CEzBx_Px
CEzBx_Px=0.5/mu_r;
% CBxHx_New_Px
CBxHx_New_Px_I=1+eaf_xE;
CBxHx_New_Px_I=repmat(CBxHx_New_Px_I,Ly+1,1);
CBxHx_New_Px_II=fliplr(CBxHx_New_Px_I);
% CBxHx_Old_Px
CBxHx_Old_Px_I=1-eaf_xE;
CBxHx_Old_Px_I=repmat(CBxHx_Old_Px_I,Ly+1,1);
CBxHx_Old_Px_II=fliplr(CBxHx_Old_Px_I);
% CHy_Px
CHy_Px_I=(1-eaf_xH)./(1+eaf_xH);
CHy_Px_I=repmat(CHy_Px_I,Ly,1);
CHy_Px_II=fliplr(CHy_Px_I);
% CEzHy_Px
CEzHy_Px_I=0.5/mu_r./(1+eaf_xH);
CEzHy_Px_I=repmat(CEzHy_Px_I,Ly,1);
CEzHy_Px_II=fliplr(CEzHy_Px_I);
%%%%%%%%%%%%%%%%   Y direction   %%%%%%%%%%%%%%%%
% CEz_Py
CEz_Py_III=(1-eaf_yE)./(1+eaf_yE);
CEz_Py_III=repmat(CEz_Py_III,1,Lx);
CEz_Py_IV=flipud(CEz_Py_III);
% CHxEz_Py
CHxEz_Py_III=0.5/ep_r./(1+eaf_yE);
CHxEz_Py_III=repmat(CHxEz_Py_III,1,Lx);
CHxEz_Py_IV=flipud(CHxEz_Py_III);