%BRILLOUIN ZONES
%VORONOI 3d
%Poms Johannes
%Based on QHULL DEMO
%http://www.mathworks.com/products/demos/shipping/matlab/qhulldemo.html

clear all;
close all;
clc;

disp('+++++++++++++++++++');
disp('  BRILLOUIN ZONES  ');
disp('(c)  Poms 2009     ');
disp('+++++++++++++++++++');

a=input(['\n','General triclinic lattice','\n','a: ']);
if isempty(a), a=1.2; end
b=input('b: ');
if isempty(b), b=2; end
c=input('c: ');
if isempty(c), c=1.5; end
alpha=input('alpha: ');
if isempty(alpha), alpha=72; end
beta=input('beta: ');
if isempty(beta), beta=80; end
gamma=input('gamma: ');
if isempty(gamma), gamma=110; end
animation=input('Make an animation and save the movie as gif yes/no: ','s');
if isempty(animation), animation='y'; end

alpha=alpha*pi/180;
beta=beta*pi/180;
gamma=gamma*pi/180;

w1 = cos(alpha) - cos(beta) * cos(gamma)/(sin(beta) * sin(gamma));
w2 = sin(gamma)^2 - cos(beta)^2 - cos(alpha)^2 ...
    + 2*cos(alpha) *cos(beta) * cos(gamma);
w2 = sqrt(w2) / (sin(beta) * sin(gamma));

T1 = [a, 0, 0];
T2 = [b * cos(gamma), b	* sin(gamma), 0];
T3 = [c * cos(beta), c * w1*sin(beta), c*w2*sin(beta)];

spat=cross(T1,T2)*T3';

G1=2*pi*cross(T2,T3)/spat;
G2=2*pi*cross(T3,T1)/spat;
G3=2*pi*cross(T1,T2)/spat;

X = zeros(125,3);

l=1;
for i=-2:2;
    for j=-2:2;
        for k=-2:2;
            X(l,:)=i*G1+j*G2+k*G3;
            l=l+1;
        end
    end
end


cla reset; hold on

% Compute Voronoi diagram.
[c,v] = voronoin(X);

nx = c(v{63},:);
tri = convhulln(nx);

% Plot the Voronoi diagram.

if (animation ~= 'y')
    fh=figure(1);
    plot3(X(:,1),X(:,2),X(:,3),'b.','markersize',10);
    for i = 1:size(tri,1)
        patch(nx(tri(i,:),1),nx(tri(i,:),2),nx(tri(i,:),3),i,'FaceAlpha',0.85);
    end
    % Modify the view.
    view(3);
    title('1st Brillouin zone');
    axis equal tight off  vis3d
    grid on;
    camzoom(2);
    rotate3d on
else
    fh=figure('Position',[10,10,400,300]);  %size of the gif!
    for i = 1:size(tri,1)
        patch(nx(tri(i,:),1),nx(tri(i,:),2),nx(tri(i,:),3),i,'FaceAlpha',0.85);
    end
    title('1st Brillouin zone');
    axis equal tight off vis3d
    camtarget([0,0,0]);
    loop=36;
    azimuth=360/loop*[0:loop-1];
    view(azimuth(1),30);
    drawnow;
    pause(0.1);
    mov=getframe(gcf);
    
    gifname = 'animation.gif';
    I = frame2im(mov);
    [X, map] = rgb2ind(I, 128);
    imwrite(X, map, gifname, 'GIF', 'WriteMode', 'overwrite', 'DelayTime', 0, 'LoopCount', Inf);
    
    for k=2:loop
     view(azimuth(k),30);
      drawnow;
      pause(0.1);
      mov=getframe(gcf);
      I = frame2im(mov);
        [X, map] = rgb2ind(I, 128);
        imwrite(X, map, gifname, 'GIF', 'WriteMode', 'append', 'DelayTime', 0);
    end
end
