% Wigner–Seitz animation for a 2D parallelogram lattice
% Author: you + MATLAB
% Requires: base MATLAB (uses voronoin)

clear; close all; clc;

% ---------- Parameters ----------
a = 1.0;                 % |a1|
b = 1.0;                 % |a2|
theta_deg = [linspace(30,150,481), linspace(150,30,481)]; % sweep ∠(a1,a2)
Ngrid = 4;               % lattice extent: m,n in [-Ngrid..Ngrid]
draw_bisectors = true;   % show perpendicular bisectors to nearest neighbors
show_vectors   = true;   % show lattice vectors a1,a2
save_mp4       = true;   % write an mp4 file 
filename       = 'wignerSeitz.mp4';
vid = VideoWriter(filename,'MPEG-4');
vid.Quality = 100;
vid.FrameRate = 30;
open(vid);


% Axis limits that comfortably contain the Wigner–Seitz cell & neighbors
ax_lim = 2.2;

% ---------- Figure setup ----------
fh = figure('Color','w','Position',[100 100 800 600],'ToolBar','none');
axis equal; axis([-ax_lim ax_lim -ax_lim ax_lim]);
hold on; box on; set(gca,'FontSize',14);
set(gca,'linewidth',2.5);
%title('Wigner–Seitz Cell (Voronoi cell of the origin)');

% ---------- Helper to draw a long line segment within bounds ----------
clipline = @(p, n, L) [p - L*n; p + L*n]; % two endpoints along normal n through point p

% ---------- Main animation loop ----------
for k = 1:numel(theta_deg)
    cla; axis equal; axis([-ax_lim ax_lim -ax_lim ax_lim]); hold on; box on;

    th = deg2rad(theta_deg(k));
    a1 = [a, 0];                            % fix a1 on x-axis
    a2 = [b*cos(th), b*sin(th)];            % rotate a2 by theta

    % Build lattice points within a rectangle of indices
    [M,N] = meshgrid(-Ngrid:Ngrid, -Ngrid:Ngrid);
    P = M(:).*a1 + N(:).*a2;                % ( (2N+1)^2 x 2 ) array

    % Ensure origin exists and find its index
    [~,i0] = min(vecnorm(P,2,2));           % closest to (0,0) should be exactly the origin

    % Voronoi (Wigner–Seitz) of the lattice
    [V,C] = voronoin(P);                    % V: vertices, C{i}: polygon indices for site i
    poly  = V(C{i0},:);                     % polygon for the origin

    % Fill cell
    patch('Faces',1:size(poly,1),'Vertices',poly,...
          'FaceColor',[0.4 0.44 0.73],'FaceAlpha',0.25,'EdgeColor',[0.4 0.44 0.73],...
          'LineWidth',2);

    % Draw lattice points
    plot(P(:,1), P(:,2), 'k.', 'MarkerSize', 11);
    plot(0,0,'ro','MarkerFaceColor','r','MarkerSize',7);   % origin

    % Optionally draw lattice vectors
    if show_vectors
        quiver(0,0,a1(1),a1(2),0,'LineWidth',2,'MaxHeadSize',0.2,'Color',[0.2 0.2 0.2]);
        quiver(0,0,a2(1),a2(2),0,'LineWidth',2,'MaxHeadSize',0.2,'Color',[0.2 0.2 0.2]);
        text(a1(1)*1.07, a1(2)*1.07, 'a_1', 'Color',[0.2 0.2 0.2], 'FontSize',11);
        text(a2(1)*1.082, a2(2)*1.082, 'a_2', 'Color',[0.2 0.2 0.2], 'FontSize',11);
    end

    % Optionally draw perpendicular bisectors to the nearest neighbors
    if draw_bisectors
        % Find a handful of nearest distinct neighbors (exclude the origin)
        d = vecnorm(P - 0, 2, 2);
        d(i0) = inf;
        [~,idx] = sort(d,'ascend');
        K = 6; idx = idx(1:K);              % take K nearest
        Lseg = 4*ax_lim;                    % long segment length for visualization
        for j = 1:numel(idx)
            r = P(idx(j),:);                % neighbor vector
            mid = r/2;                      % midpoint
            n = [-r(2), r(1)];              % normal (perpendicular bisector direction)
            n = n / norm(n);
            seg = clipline(mid, n, Lseg);
            plot([seg(1,1) seg(2,1)], [seg(1,2) seg(2,2)], '-', 'Color', [0.6 0.6 0.6], 'LineWidth', 1);
            % also show line to neighbor (faint)
            plot([0 r(1)], [0 r(2)], ':', 'Color', [0.6 0.6 0.6]);
        end
    end

    % Cosmetic: outline cell edges thicker
    plot([poly(:,1); poly(1,1)], [poly(:,2); poly(1,2)], 'color', [0.4 0.44 0.73], 'LineWidth', 2);

    % Labels
    txt = sprintf('\\theta = %3.2f^\\circ', theta_deg(k));
    text(-ax_lim+0.125, ax_lim-0.17, txt, 'FontSize',14);

    drawnow;

    % Save GIF if requested
    frame = getframe(gcf);
    writeVideo(vid,frame);
end
% Output the movie as an mp4 file
close(vid);
