平面波展开弹簧质量声子晶体带隙计算MATLAB程序

以下是平面波展开弹簧质量声子晶体带隙计算的MATLAB程序示例：

% 声子晶体参数
a = 1; % 晶格常数
r = 0.3; % 球体半径
c = 2; % 声速
M = 2; % 弹簧质量
n = 10; % 布拉格平面波展开的阶数
kx = linspace(-pi/a, pi/a, 100); % kx范围
ky = linspace(-pi/a, pi/a, 100); % ky范围
[KKX, KKY] = meshgrid(kx, ky);

% 计算频率
w = zeros(length(kx), length(ky));
for i = 1:length(kx)
    for j = 1:length(ky)
        [E,~,~,~] = dispersion(KKX(i,j),KKY(i,j),a,r,c,M,n);
        w(i,j) = E;
    end
end

% 绘制带隙图
figure
contourf(kx, ky, w, 100, 'LineColor', 'none')
xlabel('k_x')
ylabel('k_y')
colorbar
title('Band structure of phononic crystal')

function [E, dE_dkx, dE_dky, dE_dkz] = dispersion(kx, ky, a, r, c, M, n)
% 计算晶体中的频率和导数
% 输入：kx、ky、a、r、c、M、n
% 输出：E、dE_dkx、dE_dky、dE_dkz
% E为频率，dE_dkx、dE_dky、dE_dkz为对应的导数

% 计算布拉格平面波展开系数
G_vec = zeros(2*n+1, 2*n+1, 2*n+1, 3);
G_abs = zeros(2*n+1, 2*n+1, 2*n+1);
for i = -n:n
    for j = -n:n
        for k = -n:n
            if i==0 && j==0 && k==0
                continue
            end
            G_vec(i+n+1,j+n+1,k+n+1,:) = 2*pi/a*[i,j,k];
            G_abs(i+n+1,j+n+1,k+n+1) = norm(G_vec(i+n+1,j+n+1,k+n+1,:));
        end
    end
end
G_vec = reshape(G_vec, (2*n+1)^3, 3)';
G_abs = reshape(G_abs, (2*n+1)^3, 1)';
[~,I] = sort(G_abs);
G_vec = G_vec(:,I);

% 计算弹簧常数矩阵
N = (2*n+1)^3;
K = zeros(3*N, 3*N);
for i = 1:N
    for j = 1:N
        if i==j
            K(3*(i-1)+1:3*i, 3*(i-1)+1:3*i) = diag([M,M,M]);
            continue
        end
        r_ij = G_vec(:,j) - G_vec(:,i);
        r_ij_p = r_ij + a*[1,0,0]';
        r_ij_m = r_ij - a*[1,0,0]';
        r_ij_q = r_ij + a*[0,1,0]';
        r_ij_n = r_ij - a*[0,1,0]';
        r_ij_r = r_ij + a*[0,0,1]';
        r_ij_l = r_ij - a*[0,0,1]';
        if norm(r_ij_p)<=r
            K(3*(i-1)+1:3*i, 3*(j-1)+1:3*j) = [1,-1,0;-1,1,0;0,0,0];
        end
        if norm(r_ij_m)<=r
            K(3*(i-1)+1:3*i, 3*(j-1)+1:3*j) = [1,-1,0;-1,1,0;0,0,0];
        end
        if norm(r_ij_q)<=r
            K(3*(i-1)+1:3*i, 3*(j-1)+1:3*j) = [1,0,-1;0,0,0;-1,0,1];
        end
        if norm(r_ij_n)<=r
            K(3*(i-1)+1:3*i, 3*(j-1)+1:3*j) = [1,0,-1;0,0,0;-1,0,1];
        end
        if norm(r_ij_r)<=r
            K(3*(i-1)+1:3*i, 3*(j-1)+1:3*j) = [1,0,0;0,1,0;0,0,0];
        end
        if norm(r_ij_l)<=r
            K(3*(i-1)+1:3*i, 3*(j-1)+1:3*j) = [1,0,0;0,1,0;0,0,0];
        end
    end
end

% 计算本征值和本征向量
KK = kx*[1,0,0]' + ky*[0,1,0]';
KKK = repmat(KK, [N,1]) + repmat(G_vec', [1,length(kx)]);
KKK2 = sum(KKK.^2, 2);
KKK3 = sum(KKK.^3, 2);
KKK4 = sum(KKK.^4, 2);
A = K - c^2*diag(KKK2.^(-1))*diag(KKK4)*diag(KKK2.^(-1)) + diag(KKK2.^(-1))*diag(KKK3)*diag(KKK2.^(-1));
[V,D] = eig(A);
[~,I] = sort(diag(D));

% 返回频率和导数
E = sqrt(diag(D(I)));
dE_dkx = zeros(length(E), 1);
dE_dky = zeros(length(E), 1);
dE_dkz = zeros(length(E), 1);
for i = 1:length(E)
    psi = V(:,I(i));
    G_dpsi = zeros(3*N,3);
    for j = 1:3
        G_dpsi(:,j) = G_vec(:,j).*psi;
    end
    G_dpsi_dot_Kpsi = sum(G_dpsi.*K*psi, 1)';
    dE_dkx(i) = real(G_dpsi_dot_Kpsi(1)/2/E(i));
    dE_dky(i) = real(G_dpsi_dot_Kpsi(2)/2/E(i));
    dE_dkz(i) = real(G_dpsi_dot_Kpsi(3)/2/E(i));
end
end

该程序可以计算平面波展开弹簧质量声子晶体的频率和导数，并绘制带隙图。输入参数包括晶格常数、球体半径、声速、弹簧质量、布拉格平面波展开的阶数和kx、ky范围。输出包括频率和对应的导数。