function PlotFreeElectronEnergyInReducedZone

pi = 4.0*atan(1.0);
N=256;
x=(1:N);
u=0.5*x/N;
h = 0;
k = 0;
l = 0;
uu = 0.5;
LowestBandEnergy =  (uu-h+k+l).^2+(uu+h-k+l).^2+(uu+h+k-l).^2
ncnt=0;
% markers for the [hkl] 
for h=-1:1
    for k=-1:1
        for l=-1:1
            uu = 0;
            offset = (uu-h+k+l).^2+(uu+h-k+l).^2+(uu+h+k-l).^2;
            % plot only if energy is less than smallest energy unit
            if (offset <= 6*LowestBandEnergy) 
                fprintf(1,'%d %d %d\n',h,k,l);
                E=(u-h+k+l).^2+(u+h-k+l).^2+(u+h+k-l).^2;
                plot(u(ncnt*10+40),E(ncnt*10+40),'.','MarkerSize',25,'Color',...
                    hsv2rgb([ncnt/15.0 1.0 1.0]));
                    %[.4*(h+1) .4*(k+1) .5*(l+1)]);
                hold on;
                if (ncnt == 0)
                    ncnt = ncnt+1;
                    lstring = sprintf('%2d %2d %2d\n',h,k,l);
                    Key = lstring;
                else                   
                    ncnt = ncnt + 1;
                    lstring = sprintf('%2d %2d %2d\n',h,k,l);
                    Key = [Key;lstring];
                end
            end
        end
    end
end
Legend(Key,-1); 
% black lines
for h=-1:1
    for k=-1:1
        for l=-1:1
            uu = 0;
            offset = (uu-h+k+l).^2+(uu+h-k+l).^2+(uu+h+k-l).^2;
            % plot only if energy is less than smallest energy unit
            if (offset <= 6*LowestBandEnergy) 
                fprintf(1,'%d %d %d\n',h,k,l);
                E=(u-h+k+l).^2+(u+h-k+l).^2+(u+h+k-l).^2;
                plot(u,E); %,'Color',[.4*(h+1) .4*(k+1) .5*(l+1)]);
                hold on;
            end
        end
    end
end
title(['\bfFree Electron Energy for bcc Brillouin Zone' ],'fontsize',14);
ylabel(['\bfEnergy/E0 [E0=2h^{2}/ma]'],'fontsize',12);
xlabel('\bfu','fontsize',12);
% markers over lines
ncnt = 0;
for h=-1:1
    for k=-1:1
        for l=-1:1
            uu = 0;
            offset = (uu-h+k+l).^2+(uu+h-k+l).^2+(uu+h+k-l).^2;
            % plot only if energy is less than smallest energy unit
            if (offset <= 6*LowestBandEnergy) 
                fprintf(1,'%d %d %d\n',h,k,l);
                E=(u-h+k+l).^2+(u+h-k+l).^2+(u+h+k-l).^2;
                plot(u(ncnt*10+40),E(ncnt*10+40),'.','MarkerSize',25,'Color',...
                    hsv2rgb([ncnt/15.0 1.0 1.0]));
                    %[.4*(h+1) .4*(k+1) .5*(l+1)]);
                hold on;
                ncnt = ncnt+1;
            end
        end
    end
end
saveas(gcf, ['FreeElectronEnergy','.emf']);
