====== Stand der Simulation 30.1.2014 ====== 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