%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%边界条件采用Berenger分裂场PML, 计算 TM模式
%TE模式,只需要将程序中的EZ->HZ, HX->-EX,HY->-EY,SigmaE->SigmaM,EPS->Miu即可
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
clc;
WantToSeemovie=1; %是否播放动画
WantToSeeSigma=1; %是否绘制eps分布图
Ntime=1000; %总共的时间步
Interval=40; %播放动画时两次绘图间隔的时间步数
a=1.0e-6; %晶格常数,在这里只是为了归一化频率的方便。
%如果是计算晶体散射,根据需要自行设定并设置好晶格电磁参数的分布。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%以下程序段定义计算采用的物理常数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
miu0=4*pi*1.0e-7; %自由空间磁导率.
eps0=8.85*1e-12; %自由空间介电常数.
cc=1/sqrt(miu0*eps0); %自由空间光速
factor=miu0/eps0; %自由空间磁导率与自由空间介电常数比值
Omiga0=1;
Omiga=Omiga0*2*pi*cc/a;
Lumda=2*pi*cc/Omiga;
epsa = 13; %介质柱的介电常数
epsb = 1; %背景的介电常数
h_pl=6.626068*1.0e-34; %普朗克常数
OmigaP=15*1.60217646*1.0e-19*2*pi/h_pl; %铝的OmigaP
SigmaE_Al=3.57143*1.0e+7; %铝的SigmaE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%生成计算主空间离散化矩阵以及介电系数,磁导系数,电导率,磁导率的空间分布矩阵
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Nx=600; %主计算区X方向上的FDTD的格子数,注意不是晶格的格子数。
Ny=600; %主计算区Y方向上的FDTD的格子数
s1=1;
s2=Ny/2; %激励源的位置
DD=a/20; %采用均匀网格
dt=DD/cc/sqrt(2);
%定义 主空间中的Ez, H1, H2矩阵,
ez=zeros(Nx,Ny);
j_ez=ez;
ez_temp=ez;
hx=zeros(Nx,Ny+1);
hy=zeros(Nx+1,Ny);
%以下定义ez的坐标及物理参数 要注意到ez,与hx,hy空间上分别有半各步长之差
X_ez=((1:Nx)-0.5)*DD-Nx*DD/2;
Y_ez=((1:Ny)-0.5)*DD-Ny*DD/2;
[XX_ez,YY_ez]=meshgrid(X_ez, Y_ez);
XX_ez=XX_ez';
YY_ez=YY_ez';
d1=Nx/5*DD;
d2=Ny/4*DD;
SigmaE_ez=zeros(size(ez));
for i=1:length(X_ez)
for j=1:length(Y_ez)
if X_ez(i)>=-d1&&X_ez(i)<=d1&&Y_ez(j)>=-d2&&Y_ez(j)<=d2
SigmaE_ez(i,j)=SigmaE_Al; %铝的SigmaE区域
end
end
end
SigmaE_ez(:,Ny/2-Ny/8:Ny/2+Ny/8)=0;
flag=find(SigmaE_ez~=0);
Eps_ez=ones(size(ez))*epsb;
Eps_ez(flag)=1; %有损金属计算中需要的是eps0;
Eps_ez=Eps_ez*eps0;
%以下定义hx的坐标及物理参数
X_hx=((1:Nx)-0.5)*DD-Nx*DD/2;
Y_hx=((1:Ny+1)-1)*DD-Ny*DD/2;
[XX_hx,YY_hx]=meshgrid(X_hx, Y_hx);
XX_hx=XX_hx';
YY_hx=YY_hx';
SigmaM_hx=zeros(size(hx));
Miu_hx=ones(size(hx))*miu0;
%以下定义hy的坐标及物理参数
X_hy=((1:Nx+1)-1)*DD-Nx*DD/2;
Y_hy=((1:Ny)-0.5)*DD-Ny*DD/2;
[XX_hy,YY_hy]=meshgrid(X_hy, Y_hy);
XX_hy=XX_hy';
YY_hy=YY_hy';
SigmaM_hy=zeros(size(hy));
Miu_hy=ones(size(hy))*miu0;
%生成主计算区域hx hy之CP,CQ,以及EZ之CA,CB
CP_hx=(1-SigmaM_hx.*dt./2./Miu_hx)./(1+SigmaM_hx.*dt./2./Miu_hx);
CQ_hx=(dt./Miu_hx)./(1+SigmaM_hx.*dt./2./Miu_hx);
CP_hy=(1-SigmaM_hy.*dt./2./Miu_hy)./(1+SigmaM_hy.*dt./2./Miu_hy);
CQ_hy=(dt./Miu_hy)./(1+SigmaM_hy.*dt./2./Miu_hy);
Gamma_ez=zeros(size(ez));
Gamma_ez(flag)=OmigaP^2*eps0./SigmaE_ez(flag); %避免出现非金属区域电导率为0时除0错误
OmigaP_ez=zeros(size(ez));
OmigaP_ez(flag)=OmigaP;
K_ez=(1-Gamma_ez.*dt/2)./(1+Gamma_ez.*dt/2);
B_ez=eps0*OmigaP_ez.^2*dt/2./(1+Gamma_ez.*dt/2);
SigmaE_ez_eff=SigmaE_ez+B_ez;
CA=(1-SigmaE_ez_eff.*dt./2./Eps_ez)./(1+SigmaE_ez_eff.*dt./2./Eps_ez);
CB=(dt./Eps_ez)./(1+SigmaE_ez_eff.*dt./2./Eps_ez);
CD=(1+K_ez).*dt./Eps_ez./2./(1+SigmaE_ez_eff.*dt./2./Eps_ez);
%以下程序段为显示Sigma的分布
if WantToSeeSigma==1
% pcolor(XX_ez',YY_ez',real(Eps_ez')./eps0);
pcolor(XX_ez',YY_ez',SigmaE_ez');
shading interp;
view(0,90);
axis([min(X_ez), max(X_ez),min(Y_ez), max(Y_ez)])
axis off;
disp('Press any key to continue...');
pause
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%初始化PML层的场分布,以及填充PML区域
% --------------------------
% 3 | B | 4
% ----|-----------------|---
% | |
% C | MAIN GRID | D
% | |
% ----|-----------------|---
% 1 | A | 2
% --------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%定义PML参数
NPML=12; %PML的层数
delta=NPML*DD;
m=4; %采用逐层增大分布,此为阶数
R0=1e-10; %容许正入射时的反射率
SigmaEMax=-(m+1)*eps0*cc*log(R0)/(delta*2); %X,Y方向上的最大sigma,注意PML中存在SigmaE/eps0=SigmaM*miu0的关系。
% SigmaBound=SigmaEMax*(DD/2).^(m+1)/(delta^m*DD*(m+1))/eps0;
% expboundary=exp(-SigmaBound*dt);
% 另一种算法,即在场更新的部分使用的是有expboundary表达式的,这两行需要uncomment
NUM=1:2*NPML;
SigmaE_x=SigmaEMax*(NUM/2/NPML).^(m+1); %因为ezx(ezy),hx,hy采用的参数均有半各步长之差
SigmaE_y=SigmaE_x; %故定义2倍NPML,视场的需要采用。比如,hx计算需要
SigmaM_x=SigmaE_x*miu0/eps0; %SigmaM_y,而hy的计算则需要SigmaM_x,而 ezx计算
SigmaM_y=SigmaE_y*miu0/eps0; %需要SigmaE_x,而ezy计算需要SigmaE_y,因此要注意
% 区域1
ezxPML1=zeros(NPML,NPML);
ezyPML1=zeros(NPML,NPML);
hxPML1 =zeros(NPML,NPML);
hyPML1 =zeros(NPML,NPML);
SigmaE_x_PML1=repmat(SigmaE_x(NPML*2-1:-2:1)',1,NPML);
SigmaE_y_PML1=repmat(SigmaE_y(NPML*2-1:-2:1),NPML,1);
SigmaM_x_PML1=repmat(SigmaM_x(NPML*2:-2:2)',1,NPML);
SigmaM_y_PML1=repmat(SigmaM_y(NPML*2:-2:2),NPML,1);
% 区域2
ezxPML2=zeros(NPML,NPML);
ezyPML2=zeros(NPML,NPML);
hxPML2 =zeros(NPML,NPML);
hyPML2 =zeros(NPML,NPML);
SigmaE_x_PML2=flipud(SigmaE_x_PML1);
SigmaE_y_PML2=SigmaE_y_PML1;
SigmaM_x_PML2=flipud(SigmaM_x_PML1);
SigmaM_y_PML2=SigmaM_y_PML1;
% 区域3
ezxPML3=zeros(NPML,NPML);
ezyPML3=zeros(NPML,NPML);
hxPML3 =zeros(NPML,NPML);
hyPML3 =zeros(NPML,NPML);
SigmaE_x_PML3=SigmaE_x_PML1;
SigmaE_y_PML3=fliplr(SigmaE_y_PML1);
SigmaM_x_PML3=SigmaM_x_PML1;
SigmaM_y_PML3=fliplr(SigmaM_y_PML1);
% 区域4
ezxPML4=zeros(NPML,NPML);
ezyPML4=zeros(NPML,NPML);
hxPML4 =zeros(NPML,NPML);
hyPML4 =zeros(NPML,NPML);
SigmaE_x_PML4=SigmaE_x_PML2;
SigmaE_y_PML4=SigmaE_y_PML3;
SigmaM_x_PML4=SigmaM_x_PML2;
SigmaM_y_PML4=SigmaM_y_PML3;
% 区域A
ezxPMLA=zeros(Nx,NPML);
ezyPMLA=zeros(Nx,NPML);
hxPMLA =zeros(Nx,NPML);
hyPMLA =zeros(Nx+1,NPML);
SigmaE_x_PMLA=zeros(Nx,NPML);
SigmaE_y_PMLA=repmat(SigmaE_y(NPML*2-1:-2:1),Nx,1);
SigmaM_x_PMLA=zeros(Nx+1,NPML);
SigmaM_y_PMLA=repmat(SigmaM_y(NPML*2:-2:2),Nx,1);
% 区域B
ezxPMLB=zeros(Nx,NPML);
ezyPMLB=zeros(Nx,NPML);
hxPMLB =zeros(Nx,NPML);
hyPMLB =zeros(Nx+1,NPML);
SigmaE_x_PMLB=zeros(Nx,NPML);
SigmaE_y_PMLB=fliplr(SigmaE_y_PMLA);
SigmaM_x_PMLB=zeros(Nx+1,NPML);
SigmaM_y_PMLB=fliplr(SigmaM_y_PMLA);
% 区域C
ezxPMLC=zeros(NPML,Ny);
ezyPMLC=zeros(NPML,Ny);
hxPMLC =zeros(NPML,Ny+1);
hyPMLC =zeros(NPML,Ny);
SigmaE_x_PMLC=repmat(SigmaE_y(NPML*2-1:-2:1)',1, Ny);
SigmaE_y_PMLC=zeros(NPML,Ny);
SigmaM_x_PMLC=repmat(SigmaM_y(NPML*2:-2:2)',1, Ny);
SigmaM_y_PMLC=zeros(NPML,Ny+1);
% 区域D
ezxPMLD=zeros(NPML,Ny);
ezyPMLD=zeros(NPML,Ny);
hxPMLD =zeros(NPML,Ny+1);
hyPMLD =zeros(NPML,Ny);
SigmaE_x_PMLD=flipud(SigmaE_x_PMLC);
SigmaE_y_PMLD=zeros(NPML,Ny);
SigmaM_x_PMLD=flipud(SigmaM_x_PMLC);
SigmaM_y_PMLD=zeros(NPML,Ny+1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%可视化的初始化
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ntx=Nx+2*NPML;
Nty=Ny+2*NPML;
x=((1:Ntx)-0.5)*DD-Ntx*DD/2;
y=((1:Nty)-0.5)*DD-Nty*DD/2;
[X,Y]