pbs3dfcc.m

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

        Ki=LT_step+G(iG,:);
        Kj=LT_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 LT_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=LT_step+G(iG,:);
        
    Theta(i,i)=f(iG,iG)*norm(cross(Ki,ei))^2;
    end
 LT_eig(n,:)=sort(sqrt(eig(Theta))).';

 %============================================================
%  Calculate eigenfrequencies for the k-space step along the
%  TX trajectory
%============================================================
TX_step=stepsize(n)*(X-T)+T;

%=================================================================
%  Get the polarization vectors for the NG plane waves for the
%  current step along the TX trajectory
%=================================================================
for i=1:NG
    Ki=TX_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 TX_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=TX_step+G(iG,:);
        Kj=TX_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 TX_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=TX_step+G(iG,:);
        
    Theta(i,i)=f(iG,iG)*norm(cross(Ki,ei))^2;
    end
 TX_eig(n,:)=sort(sqrt(eig(Theta))).';
 
 %============================================================
%  Calculate eigenfrequencies for the k-space step along the
%  XW trajectory
%============================================================
XW_step=stepsize(n)*(W-X)+X;

%=================================================================
%  Get the polarization vectors for the NG plane waves for the
%  current step along the XW trajectory
%=================================================================
for i=1:NG
    Ki=XW_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 XW_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=XW_step+G(iG,:);
        Kj=XW_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 XW_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=XW_step+G(iG,:);
        
    Theta(i,i)=f(iG,iG)*norm(cross(Ki,ei))^2;
    end
 XW_eig(n,:)=sort(sqrt(eig(Theta))).';
 
 %============================================================
%  Calculate eigenfrequencies for the k-space step along the
%  WK trajectory
%============================================================
WK_step=stepsize(n)*(K-W)+W;

%=================================================================
%  Get the polarization vectors for the NG plane waves for the
%  current step along the WK trajectory
%=================================================================
for i=1:NG
    Ki=WK_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 WK_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=WK_step+G(iG,:);
        Kj=WK_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 WK_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=WK_step+G(iG,:);
        
    Theta(i,i)=f(iG,iG)*norm(cross(Ki,ei))^2;
    end
 
 WK_eig(n,:)=sort(sqrt(eig(Theta))).';
end
fprintf('\n  Total calculation time:  %d  sec. \n',toc);
save pbs3Ddata;

%=================================================================
%  Trajectory lengths different, account for this when plotting
%=================================================================
kaxis=0;
XUaxis=kaxis:norm(X-U)/(Nkpoints-1):(kaxis+norm(X-U));
kaxis=kaxis+norm(X-U);
ULaxis=kaxis:norm(U-L)/(Nkpoints-1):(kaxis+norm(U-L));
kaxis=kaxis+norm(U-L);
LTaxis=kaxis:norm(L-T)/(Nkpoints-1):(kaxis+norm(L-T));
kaxis=kaxis+norm(L-T);
TXaxis=kaxis:norm(T-X)/(Nkpoints-1):(kaxis+norm(T-X));
kaxis=kaxis+norm(T-X);
XWaxis=kaxis:norm(X-W)/(Nkpoints-1):(kaxis+norm(X-W));
kaxis=kaxis+norm(X-W);
WKaxis=kaxis:norm(W-K)/(Nkpoints-1):(kaxis+norm(W-K));
kaxis=kaxis+norm(W-K);

%=========================================================================
%  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)=(XU_eig(i,k))/(2*pi/a);
            EigFreq(i+1*Nk)=(UL_eig(i,k))/(2*pi/a);
            EigFreq(i+2*Nk)=(LT_eig(i,k))/(2*pi/a);
            EigFreq(i+3*Nk)=(TX_eig(i,k))/(2*pi/a);
            EigFreq(i+4*Nk)=(XW_eig(i,k))/(2*pi/a);
            EigFreq(i+5*Nk)=(WK_eig(i,k))/(2*pi/a);
        end
        if (FCC==0)
        plot(XUaxis(1:Nk),EigFreq(1+0*Nk:1*Nk),'b',...
             ULaxis(1:Nk),EigFreq(1+1*Nk:2*Nk),'b',...
             LTaxis(1:Nk),EigFreq(1+2*Nk:3*Nk),'b',...
             TXaxis(1:Nk),EigFreq(1+3*Nk:4*Nk),'b',...
             XWaxis(1:Nk),EigFreq(1+4*Nk:5*Nk),'b',...
             WKaxis(1:Nk),EigFreq(1+5*Nk:6*Nk),'b'); 
        else %(DIAMOND=1)
        plot(XUaxis(1:Nk),EigFreq(1+0*Nk:1*Nk),'b',...
             ULaxis(1:Nk),EigFreq(1+1*Nk:2*Nk),'b',...
             LTaxis(1:Nk),EigFreq(1+2*Nk:3*Nk),'b',...
             TXaxis(1:Nk),EigFreq(1+3*Nk:4*Nk),'b',...
             XWaxis(1:Nk),EigFreq(1+4*Nk:5*Nk),'b',...
             WKaxis(1:Nk),EigFreq(1+5*Nk:6*Nk),'b'); 
     end
 end
 
 grid on;
 hold off;
 
 if (FCC==0)
     titlestr=['PBS for FCC '];
 else %(DIAMOND==1)
     titlestr=['PBS for DIAMOND '];
 end     
 titlestr=[titlestr '(NG=',int2str(NG),...
                    ',Pf=',num2str(Pf),...
                    ',epsa=',num2str(epsa),...
                    ',epsb=',num2str(epsb),')'];
 title (titlestr);
 
 xlabel ('k-Space');
 ylabel ('Normalized frequency (wa/2\pic)');
 
 axis ([0 WKaxis(Nkpoints) 0 1]);
 set (gca,'XTick',[XUaxis(1)...
                   XUaxis(Nk)...
                   ULaxis(Nk)...
                   LTaxis(Nk)...
                   TXaxis(Nk)...
                   XWaxis(Nk)...
                   WKaxis(Nk)]);
 xtixlabel=strvcat('X', 'U', 'L', 'T', 'X', 'W', 'K');
 set(gca,'XTickLabel',xtixlabel);