% this program is designed for solving the two dimensional photonic structure% we are using the parameters from one of the paperclearwarning offepsa=1;epsb=13;%epsb=1.46^2; %%%Silica glass air systema=1.0;%f=0.39;%R=sqrt(sqrt(3).*f/(2*pi)*a^2);R=0.45;i=sqrt(-1);f=2*pi/sqrt(3)*R^2/a^2;%b1=2*pi/a*(1-sqrt(3)/3*i);%b2=2*pi/a*2*sqrt(3)/3*i;b1=2*pi/a*[1 -1/sqrt(3) 0];b2=2*pi/a*[0 2/sqrt(3) 0];%n=input('please input n: ');n=5;NumberofPW=(2*n+1)^2;count=1;for x=-n:n,  	for y=-n:n,     	G(count,:)=x*b1+y*b2;      count=count+1;   endendfor x=1:NumberofPW,  	for y=x+1:NumberofPW,   	eps2(x,y)=(epsa-epsb)*2*f*besselj(1,norm(G(x,:)-G(y,:))*R)./(norm(G(x,:)-G(y,:))*R);   	eps2(y,x)=eps2(x,y);      end  	eps2(x,x)=f*epsa+(1-f)*epsb;end   %k1=(0:0.1:1.0)/sqrt(3).*i*2*pi/a;%k2=((0.1:0.1:1.0)./3+1/sqrt(3)*i).*2.*pi./a;%k3=(0.9:-0.1:0).*(1.0/3.0+1/sqrt(3)*i).*2.*pi./a;%-(1/3+1/sqrt(3)*i)*2*pi/a;%k0=[k1 k2 k3];%k0=k1;%k0=(1/3+1/sqrt(3)*i).*2*pi/a;k0=zeros(length(0:0.1:1)+length(0.1:0.1:1)+length(0.9:-0.1:0),3);mm0=length(0:0.1:1);k0(1:mm0,2)=2*pi/a/sqrt(3).*(0:0.1:1)';mm0=length(0:0.1:1);mm1=length(0.1:0.1:1);k0(mm0+1:mm0+mm1,1)=2*pi/a/3.*(0.1:0.1:1)';k0(mm0+1:mm0+mm1,2)=2*pi/a/sqrt(3);mm0=mm0+mm1+1;mm1=length(0.9:-0.1:0);k0(mm0:mm0+mm1-1,1)=2*pi/a/3.*(0.9:-0.1:0)';k0(mm0:mm0+mm1-1,2)=2*pi/a/sqrt(3).*(0.9:-0.1:0)';k0(:,3)=0.0;%k0=[1/3 1/sqrt(3) 9]*2*pi/a;counter=1;eps2=inv(eps2);for ii=1:size(k0,1),   k=k0(ii,:);   K(:,1)=k(1)+G(:,1);   K(:,2)=k(2)+G(:,2);   K(:,3)=0;   %%%find the unit cell perpendicular to it in xy plane   %%%be sure to deal with the case of modulus(k)=0   %%%NaN in this case   e1=[K(:,2)./modulus(K),-K(:,1)./modulus(K),zeros(length(K),1)];   e1(isnan(e1))=1/sqrt(2); %%%when Kz is not zero,kx,ky is zero choose arbitrary e1   K(:,3)=k(3)+G(:,3);   %%%find the other perpendicular unit cell   %%%be sure to deal with the case of modulus(k)=0   %%%NaN in this case   e2=cross(e1, K);   e2=[e2(:,1)./modulus(e2),e2(:,2)./modulus(e2),e2(:,3)./modulus(e2)];   e2(isnan(e2))=0;   %%%form the equation matrix, it should be 2N by 2N   %%%in this case we will have no TE and TM decoupling   M1=([modulus(K).*e2(:,1),modulus(K).*e2(:,2),modulus(K).*e2(:,3)]*[modulus(K).*e2(:,1),modulus(K).*e2(:,2),modulus(K).*e2(:,3)]').*eps2;   M2=([modulus(K).*e1(:,1),modulus(K).*e1(:,2),modulus(K).*e1(:,3)]*[modulus(K).*e2(:,1),modulus(K).*e2(:,2),modulus(K).*e2(:,3)]').*eps2;   M3=([modulus(K).*e2(:,1),modulus(K).*e2(:,2),modulus(K).*e2(:,3)]*[modulus(K).*e1(:,1),modulus(K).*e1(:,2),modulus(K).*e1(:,3)]').*eps2;   M4=([modulus(K).*e1(:,1),modulus(K).*e1(:,2),modulus(K).*e1(:,3)]*[modulus(K).*e1(:,1),modulus(K).*e1(:,2),modulus(K).*e1(:,3)]').*eps2;   M=[M1,M2;M3,M4];	E=sort(abs(eig(M)));	freq(:,counter)=sqrt(abs(E(1:20))).*a./2./pi;	display(sprintf('calculation of k=[%f,%f,%f] is finished',k(1),k(2),k(3)));	counter=counter+1;endtmpx=1:length(k0);%plot(n,freq,'o'),hold onplot(tmpx,freq*2*pi,'linewidth',2)title('Full band structure of a 2D triangular photonic band structure')xlabel('in-plane wave vector k_/_/')ylabel('wa/2\pic')grid on%axis([0,length(k0)+1,0,1.4])