Planet shape in L1 gravity, 2D case

l1_gravity_planet.py

import numpy as np

N = 300 #size of the grid
M = 20000 #number of pixels

potential = np.zeros((N, N), dtype=np.float64)
masses = np.zeros((N, N), dtype=np.bool8)
def potential_of_mass_at(i, j):
    """Returns L1 potential of a mass at (i, j) with mass m
    potential is calcualted as:
    V(x, y) = -m/(|x-i|+|y-j|)
    """
    dx = np.abs(np.arange(N, dtype=np.float64) - j)
    dy = np.abs(np.arange(N, dtype=np.float64) - i)
    dxx, dyy = np.meshgrid(dx, dy)
    d = dxx+dyy #L1 distance
    #d = np.sqrt(dxx**2 + dyy**2) #L2 distance
    d[i,j] = 1e-5
    return -1.0/d

def add_mass_at(i, j):
    global potential, masses
    masses[i, j] = True
    potential += potential_of_mass_at(i, j)
    #at point i, j, set it to +infinity, so that it is not selected again
    potential[i, j] = np.inf

#initialize with one mass at the center
add_mass_at(int(N*0.5), int(N*0.5))
#now find the point with the lowest potential, not yet occupied by a mass
for idx in range(M-1):
    i, j = np.unravel_index(np.argmin(potential), potential.shape)
    add_mass_at(i, j)

import matplotlib.pyplot as plt
plt.contour(1.0/potential, levels=50, linewidths=0.5, colors='k')
plt.show()