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