donut orbit

from __future__ import division
from visual import *
import random

loci = []           # where the positions will be stored

# All lengths in m, all masses in kg, all times in seconds

N = 40000            # number of mass points (will decrease later)
G = 6.7E-11         # universal gravitational constant
r = 1e6             # radius of "pastry" part of donut
R = 5e6             # distance from center of donut to line running
                    # around the donut, through the centroid of the
                    # actual donut part (a drawing would help here)
C = 2*pi*R          # circumference of donut through its pastry centroid
A = pi*r**2         # cross sectional area of donut
V = C*A             # volume of donut, using Pappus' theorem

startpos = vector(5e6,5e6,0)        # starting position of asteroid
v = vector(-800,400,0)           # starting velocity of asteroid
mA = 100                        # mass of asteroid
p = mA*v                        # initial momentum of asteroid
asteroid = sphere(pos=startpos, radius=1e5, color=color.cyan,
                  make_trail=True)

dt = 100            # time increment

# Build the donut by cutting it out of a cube of random points.
for i in range(N):
    rx = random.uniform(-R-r, R+r)
    ry = random.uniform(-R-r, R+r)
    rz = random.uniform(-R-r, R+r)
    radial = sqrt( rx**2 + rz**2 )
    centroidDistance2 = (R-radial)**2 + ry**2
    if centroidDistance2 < r**2:
        loci.append(vector(rx,ry,rz))

N = len(loci)
donut = points(pos=loci, size=3, color=color.green)

m = 5500*V/N        # mass of each point particle
                    # assuming density of donut is the same as that of Earth

while True:
    rate(1000)
    F = vector(0,0,0)
    for i in range(N):
        relpos = asteroid.pos - loci[i]
        F = F - G*m*mA/relpos.mag**2 * relpos.norm()
    p = p + F*dt
    v = p/mA
    asteroid.pos = asteroid.pos + v*dt