% ASSP - Students Project 
% Calculate DOS for Tight Binding
%
% Thomas Ganner, Klaus Lang
% generate k vectors for SC, FCC, BCC
%-------------------------------------------------------------------------
function [k,A] = genKVectorsCubic(lattice,N)

% N       ... Anzahl der Vektoren
% lattice ... 1 -> SC, 2 -> FCC, 3 -> BCC
              
if nargin < 2, N = 10000; end
if nargin < 1, lattice = 2; end

switch(lattice)
  case 1 
    % SC
    k = rand(N,3)*0.5; 
    a = [1,0,0 ;0,1,0 ;0,0,1]; %komisch aber funktioniert
    A = a;
    return
              
  case 2 
    % FCC - nearest neighbors
    a = [ 1, 1,-1; 1,-1, 1; -1, 1, 1];
    ind = [1,1,1];
    
  case 3 
    % BCC - nearest neighbors
    a = [ 1, 1, 0; 0, 1, 1; 1, 0, 1];
    ind = [1,0,0; 0,1,0; 0,0,1];
    
  otherwise
    k = []; 
    return
end
A = a*0.5;
a = transpose(a)*0.5;
k = zeros(N,3);
num_ind = size(ind,1);

ind = transpose(ind);


miss = N;

h = waitbar(0,'Bitte warten...');
set(h,'Name','Status')
waitbar((N-miss)/N,h,[num2str((N-miss)/N*100,'%.1f'),' % abgeschlossen...'])

while(miss > 0)
  test_k = rand(miss,3)*0.5; 
  bool = ones(miss,1);
  
  norm_test_k = absVector(test_k(:,1),test_k(:,2),test_k(:,3));
  
  for l = 1 : num_ind
    x = a*ind(:,l);
    x1 = test_k(:,1) - x(1);
    x2 = test_k(:,2) - x(2);
    x3 = test_k(:,3) - x(3);
    bool = bool & (absVector(x1,x2,x3) > norm_test_k );
  end
  
  num_new_elements = sum(bool);
  
  k((N-miss+1):(N-miss+num_new_elements),:) = test_k(bool,:);
  miss = miss - num_new_elements;
  
  waitbar((N-miss)/N,h,[num2str((N-miss)/N*100,'%.1f'),' % abgeschlossen...'])
  pause(0.001);
end
close(h)

%----- Plot
% figure(1)
% subplot(1,3,1)
% plot(k(:,1),k(:,2),'r*');
% axis equal
% xlabel('k_X')
% ylabel('k_Y')
% 
% subplot(1,3,2)
% plot(k(:,2),k(:,3),'r*');
% axis equal
% xlabel('k_Y')
% ylabel('k_Z')
% 
% subplot(1,3,3)
% plot(k(:,3),k(:,1),'r*');
% axis equal
% xlabel('k_Z')
% ylabel('k_X')
%-----

%-------------------------------------------------------------------------
function y = absVector(x1,x2,x3)
y = sqrt(x1.^2 + x2.^2 + x3.^2);