3dpbg.m

计算3维光子晶体的能带结构

% 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,
function threeDstation
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);
%k1=interp2([0 0.5 0.5]*[b1;b2;b3],[0 0.625 0.375]*[b1;b2;b3],4,0);
%k2=interp2([0 0.625 0.375]*[b1;b2;b3],[0 0.5 0]*[b1;b2;b3],4,0);
%k3=interp2([0 0.5 0]*[b1;b2;b3],[0 0 0]*[b1;b2;b3],4,0);
%k4=interp2([0 0 0]*[b1;b2;b3],[0 0.5 0.5]*[b1;b2;b3],4,0);
%k5=interp2([0 0.5 0.5]*[b1;b2;b3],[0.25 0.5 0.75]*[b1;b2;b3],4,0);
%k6=interp2([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])
function varray=interpolate(vbegin,vend,n,inclend)
varray=zeros(n+1,3);
step=(vend-vbegin)/n;
for i=1:3,
   if step(i) ==0,
      varray(:,i)=vbegin(i);
   else
      varray(:,i)=(vbegin(i):step(i):vend(i))';
   end
end
if inclend==0
   varray(n+1,:)=[];
end

function z=modulus(x)
%z=sqrt(x(:,1).^2+x(:,2).^2+x(:,3).^2);
[m n]=size(x);
if n==3,
	z=sqrt(x(:,1).^2+x(:,2).^2+x(:,3).^2);
else
   z=sqrt(x(1,:).^2+x(2,:).^2+x(3,:).^2);
end


toc