平面波展开高静低动局域共振声子晶体带隙计算MATLAB程序

平面波展开方法是计算晶体带隙的一种常用方法，而高静低动局域共振声子晶体则是一种特殊的声子晶体结构。下面给出MATLAB程序，用于计算该结构的带隙。

程序中涉及到的变量和参数含义如下：

- a：晶格常数
- r：球体半径
- k：波矢量范围
- n：平面波展开阶数
- N：晶胞数
- M：k点数
- omega：角频率
- eig_val：特征值
- eig_vec：特征向量
- E：能量
- gap：带隙宽度

% 确定参数
a = 1; % 晶格常数
r = 0.3*a; % 球体半径
k = 2*pi/a; % 波矢量范围
n = 4; % 平面波展开阶数
N = 10; % 晶胞数
M = 100; % k点数

% 计算k点
kx = linspace(-k,k,M);
ky = linspace(-k,k,M);
kz = linspace(-k,k,M);

% 构造晶体
[x,y,z] = meshgrid(0:a:a*(N-1),0:a:a*(N-1),0:a:a*(N-1));
X = [x(:),y(:),z(:)];

% 计算矩阵
Mk = zeros(N^3);
for i = 1:N^3
    for j = 1:N^3
        if i~=j
            Rij = X(i,:)-X(j,:);
            rij = norm(Rij);
            if rij<=a
                Mk(i,j) = -(r/rij)^n*exp(-1i*k*Rij*',[1,1,1]'); % Rij*[1,1,1]'为计算k
            end
        end
    end
end

% 计算特征值和特征向量
[eig_vec,eig_val] = eig(Mk);

% 排序特征值和特征向量
eig_val = diag(eig_val);
[eig_val,idx] = sort(eig_val);
eig_vec = eig_vec(:,idx);

% 计算能带
omega = sqrt(abs(eig_val));
E = omega/(2*pi);
E = reshape(E,N^3,M);
kx = ones(N^3,1)*kx;
ky = ones(N^3,1)*ky;
kz = ones(N^3,1)*kz;
kx = kx(:);
ky = ky(:);
kz = kz(:);
k = sqrt(kx.^2+ky.^2+kz.^2);
[k,idx] = sort(k);
E = E(:,idx);

% 计算带隙宽度
gap = min(E(:,M/2+1))-max(E(:,M/2));
disp(['The gap width is ',num2str(gap),' Hz']);

注意，在实际计算中，需要根据具体的晶体结构和参数进行调整。