How to create a circle filled with equidistant points inside it?

I want to create circle filled with equidistant points inside it. I have tried making the following program but it's taking too long to give output

radius = 100;
xc = 1;
yc = xc+1;
jc= 1;
jy=jc+1;
area = pi*radius*radius;
for i=xc:yc:area
    theta = i*(2*pi);
    r = sqrt(area)*radius;
    for j=jc:jy:i
        x = xc + r.*cos(theta);
        y = yc + r.*sin(theta);
    end
    plot(x,y,'.')
end

[Stephen23]

@Rushil Asthana: what is your definition of "circle filled with equidistant points inside it" ?

[Rushil Asthana]

Points that are at equal distance from each other. These points should be inside the circle.

[John D'Errico]

There is no simple scheme to optimally solve this problem for a fully general domain. Even a simple domain can be difficult. You can think of points that are equidistant as represented by the centers of circles of a fixd radius. Now you should see this is a case of the well known circle packing problem.

https://en.wikipedia.org/wiki/Circle_packing

[Stephen23]

@Rushil Asthana: because you did not mention anything about an optimal solution or any specific boundary requirements then this task is trivially easy: first generate a square grid of points large enough to cover the required circle, then simply discard those outside the circle. Reasonably efficient and should take only a few lines of code.

That will provide you with exactly what you asked for: points with exactly equal spacing all within a circle.

[Rushil Asthana]

@John D'Errico What will the code look like essentially? Cause I need help with exactly what you mentioned.

Answer 1

[Chunru]

Here is a sub-optimal (or almost optimal) solution. Filled dots will form the hexagonal formation.

radius = 10;
d = 1;  % distance
xall=[]; yall=[];
dy = sqrt(3)/2 * d;
ny = floor(radius/dy);

for i=-ny:ny
    y = dy*i;
    if rem(i, 2)==0
        nx = floor((sqrt(radius^2 - y^2))/d);
        x = (-nx:nx)'*d;
    else
        nx = floor((sqrt(radius^2 - y^2)-d/2)/d);
        x = (-nx-0.5:nx+0.5)'*d;
    end
    xall = [xall; x];
    yall = [yall; y*ones(size(x))];
end

plot(xall(:), yall(:), '.');
hold on
theta = 0:360;
plot(radius*cosd(theta), radius*sind(theta), 'r')
axis equal

[Rushil Asthana]

Thanks man! This is exactly what i needed.

[Image Analyst]

@Rushil Asthana, then could you please click the "Accept this answer" link to grant @Chunru the reputation points he deserves for helping you. Thanks in advance. 🙂

[Chunru]

@Image Analyst Thanks for cultivating the good community culture. :-)