二维声子晶体平面波展开法计算能带

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%平面波展开法
%二维声子晶体带结构计算
%计算二维正方格子
%散射体立于基体之中
%%%%%%%%%%%%%%%%%%%%%%%%%
clear;clc;tic;epssys=1.0e-6; %设定一个最小量，避免系统截断误差或除零错误
 
%%%%%%%%%%%%%%%%%%%%%%%%%%
%定义实际的正空间格子基矢
%%%%%%%%%%%%%%%%%%%%%%%%%%
a=0.02;
a1=a*[1 0];
a2=a*[0 1];
%%%%%%%%%%%%%%%%%%%%%%%%%%
%定义晶格的参数
%%%%%%%%%%%%%%%%%%%%%%%%%%
rho1=11600;E1=4.08e10;mju1=1.49e10;lambda1=mju1*（E1-2*mju1）/（3*mju1-E1）; %散射体的材料参数
rho2=1300;E2=1.175e5;mju2=4e4;lambda2=mju2*（E2-2*mju2）/（3*mju2-E2）; %基体的材料参数
Rc=0.006; %散射体截面半径
Ac=pi*（Rc）^2; %散射体截面面积
Au=a^2; %二维格子原胞面积
Pf=Ac/Au; %填充率
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%生成倒格基矢
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b1=2*pi/a*[1 0];
b2=2*pi/a*[0 1];
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%选定参与运算的倒空间格矢量，即参与运算的平面波数量
%设定一个lm的取值范围，变化lm即可得出参与运算的平面波集合
NrSquare=10; %选定倒空间的尺度，即lm（倒格矢G=l*b1+m*b2）的取值范围。
             %NrSquare确定后，使用Bloch波数目可能为（2*NrSquare+1）^2
G=zeros（（2*NrSquare+1）^22）; %初始化可能使用的倒格矢矩阵
i=1;
for l=-NrSquare:NrSquare
    for m=-NrSquare:NrSquare
        G（i:）=l*b1+m*b2;
        i=i+1;
    end;
end;
NG=i-1; %实际使用的Bloch波数目
G=G（1:NG:）; 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%生成k空间的rho（Gi-Gj）mju（Gi-Gj）lambda（Gi-Gj）值，ij从1到NG。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rho=zeros（NGNG）;mju=zeros（NGNG）;lambda=zeros（NGNG）;
for i=1:NG
    for j=1:NG
        Gij=norm（G（j:）-G（i:））;
        if （Gij            rho（ij）=rho1*Pf+rho2*（1-Pf）;
            mju（ij）=mju1*Pf+mju2*（1-Pf）;
            lambda（ij）=lambda1*Pf+lambda2*（1-Pf）;
        else
            rho（ij）=（rho1-rho2）*2*Pf*besselj（1Gij*Rc）/（Gij*Rc）;
            mju（ij）=（mju1-mju2）*2*Pf*besselj（1Gij*Rc）/（Gij*Rc）;
            lambda（ij）=（lambda1-lambda2）*2*Pf*besselj（1Gij*Rc）/（Gij*Rc）;
        end;
    end;
end;
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%定义简约布里渊区的各高对称点
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T=（2*pi/a）*[epssys 0];
M=（2*pi/a）*[1/2 1/2];
X=（2*pi/a）*[1/2 0];
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%对于简约布里渊区边界上的每个k，求解其特征频率
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
THETA_A=zeros（NGNG）; %待解的本征方程A矩阵
THETA_B=zeros（NGNG）; %待解的本征方程B矩阵
Nkpoints=10; %每个方向上取的点数
stepsize=0:1/（Nkpoints-1）:1; %每个方向上步长
TX_eig=zeros（NkpointsNG）; %沿TX方向的波的待解的特征频率矩阵
XM_eig=zeros（NkpointsNG）; %沿XM方向的波的待解的特征频率矩阵
MT_eig=zeros（NkpointsNG）; %沿MT方向的波的待解的特征频率矩阵
for n=1:Nkpoints
    fprintf（[‘\n k-point:‘int2str（n）‘of‘int2str（Nkpoints）‘.\n‘]）;
     
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %对于TX（正方格子）方向上的每个k值，求解其特征频率
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    TX_step=stepsize（n）*（X-T）+T;
     
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %n 求本征矩阵的元素
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1:NG
        for j=1:NG
            kGi=TX_step+G（i:）;
            kGj=TX_step+G（j:）;
            THETA_A（ij）=mju（ij）*dot（kGikGj）;
            THETA_B（ij）=