function [V,dmin] = mmdlp(P,step); % function [V,dmin] = mmdlp(P,step); %------------------------------------------------------------- % Maximum Minimal Distance Lattice Partitioning % Copyright (c) 2001 by Ivan V. Bajic % This software may be freely used for research purposes, % as long as this notice stays attached. %------------------------------------------------------------- % This program finds the integer matrix V with |detV|=P, which % is a periodicity matrix for a sublettice of Z2 with maximum % minimal distance between points. This is achieved by first % finding periodicity matrix U for an optimal hexagonal lattice, % and then finding the best admissible V in the vicinity of U. % % Inputs: % P = number of partitions % step = step size (in radians) for rotating the optimal hexagonal % lattice; usually 0.01 <= step <= 0.0001 works well % % Outputs: % V = periodicity matrix for optimal sublattices of Z2 % dmin = minimal distance between the points of the lattice % defined by V %-------------------------------------------------------------- % find initial hex matrix x = sqrt(P/(2*sqrt(3))); U0 = [2*x x; 0 x*sqrt(3)]; dmin0 = 0; % (initial) minimal distance theta0 = 0; % loop for theta = 0:step:pi R = [cos(theta) sin(theta); -sin(theta) cos(theta)]; U = R*U0; % To find "nearest" integer matrix, try 16 cases: % for each coordinate, try floor and ceil. % Combination c is number between 0 and 15, binary % representation of c describes which combination % (floor/ceiling) is used for which coordinate. V = zeros(2,2); for c = 0:15 % first bit if floor(c/8) V(1,1) = ceil(U(1,1)); else V(1,1) = floor(U(1,1)); end % second bit b1 = mod(c,8); if floor(b1/4) V(2,1) = ceil(U(2,1)); else V(2,1) = floor(U(2,1)); end % third bit b2 = mod(b1,4); if floor(b2/2) V(1,2) = ceil(U(1,2)); else V(1,2) = floor(U(1,2)); end % fourth bit if mod(b2,2) V(2,2) = ceil(U(2,2)); else V(2,2) = floor(U(2,2)); end v1 = V(:,1); v2 = V(:,2); dmin = min([norm(v1) norm(v2) norm(v1-v2) norm(v1+v2)]); if (dmin > dmin0 & abs(det(V))==P) % best solution so far dmin0 = dmin; theta0 = theta; c0 = c; end end end % at this point, theta0 is the optimal angle R = [cos(theta0) sin(theta0); -sin(theta0) cos(theta0)]; U = R*U0; % rounding according to c0 % first bit if floor(c0/8) V(1,1) = ceil(U(1,1)); else V(1,1) = floor(U(1,1)); end % second bit b1 = mod(c0,8); if floor(b1/4) V(2,1) = ceil(U(2,1)); else V(2,1) = floor(U(2,1)); end % third bit b2 = mod(b1,4); if floor(b2/2) V(1,2) = ceil(U(1,2)); else V(1,2) = floor(U(1,2)); end % fourth bit if mod(b2,2) V(2,2) = ceil(U(2,2)); else V(2,2) = floor(U(2,2)); end v1 = V(:,1); v2 = V(:,2); dmin = min([norm(v1) norm(v2) norm(v1-v2) norm(v1+v2)]);