pbs3dfcc.m

Calculates photonic band structure for either an fcc or diamond lattice of dielectric spheres of die

%================================================================
% Name:           PBS3D.m - 3D Photonic Band Structure Code
% Author:         Chun-Feng Lay 
% Data:           04/15/2004
% Description:    Calculates photonic band structure for either
%                 an fcc or diamond lattice of dielectric spheres
%                 of dielectric constant epsilon_a in a dielectric 
%                 background of dielectric constant epsilon_b
%================================================================
clear;
clc;
tic;

%================================================================
% Define real lattice vectors
%================================================================
FCC=0;
%DIAMOND=1;

a=1;
a1=a*[0 0.5 0.5];
a2=a*[0.5 0 0.5];
a3=a*[0.5 0.5 0];

%================================================================
%  Dielectric lattice parameters
%================================================================
epsa=1.0;
epsb=13.0;
epsa1=1;
epsa2=1;

if (FCC==0)
   Pf=(pi/6)*sqrt(2);
else %(DIAMOND==1)
   Pf=(pi/16)*sqrt(3);
end

Vu=dot(a1,cross(a2,a3));

if (FCC==0)
   Rs=(3*Pf*Vu/(4*pi))^(1/3);
else %(DIAMOND==1)
   Rs=(3*Pf*Vu/(8*pi))^(1/3);
end

%===============================================================
%  Generate reciprocal lattice
%===============================================================
ra1=(2*pi/a)*cross(a2,a3)/dot(a1,cross(a2,a3));
ra2=(2*pi/a)*cross(a3,a1)/dot(a1,cross(a2,a3));
ra3=(2*pi/a)*cross(a1,a2)/dot(a1,cross(a2,a3));

%==============================================================
%  Determine reciprocal vectors to use
%  (Spherical region in k-space bounded by Gmax)
%==============================================================
Nrcube=20;
Gmax=0.183*Nrcube*norm(ra1);
G=zeros((2*Nrcube+1)^3,3);
i=1;

for l=-Nrcube:Nrcube
    for m=-Nrcube:Nrcube
        for n=-Nrcube:Nrcube
        Glmn=l*ra1+m*ra2+n*ra3;
        
        if (norm(Glmn)<Gmax)
            G(i,:)=Glmn;
            i=i+1;
        end
        end
    end
end

NG=i-1;
G=G(1:NG,:);

%====================================================================
%  Generate k-space f(Gi-Gj) values for every i,j=1 to NG
%====================================================================
f=zeros(NG,NG);

for i=1:NG
    for j=1:NG
        Gij=norm(G(i,:)-G(j,:));
        if (Gij==0)
            f(i,j)=(1/epsb)+Pf*(1/(2*epsa1)+1/(2*epsa2)-1/epsb);
        else
            if (FCC==0)
            f(i,j)=3*Pf*(1/epsb-1/epsa)*(Gij*Rs*cos(Gij*Rs)-sin(Gij*Rs))/((Gij*Rs)^3);
        else %(DIAMOND==1)
            tau=a*[1/8 1/8 1/8];
            f(i,j)=3*Pf*(cos(dot(G(i,:)-G(j,:),tau))/epsb...
                   -exp(-sqrt(-1)*dot(G(i,:)-G(j,:),tau))/(2*epsa1)...
                   -exp(sqrt(-1)*dot(G(i,:)-G(j,:),tau))/(2*epsa2))...
                   *(Gij*Rs*cos(Gij*Rs)-sin(Gij*Rs))/((Gij*Rs)^3);
        end
    end
    end
end

%====================================================================
%  Define reciprocal lattice vectors which k-space trajectory
%  The k-space trajectory follows the irreducible portion of the
%  Brillouin zone
%===================================================================
X=(2*pi/a)*[1 0 0];
U=(2*pi/a)*[1 0.25 0.25];
L=(2*pi/a)*[1/2 1/2 1/2];
T=(2*pi/a)*[1e-4 0 0];
W=(2*pi/a)*[1 1/2 0];
K=(2*pi/a)*[0.75 0.75 0]; 

%==================================================================
%  Get eigenfrequences for each k-space point on the trajectory
%==================================================================
Theta=zeros(2*NG,2*NG);
Nkpoints=10;
stepsize=0:1/(Nkpoints-1):1;
XU_eig=zeros(Nkpoints,2*NG);
UL_eig=zeros(Nkpoints,2*NG);
LT_eig=zeros(Nkpoints,2*NG);
TX_eig=zeros(Nkpoints,2*NG);
XW_eig=zeros(Nkpoints,2*NG);
WK_eig=zeros(Nkpoints,2*NG);
ei=zeros(1,3);
ej=zeros(1,3);

for n=1:Nkpoints
fprintf(['\n k-point:',int2str(n),'of',int2str(Nkpoints),'.\n']);

%=================================================================
%  Calculate eigenfrequencies for the k-space step along the
%  XU trajectory
%=================================================================
XU_step=stepsize(n)*(U-X)+X;

%=================================================================
%  Get the polarization vectors for the NG plane waves for the
%  current step along the XU trajectory
%=================================================================
for i=1:NG
    Ki=XU_step+G(i,:);
    [e1,e2]=get_polarization_vectors(Ki);
    e_polar_vec(1,i,:)=e1;
    e_polar_vec(2,i,:)=e2;
end

%================================================================
%  Determine the off-diag elements of Theta matrix for XU_step
%===============================================================
for i=1:(2*NG-1)
    for j=(i+1):(2*NG)
        iG=ceil(i/2);
        jG=ceil(j/2);
        ie=mod(i,2)+mod(mod(i,2)+2,3);
        je=mod(j,2)+mod(mod(j,2)+2,3);
        
        ei(1:3)=e_polar_vec(ie,iG,:);
        ej(1:3)=e_polar_vec(je,jG,:);
        Ki=XU_step+G(iG,:);
        Kj=XU_step+G(jG,:);
        
        Theta(i,j)=f(iG,jG)*dot(cross(Ki,ei),cross(Kj,ej));
        Theta(j,i)=conj(Theta(i,j));
    end
end

%===============================================================
%  Determine the diagonal elements of Theta matrix for XU_step
%===============================================================
for i=1:(2*NG)
    iG=ceil(i/2);
    ie=mod(i,2)+mod(mod(i,2)+2,3);
    
    ei(1:3)=e_polar_vec(ie,iG,:);
    Ki=XU_step+G(iG,:);
        
    Theta(i,i)=f(iG,iG)*norm(cross(Ki,ei))^2;
    end
 XU_eig(n,:)=sort(sqrt(eig(Theta))).';
 
%============================================================
%  Calculate eigenfrequencies for the k-space step along the
%  UL trajectory
%============================================================
UL_step=stepsize(n)*(L-U)+U;

%=================================================================
%  Get the polarization vectors for the NG plane waves for the
%  current step along the UL trajectory
%=================================================================
for i=1:NG
    Ki=UL_step+G(i,:);
    [e1,e2]=get_polarization_vectors(Ki);
    e_polar_vec(1,i,:)=e1;
    e_polar_vec(2,i,:)=e2;
end

%================================================================
%  Determine the off-diag elements of Theta matrix for UL_step
%===============================================================
for i=1:(2*NG-1)
    for j=(i+1):(2*NG)
        iG=ceil(i/2);
        jG=ceil(j/2);
        ie=mod(i,2)+mod(mod(i,2)+2,3);
        je=mod(j,2)+mod(mod(j,2)+2,3);
        
        ei(1:3)=e_polar_vec(ie,iG,:);
        ej(1:3)=e_polar_vec(je,jG,:);
        Ki=UL_step+G(iG,:);
        Kj=UL_step+G(jG,:);
        
        Theta(i,j)=f(iG,jG)*dot(cross(Ki,ei),cross(Kj,ej));
        Theta(j,i)=conj(Theta(i,j));
    end
end

%===============================================================
%  Determine the diagonal elements of Theta matrix for UL_step
%===============================================================
for i=1:(2*NG)
    iG=ceil(i/2);
    ie=mod(i,2)+mod(mod(i,2)+2,3);
    
    ei(1:3)=e_polar_vec(ie,iG,:);
    Ki=UL_step+G(iG,:);
        
    Theta(i,i)=f(iG,iG)*norm(cross(Ki,ei))^2;
    end
 UL_eig(n,:)=sort(sqrt(eig(Theta))).';
 
%============================================================
%  Calculate eigenfrequencies for the k-space step along the
%  LT trajectory
%============================================================
LT_step=stepsize(n)*(T-L)+L;

%=================================================================
%  Get the polarization vectors for the NG plane waves for the
%  current step along the LT trajectory
%=================================================================
for i=1:NG
    Ki=LT_step+G(i,:);
    [e1,e2]=get_polarization_vectors(Ki);
    e_polar_vec(1,i,:)=e1;
    e_polar_vec(2,i,:)=e2;
end

%================================================================
%  Determine the off-diag elements of Theta matrix for LT_step
%===============================================================
for i=1:(2*NG-1)
    for j=(i+1):(2*NG)
        iG=ceil(i/2);
        jG=ceil(j/2);
        ie=mod(i,2)+mod(mod(i,2)+2,3);
        je=mod(j,2)+mod(mod(j,2)+2,3);
        
        ei(1:3)=e_polar_vec(ie,iG,:);
        ej(1:3)=e_polar_vec(je,jG,:);