diamond3d.m

For 3d Photonic crystals with diamond structure. Calculation of photonic band gaps etc. using matlab

%   3D Diamond Lattice
%   Shangping Guo
%	 Reference data: K M Ho et al, PRL,65-25
%   Also see S G Johnson's webpage on his software
%   The k0 along the 1BZ is obtained from his website.
%   This program yields accurate results,

clear
warning off
%epsa=11.56;
epsa=3.6^2;
epsb=1; %
a=1.0;
tic;
R=sqrt(3)/8*a;%a is the length of the cubic unit cell
i=sqrt(-1);
%define the lattice vectors
a1=[0 1 1]*a/2;
a2=[1 0 1]*a/2;
a3=[1 1 0]*a/2;
%calculate the reciprocal lattice vectors
volcell=dot(a1,cross(a2,a3));
b1=2*pi*cross(a2,a3)/volcell;
b2=2*pi*cross(a3,a1)/volcell;
b3=2*pi*cross(a1,a2)/volcell;
%n=input('please input n: ');
f=2*(4*pi*R^3/3)/volcell;
n=input('Input n:');

%n=3;
NumberofPW=(2*n+1)^3;
display('Forming G vectors...');
count=1;
for x=-n:n,
   for y=-n:n,
      for z=-n:n,
     	G(count,:)=x*b1+y*b2+z*b3;
      count=count+1;
     end
   end
end
display('Calculating eps(G-G)...');
r0=[0.125 0.125 0.125];%*[a1;a2;a3]/modulus(a1);
for x=1:NumberofPW,
  	for y=x+1:NumberofPW,
      tmpG=norm(G(x,:)-G(y,:));
      eps2(x,y)=(epsa-epsb)*3*f*(sin(tmpG*R)-tmpG*R*cos(tmpG*R))/(tmpG*R)^3*cos(dot(G(x,:)-G(y,:),r0));
	%This is problemetic, for we add two spheres together,the f alrady contains it
   	eps2(y,x)=eps2(x,y);   
   end
  	eps2(x,x)=f*epsa+(1-f)*epsb;
end
eps2=inv(eps2);
%forming the list of k-vectors on the edge of the 1BZ
display('Forming k-vector array...');
k1=interpolate([0 0.5 0.5]*[b1;b2;b3],[0 0.625 0.375]*[b1;b2;b3],4,0);
k2=interpolate([0 0.625 0.375]*[b1;b2;b3],[0 0.5 0]*[b1;b2;b3],4,0);
k3=interpolate([0 0.5 0]*[b1;b2;b3],[0 0 0]*[b1;b2;b3],4,0);
k4=interpolate([0 0 0]*[b1;b2;b3],[0 0.5 0.5]*[b1;b2;b3],4,0);
k5=interpolate([0 0.5 0.5]*[b1;b2;b3],[0.25 0.5 0.75]*[b1;b2;b3],4,0);
k6=interpolate([0.25 0.5 0.75]*[b1;b2;b3],[0.375 0.375 0.75]*[b1;b2;b3],4,1);
k0=[k1;k2;k3;k4;k5;k6];
display('Calculating eigen frequency...');
counter=1;
for ii=1:size(k0,1),
   tic;
   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
   mK=modulus(K);
   M1=([mK.*e2(:,1),mK.*e2(:,2),mK.*e2(:,3)]*[mK.*e2(:,1),mK.*e2(:,2),mK.*e2(:,3)]').*eps2;
   M2=([mK.*e1(:,1),mK.*e1(:,2),mK.*e1(:,3)]*[mK.*e2(:,1),mK.*e2(:,2),mK.*e2(:,3)]').*eps2;
   M3=([mK.*e2(:,1),mK.*e2(:,2),mK.*e2(:,3)]*[mK.*e1(:,1),mK.*e1(:,2),mK.*e1(:,3)]').*eps2;
   M4=([mK.*e1(:,1),mK.*e1(:,2),mK.*e1(:,3)]*[mK.*e1(:,1),mK.*e1(:,2),mK.*e1(:,3)]').*eps2;
   M=[M1,M3;M2,M4];

	E=sort(abs(eig(M)));
	freq(:,counter)=sqrt(abs(E(1:10))).*a./2./pi;
%	display(sprintf('calculation of k=[%f,%f,%f] is finished',k(1),k(2),k(3)));
	display(sprintf('%d/%d is completed:%f s used',counter,length(k0),toc));
   counter=counter+1;
   save diafreq3 freq;%long time may needed
end
tmpx=1:length(k0);
plot(tmpx,freq,'o-','linewidth',2)
title('Full band structure of a 3D diamond photonic crystal')
xlabel('wave vector k')
ylabel('wa/2\pic')
grid on
%axis([0,length(k0)+1,0,1.4])
toc