Matlab: create close packed layer of circles

Question

I'm trying to create a close packed layer by generating my center points like a box and then shear the point in order to create a close packed layer of circles.. more info can be found here

but I'm having some difficulties.. My code so far is

%% Trying out the shear function
rad=2; n=3;
[X,Y] = meshgrid(0   :   rad*2   :   rad*(n-1)*2    ,   ...
             0   : sqrt(2*(2*rad)^2)/2 : sqrt(2*(2*rad)^2)/2*(n-1));
xyBox = [reshape(X,1,numel(X)) ; reshape(Y,1,numel(Y))];
Sh = @(m) [1 m; 0 1];       % horizontal shear - slope is qual to 1/m
slope = sqrt(3);

shearedCoordinates = Sh(1/slope) * xyBox;
figure; plot(xyBox(1,:),xyBox(2,:),'k.');
for i= 1:(numel(shearedCoordinates)/2)
    hold on;
    circle(shearedCoordinates(1,i),shearedCoordinates(2,i),rad)
    plot(shearedCoordinates(1,i),shearedCoordinates(2,i),'bx')
    hold off;
end
axis equal

Don't really understand the math behind the Sh, but can just see that it does a damn good job twisting (shearing) points. I calculated a slope of square root 3 should give the right place, but it looks lige there is something wrong with my y distances...

the circle function is

function circle(x,y,r)
    angle=0:0.01:2*pi; 
    xp=r*cos(angle);
    yp=r*sin(angle);
    plot(x+xp,y+yp,'r');
end

* SOLUTION *

%% Trying out the shear function
rad=2; n=3;
[X,Y] = meshgrid(0   :   rad*2   :   rad*(n-1)*2    , ...
    0   : sqrt(3)*rad   : sqrt(3)*rad*(n-1));
xyBox = [reshape(X,1,numel(X)) ; reshape(Y,1,numel(Y))];
Sh = @(m) [1 m; 0 1];       % horizontal shear - slope is qual to 1/m
slope = sqrt(3);

shearedCoordinates = Sh(1/slope) * xyBox;
figure; plot(xyBox(1,:),xyBox(2,:),'k.');
for i= 1:(numel(shearedCoordinates)/2)
    hold on;
    circle(shearedCoordinates(1,i),shearedCoordinates(2,i),rad)
    plot(shearedCoordinates(1,i),shearedCoordinates(2,i),'bx')
    hold off;
end
axis equal

(don't forget to vote up if you found it useful)

Answer 1

By my reckoning your mistake (or perhaps just the first one that I can see) is in the y-spacing of the lines of circle centres. You have the expression

0   : sqrt(2*(2*rad)^2)/2 : sqrt(2*(2*rad)^2)/2*(n-1)

for laying out the circle centres. I think that the vertical stride (or step value) should be rad*sqrt(3)/2 which would give you:

0   : rad*sqrt(3)/2 : sqrt(2*(2*rad)^2)/2*(n-1)

Unless I've made a mistake (entirely possible) the height of a unit equilateral triangle is sqrt(3)/2 and that ought to be your y-spacing.