平面波展开法计算二维准零刚度声子晶体带隙特性的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);

%计算频率范围
freq_min = sqrt(omega2(1));
freq_max = sqrt(omega2(end));

%计算频率和带宽
freqs = linspace(freq_min,freq_max,1000); %频率范围
bandwidths = zeros(size(freqs)); %带宽
for i = 1:length(freqs)
    freq = freqs(i);
    omega = freq^2; %角频率
    idx = find(omega2>=omega,1); %找到第一个大于等于该角频率的本征值下标
    if isempty(idx)
        bandwidths(i) = 0; %没有带隙
    else
        E1 = sqrt(omega2(idx-1)); %能带底
        E2 = sqrt(omega2(idx)); %能带顶
        bandwidths(i) = E2 - E1;
    end
end

%绘制频率和带宽曲线
plot(freqs,bandwidths);
xlabel('频率');
ylabel('带宽');

其中，计算频率范围的代码使用了声子晶体的本征频率的最小值和最大值，然后通过在这个范围内等间距取1000个点，计算每个频率对应的带宽。最后，绘制频率和带宽曲线。