平面波展开法计算二维准零刚度声子晶体带隙的matlab代码

以下是平面波展开法计算二维准零刚度声子晶体带隙的Matlab代码：

clear all;

%定义晶格常数a和波矢范围
a = 1; %晶格常数
kxrange = -pi/a:0.01:pi/a;
kyrange = -pi/a:0.01:pi/a;
Nk = length(kxrange)*length(kyrange); %波矢点数

%定义材料参数
r = 0.3*a; %圆柱半径
e1 = 1; %基底介电常数
e2 = 11.56; %圆柱介电常数
filling_ratio = pi*r^2/a^2; %填充率

%计算矩阵元
M = zeros(Nk,Nk);
for i = 1:Nk
    for j = 1:Nk
        kx1 = kxrange(mod(i-1,length(kxrange))+1);
        ky1 = kyrange(ceil(i/length(kxrange)));
        kx2 = kxrange(mod(j-1,length(kxrange))+1);
        ky2 = kyrange(ceil(j/length(kxrange)));
        if i==j
            M(i,j) = dot_product(kx1,ky1,kx2,ky2,e1,e2,filling_ratio,a);
        else
            M(i,j) = dot_product(kx1,ky1,kx2,ky2,e1,e2,filling_ratio,a)*exp(-1i*dot_product(kx1,ky1,kx2,ky2,0,0,filling_ratio,a));
        end
    end
end

%求解本征值和本征矢
[V,D] = eig(M);
omega2 = diag(D);
[omega2,idx] = sort(omega2);
V = V(:,idx);

%绘制能带图
omega = sqrt(omega2);
Nbands = 10; %能带数目
for nb = 1:Nbands
    E = omega((nb-1)*Nk+1:nb*Nk);
    E = reshape(E,length(kxrange),length(kyrange));
    contour(kxrange,kyrange,E,[E(1,1) E(1,1)],'k');
    hold on;
end
axis equal;
xlabel('k_x');
ylabel('k_y');

%计算两个波矢的点积
function dotp = dot_product(kx1,ky1,kx2,ky2,e1,e2,filling_ratio,a)
    G = 2*pi/a;
    dotp = (e1+e2)/2*(kx1-kx2)^2+(e1+e2)/2*(ky1-ky2)^2;
    if abs(kx1-kx2+G)<1e-10 && abs(ky1-ky2)<1e-10
        dotp = dotp - e2*filling_ratio;
    end
    if abs(kx1-kx2)<1e-10 && abs(ky1-ky2+G)<1e-10
        dotp = dotp - e2*filling_ratio;
    end
    if abs(kx1-kx2-G)<1e-10 && abs(ky1-ky2)<1e-10
        dotp = dotp - e2*filling_ratio;
    end
    if abs(kx1-kx2)<1e-10 && abs(ky1-ky2-G)<1e-10
        dotp = dotp - e2*filling_ratio;
    end
end

其中，`dot_product`函数计算了两个波矢的点积，`M`矩阵元素由该函数计算得到，然后通过对`M`求本征值和本征矢，得到声子晶体的能带图。