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

平面波展开法也可以用来计算二维准零刚度声子晶体的带隙特性。这里提供一个简单的Matlab程序供参考：

% 定义晶格
a = 1; % 晶格常数
r = [0 0; 0.5 0.5]; % 原子位置
crystal = Crystal({a*eye(2), r}, 'Pnma');

% 计算弹性常数矩阵
c = 1; % 声速
C = ElasticMatrix(crystal, 'c', c);

% 计算带隙
k = [0 0; 0.5 0.5; 0 0.5; 0 0]; % 波矢路径
nk = 100; % 波矢路径离散化点数
[B, kx, ky] = BandStructure(crystal, C, k, nk, [-1 1]);
plot(B, 'k-', 'LineWidth', 2);
xlabel('Wave vector k');
ylabel('Phonon frequency \omega');

该程序先定义了一个晶格，然后直接计算了弹性常数矩阵，并使用`BandStructure`函数计算并绘制了能带结构图。需要注意的是，由于准零刚度的假设，可能会出现一些误差和限制，需要根据具体情况进行评估和调整。

相关问题

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

在Matlab中，可以使用工具箱中的“Plane Wave Expansion (PWE)”函数来计算二维准零刚度声子晶体的带隙特性。具体步骤如下：

1. 定义晶格：确定二维晶格基向量和布拉吉格子矢量。

2. 定义材料：使用“Crystal”函数定义材料的密度分布。

例如，假设二维准零刚度声子晶体的晶格常数为a，定义材料的密度分布为正弦波形式：

a = 1; % 晶格常数
r = [0 0; 0.5 0.5]; % 原子位置
crystal = Crystal({a*eye(2), r}, 'Pnma', 'Density', @(x,y) sin(2*pi*x/a).*sin(2*pi*y/a));

其中，`Crystal`函数用于定义晶体结构，`Density`参数用于指定材料的密度分布，可以使用函数句柄来自定义密度函数。

3. 计算弹性常数矩阵：使用“ElasticMatrix”函数计算波矢空间中的弹性常数矩阵。

例如，假设声速为c，可以使用如下代码计算弹性常数矩阵：

c = 1; % 声速
C = ElasticMatrix(crystal, 'c', c);

其中，`ElasticMatrix`函数用于计算弹性常数矩阵，`'c'`参数用于指定声速大小。

4. 计算带隙：使用“BandStructure”函数计算声子晶体的带隙特性，并绘制能带结构图。

例如，假设带宽范围为[-0.5, 0.5]，可以使用如下代码计算并绘制能带结构图：

k = [0 0; 0.5 0; 0.5 0.5; 0 0.5]; % 波矢路径
nk = 100; % 波矢路径离散化点数
[B, kx, ky] = BandStructure(crystal, C, k, nk, [-0.5 0.5]);
plot(B, 'k-', 'LineWidth', 2);
xlabel('Wave vector k');
ylabel('Phonon frequency \omega');

其中，`BandStructure`函数用于计算声子晶体的能带结构，`k`参数用于指定波矢路径，`nk`参数用于指定波矢路径的离散化点数，最后一行代码用于绘制能带结构图。

需要注意的是，具体的计算和绘图参数可以根据需要进行调整和优化，例如使用不同的材料、晶格和波矢路径等。

平面波展开法计算二维准零刚度声子晶体带隙的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`求本征值和本征矢，得到声子晶体的能带图。