pbs3dbcc.m

Calculates photonic band structure for either the bcc lattice of dielectric spheres of dielectric co

%================================================================
% Name:           PBS3Dbcc.m - 3D Photonic Band Structure Code
% Author:         Chun-Feng Lay 
% Data:           05/20/2004
% Description:    Calculates photonic band structure for either
%                 an bcc 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
%================================================================
BCC=0;

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

%================================================================
%  Dielectric lattice parameters
%================================================================
epsa=1;
epsb=11.74;
epsa1=1;
epsa2=1;
Pf=0.74;
Vu=dot(a1,cross(a2,a3));
Rs=(3*Pf*Vu/(4*pi))^(1/3);

%===============================================================
%  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
            f(i,j)=3*Pf*(1/epsb-1/epsa)*(Gij*Rs*cos(Gij*Rs)-sin(Gij*Rs))/((Gij*Rs)^3);
    end
    end
end

%====================================================================
%  Define reciprocal lattice vectors which k-space trajectory
%  The k-space trajectory follows the irreducible portion of the
%  Brillouin zone
%===================================================================
P=(2*pi/a)*[0.5 0.5 0.5];
T=(2*pi/a)*[1e-4 0 0];
H=(2*pi/a)*[0 0 1];
N=(2*pi/a)*[0 0.5 0.5];
P=(2*pi/a)*[0.5 0.5 0.5];

%==================================================================
%  Get eigenfrequences for each k-space point on the trajectory
%==================================================================
Theta=zeros(2*NG,2*NG);
Nkpoints=10;
stepsize=0:1/(Nkpoints-1):1;
PT_eig=zeros(Nkpoints,2*NG);
TH_eig=zeros(Nkpoints,2*NG);
HN_eig=zeros(Nkpoints,2*NG);
NP_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
%  PT trajectory
%=================================================================
PT_step=stepsize(n)*(T-P)+P;

%=================================================================
%  Get the polarization vectors for the NG plane waves for the
%  current step along the XU trajectory
%=================================================================
for i=1:NG
    Ki=PT_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 PT_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=PT_step+G(iG,:);
        Kj=PT_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 PT_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=PT_step+G(iG,:);
        
    Theta(i,i)=f(iG,iG)*norm(cross(Ki,ei))^2;
    end
 PT_eig(n,:)=sort(sqrt(eig(Theta))).';
 
%============================================================
%  Calculate eigenfrequencies for the k-space step along the
%  TH trajectory
%============================================================
TH_step=stepsize(n)*(H-T)+T;

%=================================================================
%  Get the polarization vectors for the NG plane waves for the
%  current step along the TH trajectory
%=================================================================
for i=1:NG
    Ki=TH_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 TH_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=TH_step+G(iG,:);
        Kj=TH_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 TH_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=TH_step+G(iG,:);
        
    Theta(i,i)=f(iG,iG)*norm(cross(Ki,ei))^2;
    end
TH_eig(n,:)=sort(sqrt(eig(Theta))).';
 
%============================================================
%  Calculate eigenfrequencies for the k-space step along the
%  HN trajectory
%============================================================
HN_step=stepsize(n)*(N-H)+H;

%=================================================================
%  Get the polarization vectors for the NG plane waves for the
%  current step along the HN trajectory
%=================================================================
for i=1:NG
    Ki=HN_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 HN_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=HN_step+G(iG,:);
        Kj=HN_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 HN_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=HN_step+G(iG,:);
        
    Theta(i,i)=f(iG,iG)*norm(cross(Ki,ei))^2;
    end
 HN_eig(n,:)=sort(sqrt(eig(Theta))).';

 %============================================================
%  Calculate eigenfrequencies for the k-space step along the
%  NP trajectory
%============================================================
NP_step=stepsize(n)*(P-N)+N;

%=================================================================
%  Get the polarization vectors for the NG plane waves for the
%  current step along the NP trajectory
%=================================================================
for i=1:NG
    Ki=NP_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 NP_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=NP_step+G(iG,:);
        Kj=NP_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 NP_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=NP_step+G(iG,:);
        
    Theta(i,i)=f(iG,iG)*norm(cross(Ki,ei))^2;
 end
 NP_eig(n,:)=sort(sqrt(eig(Theta))).';
end

fprintf('\n  Total calculation time:  %d  sec. \n',toc);
save pbs3Dbccdata;

%=================================================================
%  Trajectory lengths different, account for this when plotting
%=================================================================
kaxis=0;
PTaxis=kaxis:norm(P-T)/(Nkpoints-1):(kaxis+norm(P-T));
kaxis=kaxis+norm(P-T);
THaxis=kaxis:norm(T-H)/(Nkpoints-1):(kaxis+norm(T-H));
kaxis=kaxis+norm(T-H);
HNaxis=kaxis:norm(H-N)/(Nkpoints-1):(kaxis+norm(H-N));
kaxis=kaxis+norm(H-N);
NPaxis=kaxis:norm(N-P)/(Nkpoints-1):(kaxis+norm(N-P));
kaxis=kaxis+norm(N-P);

%=========================================================================
%  Plot 3D Photonic Band Structure (PBS)
%=========================================================================
Ntraject=6;
EigFreq=zeros(Ntraject*Nkpoints,1);

figure(1);
hold on;

Nk=Nkpoints;

for k=1:NG
    for i=1:Nkpoints
            EigFreq(i+0*Nk)=(PT_eig(i,k))/(2*pi/a);
            EigFreq(i+1*Nk)=(TH_eig(i,k))/(2*pi/a);
            EigFreq(i+2*Nk)=(HN_eig(i,k))/(2*pi/a);
            EigFreq(i+3*Nk)=(NP_eig(i,k))/(2*pi/a);
        end
        if (BCC==0)
        plot(PTaxis(1:Nk),EigFreq(1+0*Nk:1*Nk),'b',...
             THaxis(1:Nk),EigFreq(1+1*Nk:2*Nk),'b',...
             HNaxis(1:Nk),EigFreq(1+2*Nk:3*Nk),'b',...
             NPaxis(1:Nk),EigFreq(1+3*Nk:4*Nk),'b');
     end
 end
 
 grid on;
 hold off;
 
 titlestr=['PBS for BCC '];
     
 titlestr=[titlestr '(NG=',int2str(NG),...
                    ',Pf=',num2str(Pf),...
                    ',epsa=',num2str(epsa),...
                    ',epsb=',num2str(epsb),')'];
 title (titlestr);
 
 xlabel ('k-Space');
 ylabel ('Frequency (wa/c2\pi)');
 
 axis ([0 NPaxis(Nkpoints) 0 1]);
 set (gca,'XTick',[PTaxis(1)...
                   PTaxis(Nk)...
                   THaxis(Nk)...
                   HNaxis(Nk)...
                   NPaxis(Nk)]);
 xtixlabel=strvcat('P', 'T', 'H', 'N', 'P');
 set(gca,'XTickLabel',xtixlabel);