This section is devoted to finding the gravitational potential field of odd shaped constant density 3D things to support this higher goal. First and second moments are also required. The general plan is to express those shapes as the disjoint union of simple shapes whose potential field is given by some formula. The field of the union is the sum of the fields of the simple shapes. This paper which gives a somewhat complex formula for the cube’s field, stimulated this quest.

Here are some schemes:

Monte Carlo versions of some of these should be considered!

Spherical Harmonics may facilitate some of these schemes.

As we move mass on the surface to make the potential at the surface more uniform, the vector angular momentum will change. We compute the new covariance tensor and modify the angular velocity to preserve the old vector angular momentum. This preserves the invariant that the spin axis remains the principle axis of the tensor inertia; the axis of rotation is aligned with the angular momentum.
I can use Jacoby’s routine to compute an eigenvector for the inertia tensor of the blob, and spin about that. The covariance tensor has the same eigenvectors as the inertia tensor.

notes on that code, Algebra with computer help, Initial configuration
earlier notes, These pages note that the spheroid is too simple.
Solution of new spin axis: Compute biggest eigenvector after moving mass.
MP clues