% this program is designed for solving the two dimensional photonic structure
clear
warning off
epsa=1;
epsb=11.56;
%epsb=11.58;
a=1.0;
%f=0.8358;
%R=sqrt(sqrt(3).*f/(2*pi)*a^2);
R=0.48*a;
i=sqrt(-1);
f=pi*R^2/a^2;
a1=a;
a2=a*i;
b1=2*pi/a;
b2=2*pi/a*i;
n=input('please input n: ');
NumberofPW=(2*n+1)^2;

count=1;
for x=-n:n,
  	for y=-n:n,
     	G(count)=x*b1+y*b2;
      count=count+1;
   end
end

for x=1:NumberofPW,
  	for y=x+1:NumberofPW,
   	eps2(x,y)=(epsa-epsb)*2*f*besselj(1,abs(G(x)-G(y))*R)./(abs(G(x)-G(y))*R);
   	eps2(y,x)=eps2(x,y);   
   end
  	eps2(x,x)=f*epsa+(1-f)*epsb;
end
   
k1=(0.0:0.1:1.0)*pi./a;
k2=(0.1:0.1:1.0)*pi./a*i+pi/a;
k3=(0.9:-0.1:0).*(1.0+1*i).*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;
counter=1;
eps2=inv(eps2);
for ii=1:length(k0),
   k=k0(ii);
   M=abs(k+G.')*abs(k+G).*(eps2);%TM mode
   %M=(real(k+G.')*real(k+G)+imag(k+G.')*imag(k+G)).*(eps2);%TE mode

	E=sort(abs(eig(M)));
	freq(:,counter)=sqrt(abs(E(1:20))).*a./2./pi;
	%display(sprintf('calculation of k=%f+%fi is finished',real(k),imag(k)));
	counter=counter+1;
end
tmpx=0:length(k0)-1;
%plot(n,freq,'o'),hold on
plot(tmpx,freq,'linewidth',1)
title('TM:Band structure of a 2D square photonic band structure')
xlabel('wave vector')
ylabel('\omegaa/2\pic')
grid on
axis([0,length(k0)-1,0,0.8])