平面波展开法计算二维弹簧质量声子晶体超材料板带隙matlab程序

平面波展开法（Plane Wave Expansion Method，PWEM）可以用于计算二维弹簧质量声子晶体超材料板的带隙。Matlab程序如下：

%% 声子晶体超材料板平面波展开法计算程序

%% 参数设置
a = 1; % 基本元胞尺寸
h = a/2; % 板厚
r = 0.2*a; % 弹簧半径
kmax = 10; % 最大平面波波矢数
n = 20; % 基本元胞个数
eps1 = 1; % 材料1介电常数
eps2 = 4; % 材料2介电常数

%% 构造二维弹簧质量声子晶体超材料板
% 构造二维弹簧质量声子晶体超材料板的基本元胞
[x,y] = meshgrid(0:a/(n-1):a,0:a/(n-1):a); % 坐标网格
x = x'; y = y';
lattice1 = ones(size(x)); % 材料1
lattice2 = ones(size(x)); % 材料2
for i = 1:n
    for j = 1:n
        if mod(i+j,2) == 0
            lattice1(i,j) = 0; % 材料1
        else
            lattice2(i,j) = 0; % 材料2
        end
    end
end
% 构造弹簧
spring = (x-a/2).^2 + (y-a/2).^2 < r^2;
% 构造弹簧质量声子晶体超材料板
structure = lattice1.*(~spring)*eps1 + lattice2.*(spring)*eps2;

%% 计算平面波展开系数
[Kx,Ky] = meshgrid(-kmax:kmax,-kmax:kmax);
Kx = Kx(:); Ky = Ky(:);
K = [Kx Ky]*2*pi/a;
num_k = length(K)
% 计算平面波展开系数
C = zeros(num_k,num_k);
for i = 1:num_k
    for j = 1:num_k
        C(i,j) = sum(sum(exp(1i*sum([K(i,:);-K(j,:)]*...
            [x(:)';y(:)'])')));
    end
end

%% 计算频率-波矢关系曲线
w0 = 2*pi*1e10; % 基频角频率
w = linspace(0.1*w0,2*w0,100); % 计算频率范围
num_w = length(w);
kx = zeros(num_k,num_w);
ky = zeros(num_k,num_w);
for i = 1:num_w
    H = zeros(size(x));
    for j = 1:num_k
        H = H + C(:,j)*exp(1i*K(j,:)*[x(:)';y(:)'])';
    end
    [Ex,Ey] = gradient(H);
    kx(:,i) = w(i)/3e8*Ex(:);
    ky(:,i) = w(i)/3e8*Ey(:);
end

%% 绘制频率-波矢关系曲线
figure;
plot(kx,ky,'.','MarkerSize',10);
xlabel('k_x (1/m)');
ylabel('k_y (1/m)');
title('Phononic Crystal Slab');

%% 计算波导模式
% 以k=(0,0)为中心切出一个正方形，作为波导
kx_wg = (abs(kx)<0.2*max(abs(kx(:)))) & (abs(ky)<0.2*max(abs(ky(:))));
kx_wg_center = round(num_k/2);
kx_wg(:,:,kx_wg_center) = 1;
% 计算波导模式
H = zeros(size(x));
for i = 1:num_k
    for j = 1:num_w
        if kx_wg(i,j,kx_wg_center)
            H = H + C(:,i)*exp(1i*K(i,:)*[x(:)';y(:)'])*...
                exp(-1i*w(j)*h);
        end
    end
end

%% 绘制波导模式
figure;
imagesc(abs(H));
axis equal off;
title('Waveguide Mode');

程序中使用了二维弹簧质量声子晶体超材料板的基本元胞构造了结构，然后计算了平面波展开系数，进而计算了频率-波矢关系曲线。最后，以k=(0,0)为中心切出一个正方形，作为波导，并计算了波导模式。运行程序后，会分别绘制频率-波矢关系曲线和波导模式的图像。