Code that outputs the center coordinates of a circle in the method of 2D packing the same circle of a regular hexagon

Question

I'm writing a matlab code that outputs the center coordinates of the circle in how to 2D pack the same circle of regular hexagon. It's a picture of the center coordinates of the circle that I want to output:

The number of circles increases by 6*n each time a layer is added from the center of the circle. It has 1 circle in layer 1, 6 circles in layer 2, and 12 circles in layer 3.

I want to make this code simple using the Polar coordinate system.

The code I wrote is a Cartesian coordinate system that outputs the coordinates of all the corresponding circles when there are n layers.

Diameter=2;
r=Diameter/2;

x_values = {};
y_values = {};

layer = 5;
Tcircle = layer * 6;
n = layer;

a = 1;
number=1;
while a <= n

    temp_x_values = [];
    temp_y_values = [];
    
    x = (-r) * a;
    y = a * 2*r;
    z=x;

    temp_x_values = [temp_x_values, x];
    temp_y_values = [temp_y_values, y];

    temp_x_values = [temp_x_values, x];
    temp_y_values = [temp_y_values, -y];

    b = 0;
    while b < a  % 
        x = x + 2*r;

        temp_x_values = [temp_x_values, x];
        temp_y_values = [temp_y_values, y];


        temp_x_values = [temp_x_values, x];
        temp_y_values = [temp_y_values, -y];
        
        b = b + 1;
    end

    diagonal_x = x;
    diagonal_y = y; % 

    i = 0;
    diagonal_x2 = z ;
    while i < a-1  % 
        diagonal_x = diagonal_x + 1*r;
        diagonal_y = diagonal_y - 2*r;
        diagonal_x2 = diagonal_x2 - 1*r;
        diagonal_y2 = diagonal_y;     

        temp_x_values = [temp_x_values, diagonal_x];
        temp_y_values = [temp_y_values, diagonal_y];
        
        temp_x_values = [temp_x_values, diagonal_x2];
        temp_y_values = [temp_y_values, diagonal_y2];

        temp_x_values = [temp_x_values, diagonal_x];
        temp_y_values = [temp_y_values, -diagonal_y];
        
        temp_x_values = [temp_x_values, diagonal_x2];
        temp_y_values = [temp_y_values, -diagonal_y2];
        
        i = i + 1;
    end
    
    zero_x= a*2*r;
    zero_y= 0;
    temp_x_values = [temp_x_values, zero_x];
    temp_y_values = [temp_y_values, zero_y];

    temp_x_values = [temp_x_values, -zero_x];
    temp_y_values = [temp_y_values, zero_y];

    x_values{a} = temp_x_values;
    y_values{a} = temp_y_values;
    
    fprintf('\n%d layer x : ', a);
    fprintf('%s\n', mat2str(x_values{a}));
    fprintf('%d layer y : ', a);
    fprintf('%s\n', mat2str(y_values{a}));      

    for aa=1:a  
        aa;
        B=x_values{aa};
        C=y_values{aa};
        b=size(B);
        for bb=1:b(1,2)
            x_location(number)=B(bb);
            y_location(number)=C(bb);
            number=number+1;  
        end
    end

    a = a + 1;
end
x_location(number+1)=0;
y_location(number+1)=0;
scatter(x_location, y_location)

This is the result.

Layer 1:

x : [-1 -1 1  1 2 -2]
y : [ 2 -2 2 -2 0  0]

Layer 2:

x : [-2 -2 0  0 2  2 3 -3  3 -3 4 -4]
y : [ 4 -4 4 -4 4 -4 2  2 -2 -2 0  0]

Layer 3:

x : [-3 -3 -1 -1 1  1 3  3 4 -4  4 -4 5 -5  5 -5 6 -6]
y : [ 6 -6  6 -6 6 -6 6 -6 4  4 -4 -4 2  2 -2 -2 0  0]

Layer 4:

x : [-4 -4 -2 -2 0  0 2  2 4  4 5 -5  5 -5 6 -6  6 -6 7 -7  7 -7 8 -8]
y : [ 8 -8  8 -8 8 -8 8 -8 8 -8 6  6 -6 -6 4  4 -4 -4 2  2 -2 -2 0  0]

Layer 5:

x : [-5  -5 -3  -3 -1  -1  1   1  3   3  5   5 6 -6  6 -6 7 -7  7 -7 8 -8  8 -8 9 -9  9 -9 10 -10]
y : [10 -10 10 -10 10 -10 10 -10 10 -10 10 -10 8  8 -8 -8 6  6 -6 -6 4  4 -4 -4 2  2 -2 -2  0   0]

Answer 1

First you should notice that your current code, as it is, does not "pack" the circles tightly on the Y axis. If I run your code (with only 2 layers for easier visualisation), and fully draw the circles you'll notice the gap between the Y slices:

---

To simplify the calculation of all the centers, I did not use a polar coordinate approach, but I took advantage of the many symetries such an arrangment of circle will have:

• The middle slice is constructed easily: all circles stacked on a line y=0, we know the number of circles, and each center is separated from the next one by a distance equal to the diameter of a circle.
• The slice above: There is one less circle in total. Each x center is in the middle of the position of the centers in the line below. Another way to say that is that the x centers for slice 1 are the x center for slice 0 + an offset equal to the radius of one circle. Once again, all the y centers for this line are equal, with the value yc = Diameter * sqrt(3)/2
• The slice below: Nothing new to calculate here. xc are the same than for the slice above, and yc for this slice is the negative of the above slice (symmetry)
• Then for each additional slice above, just repeat the procedure, and mirror it for the corresponding slice under.

In code, it would look like this:

%%
D = 2 ;     % Circle diameter
R = D/2 ;   % Radius
N = 6 ;     % Number of layers

% Total number of circles
Ncircles = 1 + sum(6*(1:N)) ;

% Preallocate output vectors
xc = NaN(Ncircles,1) ;
yc = NaN(Ncircles,1) ;

% maximum |X| center value, and Y base step
xmax = N * D ;
baseYstep = sqrt(3)/2 ;

% X Centers for the middle line/slice
xslice = (-xmax:D:xmax) ;

% Save values in final arrays
xptr = numel(xslice) ;
xc(1:xptr) = xslice ;
yc(1:xptr) = zeros(1,xptr) ; % Y Centers for the middle line (All zeros obviously)

% now for each line above (and then mirrored below)
for slice =1:N
    % To get X center value, remove one element and offset all the
    % remaining element by the radius:
    xslice = xslice(2:end) - R ;
    % Y centers for this slice all have same value
    nps = numel(xslice) ;
    yslice = ones(1,nps) * (slice*D*baseYstep) ;
    
    % Save centers in output arrays (Once for the slice "above" and one for
    % the slice "below"
    for side = [1,-1]
        istart = xptr + 1 ;
        istop  = xptr + nps ;
        xc(istart:istop) = xslice ;
        yc(istart:istop) = yslice * side ;
        xptr = istop ;
    end
end

hold on, grid on, axis equal
hp = plot(xc,yc,'+r') ;

This will give you neatly stacked circles in an hexagonal pattern:

And it works for higher orders of N, here with N=6:

---

Code to add to the example for full circle drawing:

%% Full visualisation to make sure the circles are tangent
% Just to preallocate an array of graphic handles
hc(Ncircles,1) = plot(NaN,NaN) ; delete(hc(Ncircles)) ;

% Now plot all circles
tt = linspace(0,2*pi,36) ;
for ic=1:Ncircles
    xp = R * cos(tt) + xc(ic) ;
    yp = R * sin(tt) + yc(ic) ;
    hc(ic) = plot(xp,yp,'b') ;
end