First, let's recall the Newtonian theory of gravity: the gravitational force exerted by a point-like object on another point-like object is

where m1 and m2 are the two object masses, G is a cosmological constant, and d is the distance between the two objects. (By point-like, we mean that the two are far enough apart that ||d|| is much bigger than the object sizes.)

To make simulation simpler, it would be nice to separate the x, y, and z dimensions. After all, the force F can be described as the combination of its x, y, and z components.

And the distance d between an object pair can be written as

We also know that the force F points in the same direction as the vector d, which means that

therefore

(The same holds true for F_{y}, F_{z}.)

What's more, the gravitation force adds up nicely when there are more objects present. Of course, in the real world, things are continuous, and solving this means going through double, even triple integrals. Here comes discretization to the rescue. Given a high enough resolution, if we model our space and disc as a bunch of voxels (that is, 3d pixels), then we can treat each voxel as a constant-mass-and-density object, and easily sum up their interactions. This gives a rather good approximation of the continuous field.

here, the force felt by our observer (o) is the sum of forces exerted by each voxel (v1, v2, ..., vn),
d^{o,v1} is the distance vector from the observer to the voxel v1,
m_{v1}, m_{v2}, etc... are the masses of each voxel.

given the constant-mass-and-density approximation, m_{v1}, m_{v2}, etc... are all equal to some mass m_{vox}
thus we can write F_{x} as

Perfect! This equation makes it easy to simulate the force felt by any observer, at any point:
calculate d_{x} d_{y}, d_{z} between every point in the simulation grid and every occupied voxels.

dx = ox[None,None,None,:] - xx[:,:,:,None] dy = oy[None,None,None,:] - yy[:,:,:,None] dz = oz[None,None,None,:] - zz[:,:,:,None]

this gives you three big arrays, of shape (X, Y, Z, N_{occupied voxels})
We can combine them to get the norm of d at every point.

d = np.sqrt(dx * dx + dy * dy + dz * dz)

Now we can obtain the force at every point caused by every voxel by using the equation.

fx = dx / d * G / d**2 fy = dy / d * G / d**2 fz = dz / d * G / d**2

Finally, we sum over the occupied voxels to get the total force experienced at every point.

fx = np.sum(fx, axis=-1) fy = np.sum(fy, axis=-1) fz = np.sum(fz, axis=-1)

**Note**
given these simplifications, we can already notice an interesting peculiarity about the angle of the gravity field
around some massive object.
For example, let's fix the y coordinate, taking a cross-cut in the middle of our object,
which implies F_{y} = 0.
Now say we care about the angle theta of gravity at an arbitrary point on this cross-cut.
According to trigonometry,

which is

We see that the constant terms outside the sum cancel out! This leads us to our first discovery:
Given constant density throughout the disc, the angle of gravity at some latitude on our disc
depends **only on the geometric shape**, not the scale (bigger disc, same result), nor the mass (denser soil everywhere, same result),
nor the constant of gravity G.
Pretty neat!