%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%平面波展開法
%二維聲子晶體帶結構計算
%計算二維正方格子
%散射體立于基體之中
%%%%%%%%%%%%%%%%%%%%%%%%%
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）=