% 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);