How to add Objects to a scene AND to a list? (Sphere Packing)
Creating an empty list, and adding objects to it

# Create empty list
lobjects = []
# adding an object to it
lobjects.append(rock_object)


If you create more objects, you keep appending them to lobjects.

Creating objects

If you already have a template rock object that you want to use to create all of the other rocks:


Duplicating the mesh data of the template rock object

original_data = rock_object.data
new_rock_obj = bpy.data.objects.new('CopiedMeshRock', original_data.copy())
Linking the original object's mesh data

new_rock_obj = bpy.data.objects.new('LinkedMeshRock', original_data)


Distribution of the rocks:

I think you could start by creating a grid, so an object is placed within each cell of this grid (I am using spheres for this example)



Then if you distort this grid, using any (radial?) distortion function, you can change the size of each cell of the grid, and fit a sphere at the centers with a radius less or equal to the smallest side of the cell.



You would only do this for cells that fall within the desired target pattern (a circle in your case).

I have used the following code to generate those images. If you set $a_x = a_y = 0.0$ you get the first image (no distortion); if you set $a_x$ and $a_y$ to positive values you get barrell distortion.

Blender

import bpy          # bpy.data.filepath (path to the current blend file)
import mathutils    # mathutils.Vector

class ExitError(Exception):
    pass

def addSphere(center, radius):
    '''
    Add a sphere at the given location, with the given radius

    @param center (mathutils.Vector) 3D location of the center of the sphere
    @param radius   (float) Radius of the sphere
    #return (bpy.types.object) newly created object
    '''
    bpy.ops.mesh.primitive_uv_sphere_add(radius = radius, location = center)

    return bpy.context.active_object

def distortPoint(p, ax, ay):
    '''
    Distort point p

    @param p    (mathutils.Vector) 2D point to be distorted
    @param ax   (float) Distortion coefficient in X
    @param ay   (float) Distortion coefficient in Y
    @return (mathutils.Vector) Distorted 2D point
    '''
    px = p.x * ( 1.0 - ax * p.length * p.length)
    py = p.y * ( 1.0 - ay * p.length * p.length)

    return mathutils.Vector((px, py))

def isPointWithinCircle(p, center, radius):
    '''
    Is point p within circle

    @param p        (mathutils.Vector) 2D point
    @param center   (mathutils.Vector) Center of circle
    @param radius   (float) radius
    @return         (boolean)
    '''
    return (p - center).length &lt; radius


def main():
    #fireworks.testDrawCubes(5, mathutils.Vector((2, 1, 0)))
    #fireworklib.drawCubes(2)
    scene = bpy.data.scenes.get('Scene')
    if not scene:
        raise ExitError('Could not retrieve scene')

    # Dimensions of the grid
    nrows = 50
    ncols = 50
    # Original size of each cell of the grid
    width = 0.8
    height = 0.8
    center = mathutils.Vector(( (nrows / 2.0) * width, (ncols / 2.0) * height ))
    # Normalized grid cell size
    nwidth = 1.0 / nrows
    nheight = 1.0 / ncols
    nwidth_half = nwidth / 2.0
    nheight_half = nheight / 2.0
    # Center of the normalized grid
    ncenter = mathutils.Vector(( (nrows / 2.0) * nwidth, (ncols / 2.0) * nheight ))
    # distortion coefficients
    ax = 0.9
    ay = 1.5
    # Circle: Center at center, and with radius
    radius = width * nrows / 3.0

    # List of created objects
    lobjects = []

    for i in range(nrows):
        for j in range(ncols):
            # Calculate center of the cell
            pc = mathutils.Vector(( i * nwidth + nwidth_half, j * nheight + nheight_half)) - ncenter
            # distort it
            dpc = distortPoint(pc, ax, ay) + ncenter
            # offset by center
            unnormalized_center = mathutils.Vector((dpc.x * width / nwidth, dpc.y * height / nheight))

            # only proceed if the center of the distorted cell falls within the circle
            if isPointWithinCircle(unnormalized_center, center, radius):
                p0 = mathutils.Vector((    i * nwidth,     j * nheight)) - ncenter
                p1 = mathutils.Vector((    i * nwidth, (j+1) * nheight)) - ncenter
                p2 = mathutils.Vector(((i+1) * nwidth,     j * nheight)) - ncenter
                p3 = mathutils.Vector(((i+1) * nwidth, (j+1) * nheight)) - ncenter

                # compute distorted points
                dpc = distortPoint(pc, ax, ay) + ncenter
                d0 = distortPoint(p0, ax, ay)
                d1 = distortPoint(p1, ax, ay)
                d2 = distortPoint(p2, ax, ay)
                d3 = distortPoint(p3, ax, ay)

                # approximate distorted bounding box
                dwidth = (d1 - d0).length * width / nwidth
                dheight = (d2 - d0).length * height / nheight
                # smallest side of the bounding box approximation
                dmin = min(dwidth, dheight) / 2.0

                # sphere's radio cannot be larger than smallest side of the approximate bounding box
                r = dmin

                lobjects.append(addSphere(mathutils.Vector((dpc.x * width / nwidth, dpc.y * height / nheight, 0.0)), r))

if __name__ == '__main__':
    try:
        main()
    except ExitError as e:
        print('Failed: {}'.format(e))


Note:* The example above creates a new sphere for each cell of the grid, which in this case is very inefficient: we could just reuse the same mesh and have the subsequent spheres link to it.

Another idea, if you use objects that are not spherical, is to:


give the objects (rocks) a small value to their location Z coordinate
enable rigid body dynamics
add a plane right below the given Z coordinate
Animate the scene so the objects fall by that small z, colliding with the plane and with each other, naturally randomizing them (without having to do that from code).


You could add rigid body dynamics for an object with

sphere = addSphere(mathutils.Vector((dpc.x * width / nwidth, dpc.y * height / nheight, 0.0)), r)
scene.rigidbody_world.collection.objects.link(sphere)                    

# add to the list of objects
lobjects.append(sphere)


Note: make sure that the target scene has rigid body world enabled:

if not scene.rigidbody_world:
    bpy.ops.rigidbody.world_add()

if not scene.rigidbody_world.collection:
    scene.rigidbody_world.collection = bpy.data.collections.new('RigidBodyWorld')