Packing Dots

b = 0.3;
(*
	can take a while to calculate points[i] i=1:100*maxiter - could be sped up but doesn't really matter
*)
maxiter = 3;
Density[points_, x_, y_] :=
  Total[Map[Exp[-Norm[# - {x, y}]^2/b^2] &, points]] -
   3 Exp[-x^2 - y^2];
points[0] = RandomReal[NormalDistribution[0, 1], {100, 2}];
j = 0;
xs = Table[xt, {xt, -2, 2, 4/100}];
Do[
  Do[

    {x, y} = points[j][[i]];
    ps = Join[ points[j][[;; i - 1]], points[j][[i + 1 ;;]]];
    vals = Table[Density[ps, xt, y], {xt, xs}];

    {x, y} = {xs[[
        Position[vals, Min[vals]][[1]][[1]]
        ]], y};

    vals = Table[Density[ps, x, yt], {yt, xs}];

    {x, y} = {x, xs[[
        Position[vals, Min[vals]][[1]][[1]]
        ]]};

    j += 1;
    points[j] =
     Join[ points[j - 1][[;; i - 1]], {{x, y}},
      points[j - 1][[i + 1 ;;]]];
(* rotate to get rid of coordinate effects*)
    points[j] = Map[N[RotationMatrix[1].#] &, points[j]];
    , {i, Length[points[0]]}];
  , {iter, maxiter}];
Manipulate[
 Graphics[Point[ Map[N[RotationMatrix[-i].#] &, points[i]]],
  PlotRange -> 4], {i, 0, 100*maxiter, 1}]