Hier werden die Unterschiede zwischen zwei Versionen gezeigt.
Nächste Überarbeitung | Vorhergehende Überarbeitung | ||
ws1314:stand_der_simulation_30.1.2014 [2014/02/05 22:23] stefanborn angelegt |
ws1314:stand_der_simulation_30.1.2014 [2016/05/10 14:46] (aktuell) |
||
---|---|---|---|
Zeile 1: | Zeile 1: | ||
- | ====== Stand der Simulation 30.1. ====== | + | ====== 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. | + | 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. | ||
+ | |||
+ | <code python> | ||
+ | # -*- 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 | ||
+ | |||
+ | |||
+ | </code> |