Um die Grundlage für weitere Arbeiten und Ideen zu legen, habe ich unter Verwendung eures Programmcodes das Programm umstrukturiert.

# -*- coding: utf-8 -*-
from __future__ import division
from visual import *
import numpy as np
 
 
 
#  Einheit der Länge 1.5 * 10**11 Meter (eine 'Astronomische Einheit'
#  Entsprechend  Graviationskonstante   6.672 10**(-11) / (1.5*10**11)**3
 
G= 6.672 *10**(-11) / (1.5*10**11)**3
 
#  Alle anderen Einheiten sind SI-Einheiten (kg,s)
 
 
def beschleunigung(p1,p2):
	'''gibt die Beschleunigung des Körpers p1 durch die von p2 wirkende
	Kraft zurück.'''
	if (np.linalg.norm(p1.pos-p2.pos)>0.000000001):
		beschleunigungsvektor=G*p2.mass*(p2.pos-p1.pos)/(np.linalg.norm (p2.pos-p1.pos)**3)
	else: # ACHTUNG: Das ist kein 'elastischer Stoß', sondern etwas eher Merkwürdiges
		beschleunigungsvektor=p1.pos-p2.pos
	return beschleunigungsvektor
 
 
def schritt(dt, mode): 
	'''Berechnet einen Zeitschritt des Systems. 
	mode bezeichnet die Integrationsmethode. Bisher
	nur implementiert:   mode='Euler' '''
	if mode=='Euler':
		for planet in planets:
			planet.pos = planet.pos +  planet.velocity*dt		
			planet.trail.append(pos=planet.pos)
			beschl=vector(0,0,0)
			for planetplanet in planets:
				if not planet is planetplanet:
					beschl+=beschleunigung(planet,planetplanet)
			#print beschl
			planet.velocity = planet.velocity + beschl*dt
	elif mode=='RK':  #Runge-Kutta
		pass
	else:
		print 'Unzulässiger Modus'
		raise Exception("Unzulässiger Integrationsmodus")
 
 
 
##### Hauptprogramm ##########################################
 
 
scene.autoscale=True
 
# Initialisieren der Planeten, alle werden eine Liste 'planets' zugefügt
 
planets=[]
 
# Stelle Planeten 500* größer dar als sie sind, Sonne 20*
 
 
 
# Erde
planet1 = sphere (pos=(1,0, 0), radius=500*6.4*10**6/(1.5*10**11), color=color.blue)
planet1.velocity = vector(0.0,1.*2*np.pi/(3.15*10**7),0.0)
planet1.trail=curve(color =color.blue)
planet1.mass=5.97*10**24
 
planets.append(planet1)
 
#Sonne
planet2 = sphere (pos=(0,0,0), radius=20*7*10**8/(1.5*10**11), color=color.orange)
planet2.velocity = vector(0,0,0)
planet2.trail=curve(color =color.orange)
planet2.mass=1.99*10**30
planets.append(planet2)
 
#Mars
planet3= sphere (pos=(0,1.524,0), radius=500*3.85*10**6/(1.5*10**11), color=color.red)
planet3.velocity = vector(-1.524*2*np.pi/(687.*86400),0.0,0.0)
planet3.trail=curve(color =color.red)
planet3.mass=6.4*10**23
planets.append(planet3)
 
# Kleinstplaneten, Größe auf 2/5 der Erddarstellung gesetzt
 
for i in range(20):
	alpha=np.random.rand()*2*np.pi
	r=(np.random.rand()*2.+0.3)
	umlaufzeit=(4*np.pi**2/(G*planet2.mass) *r**3)**0.5   # Das wäre die Umlaufzeit bei Kreisbahn nach den Keplerschen Gesetzen
 
	planet=sphere(pos=(r*np.cos(alpha),0.+r*np.sin(alpha),0),radius=200*6.4*10**6/(1.5*10**11),color=color.yellow)
	planet.velocity=vector(-r*np.sin(alpha), r*np.cos(alpha),0)*2*np.pi/umlaufzeit
	# Durch Multiplizieren der velocity mit  *(1+0.2*np.random.rand()) wird eine Zufällige Exzentrizität bewirkt.
	planet.mass=np.random.rand()*10**10+10**9
	planet.trail=curve(color=color.yellow)
	print "umlaufzeit ",i, "r", r,"Um",  umlaufzeit, "alpha", alpha, "v ", planet.velocity
	planets.append(planet)
 
#time.sleep(59)
 
######################################################################################
 
# Eigentliche Simulation
 
 
dt = 43300.    # halber Tag in Sekunden
TRAILLENGTH=1   # Drückt man die Taste 'd', wird alles bis auf die letzten TRAILLENGTH Stücke der Trails gelöscht
rate(10)
n=0
 
while 1:
	n=n+1
 
	# Berechnung:
 
	schritt(dt,'Euler')
 
	# Interaktion mit Benutzer:
 
	if scene.kb.keys: # event waiting to be processed?
		s = scene.kb.getkey() # hole ein Zeichen aus der Warteschlange
		if s=='d':
			for planet in planets:
				planet.trail.x=planet.trail.x[-TRAILLENGTH:]
				planet.trail.y=planet.trail.y[-TRAILLENGTH:]
				planet.trail.z=planet.trail.z[-TRAILLENGTH:]
		if s=='a':
			scene.autoscale=not scene.autoscale