% ASSP - Students Project 
% Calculate DOS for Tight Binding
%
% Thomas Ganner, Klaus Lang
%-------------------------------------------------------------------------
function dosCubic


N = 4*10^6;
lattice = 2;
[k,a] = genKVectorsCubic(lattice,N);
k = k*2*pi; 
epsilon = 0;
t = 1;

switch(lattice)
  case 1 
    % SC
    disp = 2 .*( cos(k(:,1)) + cos(k(:,2)) + cos(k(:,3)) );
    name = 'SC';
    vertices = [1,1,1;-1,1,1;1,-1,1;1,1,-1;1,-1,-1;-1,1,-1;-1,-1,1;-1,-1,-1];
    vertices = vertices*pi;
    sample_space = [0,0,0;1,0,0;0,1,0;0,0,1;1,1,1;1,1,0;1,0,1;0,1,1]*pi;
    sample_faces = [1,2,7,4;1,3,8,4;1,2,6,3;2,7,5,6;3,8,5,6;4,8,5,7;1,2,6,3];

    
  case 2 
    % FCC
    disp = 4 .*( cos(k(:,1)).*cos(k(:,2)) + cos(k(:,2)).*cos(k(:,3)) + cos(k(:,3)).*cos(k(:,1)));
    name = 'FCC';
    vertices = [perms([0,1,2]);perms([0,-1,2]);perms([0,1,-2]);perms([0,-1,-2])];
    vertices = vertices*0.5*pi;
    sample_space = [2,0,0;0,2,0;0,0,2;perms([0,1,2]);0,0,0]*0.5*pi;
    sample_faces = [1,4,5,nan,nan,nan;5,4,6,9,8,7;2,6,9,nan,nan,nan;...
        3,7,8,nan,nan,nan;10,1,5,7,3,nan;10,2,9,8,3,nan;10,1,4,6,2,nan];
    
  case 3 
    % BCC
    disp = 8 .* cos(k(:,1)) .* cos(k(:,2)) .* cos(k(:,3));     
    name = 'BCC';
    vertices = [[1,1,1];[-1,-1,-1];[-1,1,1];[1,-1,1];[1,1,-1];[-1,-1,1];[-1,1,-1];[1,-1,-1];[2,0,0];[0,2,0];[0,0,2];[-2,0,0];[0,-2,0];[0,0,-2]];
    vertices = vertices*0.5*pi;
    sample_space = [2,0,0;0,2,0;0,0,2;1,1,1;0,0,0]*0.5*pi;
    sample_faces = [1,2,4;1,3,4;2,3,4;1,2,5;1,3,5;2,3,5];
end  

dE = linspace(min(disp),max(disp),200);

DOS = hist(disp,dE);
DOS = DOS/(N/2*(dE(2)-dE(1)));

%----- write to File
fid = fopen(['DOS_',name,'_Tight_Binding_1.txt'],'wt');
fprintf(fid,'%f %f\n',[dE;DOS]);
fclose(fid);
%----- 

%----- Plot Density of states
figure(1)
% Strecken der Kurve
dE = epsilon - t .* dE;
plot(dE,DOS,'r-','linewidth',1.5);
xlabel('E [eV]','interpreter','latex')
ylabel('D(E) [1/(eV atom)]','interpreter','latex');
title(['Density of States for ',name,' lattice: \epsilon = ',num2str(epsilon),', t = ',num2str(t)]);

%plot first brillouin zone-------------------------------------------------
figure(3)
%draw vertice points of the brillouin zone

%norm to 1/a*pi. looks better in the plot
vertices = vertices/pi;
sample_space = sample_space/pi;
k = k/pi;

plot3(vertices(:,1),vertices(:,2),vertices(:,3),'ko','markersize',5,...
'markerfacecolor','r');
hold on
%draw lines between vertice points

min_dist_matrix = [];

for i = 1:size(vertices,1)
    for j = 1:size(vertices,1)
        dist(j) = norm(vertices(j,:)-vertices(i,:));
    end
    dist(dist == 0) = nan;
    dist_min = dist == min(dist);
    positions = find(dist_min);
    dm_size = size(min_dist_matrix,2);
    ppoints = numel(positions);
    if ((ppoints+1)>dm_size && i > 1)
        for rp = 1:i
            min_dist_matrix(rp,dm_size:ppoints+1) = nan;
        end
    end
    min_dist_matrix(i,:) = [i,positions];
    for count = 1:numel(positions)
        P(1,:) = vertices(i,:);
        P(2,:) = vertices(positions(count),:);
        line(P(:,1),P(:,2),P(:,3),'color',[0.3,0.3,0.3],'linewidth',1);
    end
end

%create G-vectors or the reciprocal lattice up to order of one
G = getGVectors(a,1)*2; %not 2*pi only 2 (normalized to 1/a*pi)
%get Faces
Faces = getFaces(min_dist_matrix,G,vertices);

%draw patched version

patch('vertices',vertices,'faces',Faces,'facecolor','k');
patch('vertices',sample_space,'faces',sample_faces,'facecolor','r','linewidth',2);
alpha(0.1);
ks = size(k,1);
dk = N/1000; %draw only 1000 points 
plot3(k([1:dk:ks],1),k([1:dk:ks],2),k([1:dk:ks],3),'.');
xlabel('k$_x$ [1/(a$\cdot$pi)]','interpreter','latex')
ylabel('k$_y$ [1/(a$\cdot$pi)]','interpreter','latex')
zlabel('k$_z$ [1/(a$\cdot$pi)]','interpreter','latex')
title(['First brillouin zone of a ',name,' Crystal']);
axis equal

end

% This functions are needed for the visualisation of the brilloin zone. 
function G = getGVectors(a,order)
    pot = [0:order,-order:-1];
    count = 1;
    for i=1:size(pot,2)
        for j = 1:size(pot,2)
            for k=1:size(pot,2)
                G(count,:) = a(1,:)*pot(i) + a(2,:)*pot(j) + a(3,:)*pot(k);
                count = count+1;
            end
        end
    end
    G(1,:) = [];
end


function Faces = getFaces(min_dist_matrix,G_vectors,vertices)
face_id = 0;
Faces = [];

for G_c = 1:size(G_vectors,1)
    face_set = [];
    for v_c = 1:size(vertices,1)
        Bragg1 = (G_vectors(G_c,:)/2)*transpose(vertices(v_c,:));
        Bragg2 = (G_vectors(G_c,:)/2)*transpose(G_vectors(G_c,:)/2);
        if Bragg1 == Bragg2
            face_set = vertcat(face_set,v_c);
        end
    end
    N_fs = numel(face_set);
    N_F = size(Faces,2);
    missing = 0;
    if (N_fs > N_F && face_id ~= 1)
        for rec = 1:face_id
            Faces(rec,N_F:N_fs) = nan; 
        end
    end
    missing = N_F-N_fs;
    if numel(face_set)
        for pch = 1:(numel(face_set)-1)
            vertice_id = face_set(pch);
            position2 = [];
            pch2 = 1;
            while( pch2 <= size(min_dist_matrix,2) && isempty(position2))
                position2 = find(face_set(pch+1:end) == min_dist_matrix(vertice_id,pch2));
                position2 = position2+pch;
                pch2 = pch2 + 1;
            end
            position1 = pch;
            face_set = changePosition(face_set,position1,position2);
        end
        face_id = face_id + 1;
        Faces(face_id,:) = horzcat(transpose(face_set),nan(1,missing));
    end
end
end

function out = changePosition(in,p1,p2)
    if p1 ~= p2
        dummy = in(p2);
        in(p2) = in(p1+1);
        in(p1+1) = dummy;
        out = in;
    else
        out = in; 
    end

end

