Um die Grundlage für weitere Arbeiten und Ideen zu legen, habe ich unter Verwendung eures Programmcodes das Programm umstrukturiert. Ich habe einige Kommentare hinzugefügt, glaube aber, dass es sinnvoll ist, wenn ihr selbst den Code analysiert und (in der Wiki) kommentiert.

Wichtige Neuerungen:

• Der Integrationsschritt ist jetzt eine Funktion. Diese kann man in der Zukunft um andere Integrationsmethoden erweitern. Bisher wird ein Euler-Schema verwendet.
• Alle erzeugten Planeten werden einer Liste planets hinzugefügt, so dass man die Berechnung der Bewegung durch einfache Schleifen aufschreiben kann. Um die Allgemeinheit des Ansatzes zu zeigen, werden zufällig Kleinplaneten erzeugt.
• Es gibt eine Benutzerinteraktion, die sich beliebig erweitern lässt: Wird die Taste 'd' gedrückt, so werden die Spuren der Planetenbahnen gelöscht. Mit der Taste 'a' wird das Autoscaling ein- und ausgeschaltet.
• Als Längeneinheit wird 1 Astronomische Einheit AE verwendet (der mittlere Abstand Erde-Sonne, hier wird nur ein ungefährer Wert verwendet, ihr könnt aber den genauen nachschlagen und einsetzen.) Das Newton'sche Kraftgesetz bleibt gültig, wenn man auch die Gravitationskonstante auf diese Längeneinheit umrechnet. Für die Masse wird weiter die Einheit kg, für die Zeit s verwendet. Als Beispiele wurden die echten (d.h. im Rahmen der gegenwärtigen Messgenauigkeit bestimmten) Massen von Erde, Sonne und Mars verwendet, sowie die zugehörigen Radien für Kreisbahnen nach dem 3. Kepler'schen Gesetz: $T^2=\frac{4\pi^2}{G(m+M)}$. In Wahrheit sind die Bahnen aber etwas exzentrisch (Ellipsen). Indem man die jeweils für eine Kreisbahn berechneten Geschwindigkeitsvektoren mit einer Zahl $\not=1$ multipliziert, erhält man exzentrische Bahnen.

# -*- 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