用平面波展开法计算准零刚度超材料板带隙MATLAB程序

准零刚度超材料板的平面波展开法计算程序如下：

%% 参数设置
a = 1;              % 基元尺寸
r = 0.3*a;          % 孔半径
h = 0.5*a;          % 基板厚度
e_r = 2;            % 介电常数
mu_r = 1;           % 磁导率
d = 0.4*a;          % 波导宽度
lambda = 2*a;       % 波长
k = 2*pi/lambda;    % 波数
N = 15;             % 平面波展开阶数

%% 构造Brillouin区
G1 = 2*pi/a*[0,1,1];
G2 = 2*pi/a*[1,0,1];
G3 = 2*pi/a*[1,1,0];
K1 = 2*pi/a*[1,-1,0];
K2 = 2*pi/a*[1,0,-1];
K3 = 2*pi/a*[0,1,-1];
M = 2*pi/a*[0.5,0.5,0];
BZ = {G1,G2,G3,K1,K2,K3,M};

%% 定义计算函数
function omega2 = calc_omega2(kx,ky,kz,a,r,h,e_r,mu_r,d)
    % 计算平面波展开系数
    G = zeros(2*N+1,2*N+1,2*N+1,3);
    for l = -N:N
        for m = -N:N
            for n = -N:N
                G(l+N+1,m+N+1,n+N+1,:) = k + l*G1 + m*G2 + n*G3;
            end
        end
    end

    % 计算传播常数和反射系数
    omega2 = zeros(size(kx));
    for i = 1:numel(kx)
        K = [kx(i),ky(i),kz];
        S = zeros(6,6);
        for j = -N:N
            for k = -N:N
                for l = -N:N
                    if j==0 && k==0 && l==0
                        continue
                    end
                    Gjkl = G(j+N+1,k+N+1,l+N+1,:);
                    Gjkl2 = dot(Gjkl,Gjkl);
                    if Gjkl2 > (2*pi/lambda)^2
                        continue
                    end
                    R = [1,0,0,0,0,0;0,1,0,0,0,0;0,0,1,0,0,0;
                         0,0,0,1,0,0;0,0,0,0,1,0;0,0,0,0,0,1];
                    Kjkl = K + Gjkl;
                    Kjkl2 = dot(Kjkl,Kjkl);
                    if Kjkl2 > (2*pi/lambda)^2
                        R = zeros(6,6);
                    else
                        p = sqrt((2*pi/lambda)^2-Kjkl2);
                        A = 1/mu_r*[Kjkl(2)*Gjkl(3)-Kjkl(3)*Gjkl(2),...
                                    Kjkl(3)*Gjkl(1)-Kjkl(1)*Gjkl(3),...
                                    Kjkl(1)*Gjkl(2)-Kjkl(2)*Gjkl(1)];
                        B = 1/e_r*[Gjkl(2)*Kjkl(3)-Gjkl(3)*Kjkl(2),...
                                   Gjkl(3)*Kjkl(1)-Gjkl(1)*Kjkl(3),...
                                   Gjkl(1)*Kjkl(2)-Gjkl(2)*Kjkl(1)];
                        R = (1/(2*p))*[A,B;...
                                       Gjkl(1)*p/Kjkl(1),-Kjkl(2)/Kjkl(1)*A;...
                                       Gjkl(1)*p/Gjkl(1),-Gjkl(2)/Gjkl(1)*B];
                    end
                    S = S + exp(-1i*dot(Kjkl,d))*R;
                end
            end
        end
        omega2(i) = det(S);
    end
end

%% 计算频率
nk = 21;
kx = linspace(-pi/a,pi/a,nk);
ky = linspace(-pi/a,pi/a,nk);
[KKX,KKY] = meshgrid(kx,ky);
kz = sqrt((k^2-e_r*kx.^2-e_r*ky.^2)/e_r);
omega2 = calc_omega2(KKX(:),KKY(:),kz(:),a,r,h,e_r,mu_r,d);
omega2 = reshape(omega2,nk,nk);
figure,contour(KKX,KKY,omega2,[0,0],'LineWidth',2);
hold on,plot([-pi/a,pi/a],[0,0],'k--',[0,0],[-pi/a,pi/a],'k--','LineWidth',1);
hold off,axis equal,axis([-pi/a,pi/a,-pi/a,pi/a]*1.1);
title('准零刚度超材料板带隙');
xlabel('k_xa/\pi'),ylabel('k_ya/\pi');

其中，输入参数为基元尺寸$a$、孔半径$r$、基板厚度$h$、介电常数$e_r$、磁导率$\mu_r$、波导宽度$d$、波长$\lambda$、平面波展开阶数$N$。输出为频率的平方$\omega^2$，通过在$k_x-k_y$平面上绘制等频率线来得到带隙。