9307785PWE_pononic.rar_Phononic Crystals_bluecoh_二维声子晶体_声子晶体_声子晶

%正方格子圆形夹杂二维声子晶体XY模态禁带图 clc; clear; %硬散射体材料参数-铅 rou1=2730; lbd1=6.89e10; mu1=2.87e10; %基体材料参数-环氧树脂 rou2=1180; lbd2=4.43e9; mu2=1.59e9; %晶格常数 a=0.02; R=0.008; %硬夹杂散射体的填充比 f=pi*R^2/a^2; %计算、禁带图设置 m=10; %取的晶格周期数：2*m+1 n=10; %每段波失的取值数：n b=8; %画出的频率曲线数 %倒格矢取值 Gx=2*pi/a*[-m:m]; Gy=2*pi/a*[-m:m]; !G=[Gx;Gy]; G=outG(Gx,Gy); %倒格矢组合，组合个数即是计算使用的平面波数（2*m+1)^2 !G3x=Gx*2*[-2*m:2*m]; !G3y=Gy*2*[-2*m:2*m]; !G3=[G3x;G3y]; !G=(Gx.^2+Gy.^2)^(1/2); %形成计算禁带用的倒格矢参数矩阵 rou=GG(rou1,rou2,f,R,G); lbd=GG(lbd1,lbd2,f,R,G); mu=GG(mu1,mu2,f,R,G); !nrou=length(rou); !for i=1:nrou; !N((i-1)*2+1,:)=[rou(i,:) zeros(1,nrou)]; !N((i-1)*2+2,:)=[zeros(1,nrou) rou(i,:)]; !end %形成矩阵N N=[rou zeros((2*m+1)^2,(2*m+1)^2);zeros((2*m+1)^2,(2*m+1)^2) rou]; %生成M点--Gama点方向变化的波失 kx=linspace(pi/a,0,n); ky=linspace(pi/a,0,n); k=[kx;ky]; kmo1=(kx.^2+ky.^2).^(1/2); %计算波失沿M点--Gama点方向变化时的特征值 omega1=omega(k,G,lbd,mu,N,b); %生成Gama点--X点方向变化的波失 kx=linspace(0,pi/a,n); ky=linspace(0,0,n); k=[kx;ky]; kmo2=(kx.^2+ky.^2).^(1/2); %计算波失沿Gama点--X点方向变化时的特征值 omega2=omega(k,G,lbd,mu,N,b); %生成X点--M点方向变化的波失 kx=linspace(pi/a,pi/a,n); ky=linspace(0,pi/a,n); k=[kx;ky]; kmo3=(kx.^2+ky.^2).^(1/2); %计算波失沿X点--M点方向变化时的特征值 omega3=omega(k,G,lbd,mu,N,b); %连接M-Gama-X-M方向变化的波失 kmo=[kmo1 kmo2(2:n) kmo3(2:n)]; %连接计算出的M-Gama-X-M方向波失变化对应的特征值 omega=[omega1;omega2(2:n,:);omega3(2:n,:)]; !plot(kmo,real(omega)); %特征值肯定是实数！ %波的频率 f=omega/2/pi; %画出禁带图！ plot(0:length(kmo)-1,f);