findmode.m

This will solve 2d photonic crystal modes by plane wave method

% this program is designed for solving the two dimensional photonic structure
% we are using the parameters from one of the paper
clear
warning off
epsa=1;
epsb=13;
%epsb=1.46^2; %%%Silica glass air system
a=1.0;
%f=0.7;
%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;
   end
end

for 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)=9.0;

%band=zeros(3,100);
bandcount=1;
eps2=inv(eps2);
for xx=0:0.2:6,
   k0(:,3)=xx;
	counter=1;
	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
   %%%to save computing time, use absk=modulus(K)
   absk=modulus(K);
   e1=[K(:,2)./absk,-K(:,1)./absk,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
   absk=modulus(K);
   e2=cross(e1, K);
   e2=[e2(:,1)./absk,e2(:,2)./absk,e2(:,3)./absk];
   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=([absk.*e2(:,1),absk.*e2(:,2),absk.*e2(:,3)]*[absk.*e2(:,1),absk.*e2(:,2),absk.*e2(:,3)]').*eps2;
   M2=([absk.*e1(:,1),absk.*e1(:,2),absk.*e1(:,3)]*[absk.*e2(:,1),absk.*e2(:,2),absk.*e2(:,3)]').*eps2;
   M3=([absk.*e2(:,1),absk.*e2(:,2),absk.*e2(:,3)]*[absk.*e1(:,1),absk.*e1(:,2),absk.*e1(:,3)]').*eps2;
   M4=([absk.*e1(:,1),absk.*e1(:,2),absk.*e1(:,3)]*[absk.*e1(:,1),absk.*e1(:,2),absk.*e1(:,3)]').*eps2;
   M=[M1,M3;M2,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;
   end
   %%%looking for the band gap and save them into band
   display(sprintf('Searching bandgap for beta=%f',xx));
   for ii=1:5,
      if min(freq(ii+1,:)) > max(freq(ii,:))+0.005   %%find a band gap, save these values
         band(1,bandcount)=xx;
         band(2,bandcount)=max(freq(ii,:));
         band(3,bandcount)=min(freq(ii+1,:));
         band(4,bandcount)=ii;
         display(sprintf('find a band (%f %f %f):',band(1,bandcount),band(2,bandcount),band(3,bandcount)));
         bandcount=bandcount+1;
      end
   end
end   
%%%now plotting the variation of the band gap, divide the array data according to the band number
ab1=band;
count1=1;
while length(band)>0,
   ab(:,1,count1)=band(:,1);
   band(:,1)=[];
   count2=1;
   ii=1;
   while ii<=size(band,2),
      %display(sprintf('ii=%f band(4,ii)=%f',ii,band(4,ii)));
      if band(4,ii)==ab(4,1,count1),
         count2=count2+1;
         ab(:,count2,count1)=band(:,ii);
         band(:,ii)=[];
		else
   		ii=ii+1;
		end
   end
   count1=count1+1;
end

plot(ab(1,:,1),ab(2,:,1)*2*pi,'+-',ab(1,:,1),ab(3,:,1)*2*pi,'o-')
hold on
plot(ab(1,1:6,2),ab(2,1:6,2)*2*pi,'+-',ab(1,1:6,2),ab(3,1:6,2)*2*pi,'o-')
plot(ab(1,1:6,3),ab(2,1:6,3)*2*pi,'+-',ab(1,1:6,3),ab(3,1:6,3)*2*pi,'o-')
plot(ab(1,1:5,4),ab(2,1:5,4)*2*pi,'+-',ab(1,1:5,4),ab(3,1:5,4)*2*pi,'o-')
plot(ab(1,1:4,5),ab(2,1:4,5)*2*pi,'+-',ab(1,1:4,5),ab(3,1:4,5)*2*pi,'o-')

%plot(ab(1,1:3,3),ab(2,1:3,3)*2*pi,'+-',ab(1,1:3,3),ab(3,1:3,3)*2*pi,'o-')
xx=7:15;
plot(xx,xx,'r--','linewidth',1);
nav=sqrt(epsb)+f*(sqrt(epsa)-sqrt(epsb));
plot(k0(3,:),k0(3,:)/nav,'r--','linewidth',1);
plot(xx,xx/sqrt(epsb),'r:','linewidth',1);
hold off

title('band gap variation of a 2D triangular SiO2-Air PBG')
xlabel('wave vector k_za')
ylabel('ka or wa/c')
grid on
%axis([0,length(k0)+1,0,1.4])
neff=sqrt(f*epsa+(1-f)*epsb);
alpha1=asin(ab1(1,:)/2/pi/neff./ab1(2,:));
alpha2=asin(ab1(1,:)/2/pi/neff./ab1(3,:));
plot(alpha1*180/pi,ab1(2,:),'o-',alpha2*180/pi,ab1(3,:),'+-')
xlabel('\alpha [degree]')
ylabel('Normalized f: a/\lambda')
axis([0 60 0.38 0.52])