states_densityofphc.m

计算光子晶体态密度的程序

clc
clear
tic
a=2.32*10;

N=7;%(N的选择不应该影响到最终能带图）

ea=14;%ea=1; 

eb=1;%eb=12.5;
f=0.431;
tic;

disp('--------------------------------------------------')

r=sqrt(f*sqrt(3)*a*a/(2*pi));
stable=(4*pi/a)/3;
b1=2*pi/a-i*2*pi/(sqrt(3)*a);b2=i*4*pi/(sqrt(3)*a);
h=-N:N;
G1=repmat(h,(2*N+1),1);
G2=G1';
G=G1*b1+G2*b2
G3=G(:);
G4=repmat(G3,1,(2*N+1).^2);
G5=conj(G4');

Gr=r*abs(G4-G5);
na=find(Gr==0);
   Gr(na)=1;   
    kf=(1/ea-1/eb)*f*2*Besselj(1,Gr)./Gr;
    kf(na)=f/ea+(1-f)/eb;

  
   k1=linspace(0,stable*sqrt(3)/2,360); 
   k11=-sqrt(3)*k1/2-i*k1/2;
   k2=linspace(0,stable,400);
   k22=k2;
   k3=linspace(0,stable/2,200);
   k33=(stable-k3/2)+i*sqrt(3)*k3/2;
   k4=[k11,k22,k33];
   
   s=[];
   for m=1:length(k4);
       k=k4(m);
    Q=(real(k+G4).*real(k+G5)+imag(k+G4).*imag(k+G5)).*kf;%%%%%%%%%
    C=eig(Q);
 
        for l=1:(2*N+1).^2;
          if imag(C(l))==0&real(C(l))>0;
            E(l)=a*sqrt(C(l))/(2*pi);
          else
           E(l)=0 ;
          end
        end
       s=[s;E];
end
w=linspace(0,0.9,100);
w2=[];
num2=[];
for m=1:length(w);
    w1=w(m);
   a=find(s>w1-0.45/100&s<w1+0.45/100);
   num1=length(a);
   w2=[w2,w1];
   num2=[num2,num1];
end
plot(num2,w2)
toc