clc
clear
tic

%a=2.32*10;%參數另一組：ea=1;eb=12.5;f=0.8;縱坐標[0,1.1]
a=1;
l=a/1.63;
bet=0.35;
d=bet*l;
f=l*l/a/a*(1-pi*bet*bet/4);
gamma=pi/4;
ss=0.11;
s=ss*exp(i*gamma);
ea=13.6;%ea=1; 

eb=1;%eb=12.5;


M=[pi/a,pi/a];
U=[0,pi/a];
T=[0,0];
X=[pi/a,0];
X1=[-pi/a,0];
M1=[-pi/a,-pi/a];


N=7;%(N的選擇不應該影響到最終能帶圖）
b1=2*pi/a;
b2=2*i*pi/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');



Gr21=d/2*abs(G4-G5);%貝塞爾函數的自變量
na=find(Gr21==0);
   Gr21(na)=1; 
   Gr22=Besselj(1,Gr21)./Gr21;
   SG23=l*l/a/a*pi*bet*bet/2*Gr22;
   Gr1=(G4-G5);
   SG24=exp(-i*real(Gr1)*real(s)-i*imag(Gr1)*imag(s));
   SG2=SG24.*SG23;
   SG1=l*l/a/a*sinc(l/2/pi*real(Gr1)).*sinc(l/2/pi*imag(Gr1));
  
   kf=(ea-eb)*(SG1-SG2);
   kf(na)=f*ea+(1-f)*eb;
    kf1=inv(kf);
    k55=linspace(0,pi/a,12);
    k5=k55+i*k55;
   
   
      s1=[];
   for m=1:length(k5);
       k=k5(m);
  Q=(real(k+G4).*real(k+G5)+imag(k+G4).*imag(k+G5)).*kf1;%%%%%%%%%
     %Q=abs(k+G4).*kf1.*abs(k+G5);
    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
       s1=[s1;sort(E)];
end
 %k=[stable*sqrt(3)/2-k1,k2+stable*sqrt(3)/2,k3+stable*sqrt(3)/2+stable];
 s1=s1(:,1:9);
 plot(sqrt(2)*k55,s1)