import matplotlib as mpl from mpl_toolkits.mplot3d import Axes3D import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation import time # Einstellungen fuer die Grafik: fig = plt.figure() ax = fig.add_axes([0,0,1,0.9], projection='3d',label='axes3') # 3D Grafik ax_head = fig.add_axes([0,0.9,1,1],label='axes_head') # Kopfgrafik, in der Zeit dargestellt wird ax.set_aspect('equal') class Himmelskoerper(object): Liste_der_Kraefte = [] # die Liste der Kraefte m = 1.0 # Masse des Himmelskoerpers, bei der Rakete sogar variabel v = np.array([0.0,0.0,0.0]) # Geschwindigkeitsvektor ort = np.array([0.0,0.0,0.0]) # Ortsvektor name = "" # Name des Himmelskoerpers def __repr__(self): return self.name def __init__(self,name,m): self.name = name self.m = m # Ab hier beginnen die Methoden zur Ortsberechnung. Damit man nicht jede Methode einzeln aufrufen muss, # greift jede Methode auf die jeweils dafuer benoetigte Methode zu, sodass am Ende im Programm nur die letzte Methode benutzt werden muss def berechne_gravitationskraft(self,objekt1): """ Gibt die Gravitationskraft aus, die zwischen den Objekten "objekt1" und "objekt2" herrscht Das Ergebnis wird als Vektor ausgegeben, immer in Richtung von objekt1 bis objekt2 """ verbindungsvektor = self.ort-objekt1.ort betrag = np.linalg.norm(verbindungsvektor) F_betrag = (welt.gravitationskonstante * objekt1.m * self.m)/betrag**2 F = (F_betrag/betrag)*verbindungsvektor return F def berechne_kraefte_neu(self,welt): """ Aktualisiert die Werte in der Liste der Kraefte, indem alle Werte neu berechnet werden """ self.Liste_der_Kraefte = [] for i in welt.Liste_der_Objekte: if i.name != self.name: k = i.berechne_gravitationskraft(self) self.Liste_der_Kraefte.append(k) def berechne_resultierende_Kraft(self,welt): """ Aus der Liste der Kraefte wird eine resultierende Kraft in Form eines Vektors berechnet """ self.berechne_kraefte_neu(welt) r = np.array([0.0,0.0,0.0]) #resultierende Kraft for i in self.Liste_der_Kraefte: r = r + i return r def berechne_beschleunigung(self,welt): """ Aus der resultierende Kraft wird mit der Masse die Beschleunigung berechnet """ r = self.berechne_resultierende_Kraft(welt) a = r/self.m return a def berechne_geschwindigkeitsvektor(self,zeitschritt,welt): """ Aus der Beschleunigung wird mit dem uebergebenen Zeitschritt der Geschwindigkeitsvektor berechnet """ a = self.berechne_beschleunigung(welt) v = a*zeitschritt return self.v + v #alte und neue Geschwindigkeit addiert def berechne_ort(self,zeitschritt,welt): """ Nun wird wieder mithilfe des Zeitschrittes aus der Geschwindigkeit die zurueckgelegte Strecke berechnet und daraus folgt dann der neue Ort """ v = self.berechne_geschwindigkeitsvektor(zeitschritt,welt) self.v = v s = v*zeitschritt #zurueckgelegte Strecke im Zeitintervall ort = self.ort + s return ort def leapfrog(self,zeitschritt,welt): """ Die Methode berechne_ort ist eigentlich schon ausreichend, allerdings gibt es Probleme mit der Genauigkeit. Hier wird das geloest, indem das Leapfrog Verfahren angewendet wird. """ ort_zwischen = self.ort + self.v * (zeitschritt/2.0) self.ort = ort_zwischen a_zwischen = self.berechne_beschleunigung(welt) self.v = self.v + a_zwischen * zeitschritt self.ort = ort_zwischen + self.v * (zeitschritt/2.0) return self.ort #------------------------------------------------------------------------------------------------------------------------------------ class Rakete(Himmelskoerper): richtungsvektor = np.array([0.0,0.0,0.0]) # gibt Richtung an, in die die Spitze der Rakete zeigt, in diese Richtung wirkt auch der Schub fuellmenge = 0.8 # Gibt relative Treibstoffmenge in gefuellter Stufe an [gesamtmasse_der_stufe = leermasse*fuellmenge/(1-fuellmenge)] leermasse_stufen = [] # Die einzelnen Zahlen stehen in der Liste, angefangen bei der ersten Stufe brennrate = [] # Die einzelnen Zahlen stehen in der Liste, angefangen bei der ersten Stufe momentan_brennende_stufe = 0 # gibt an, welche Stufe gerade brennt, wenn Null brennt keine Stufe anzahl_stufen = len(leermasse_stufen) schubkraft = 0 # maximale Schubkraft der Rakete schub = 0 # momentaner Schub der Rakete aufenthaltsort = "space" # gibt den Planeten an, auf dem sich die Rakete gerade befindet. Befindet sie sich auf keinem Planeten, so steht dieser Wert auf "space" momentane_gesamtmasse_der_brennenden_stufe = 0.0 def __init__(self, richtungsvektor, m, fuellmenge, leermasse_stufen, brennrate, name, schubkraft): self.richtungsvektor = richtungsvektor self.m = m self.fuellmenge = fuellmenge self.leermasse_stufen = leermasse_stufen self.brennrate = brennrate self.name = name self.schubkraft = schubkraft def massenberechnung(self,zeitschritt): """ Berechnet den Verlust der Masse durch die Verbrennung des Treibstoffes. Dies ist natuerlich nur noetig wenn eine Stufe brennt """ if self.momentan_brennende_stufe != 0: self.m = self.m - (self.brennrate[self.momentan_brennende_stufe-1]*zeitschritt) def berechne_raketenbeschleunigung(self,welt,zeitschritt): """ Berechnet Beschleunigung der Rakete aus zwei Teilen: der Gravitationskraft der anderen Himmelskoerper und der Beschleunigung der Rakete selbst, inklusive Massenverlust """ r = self.berechne_resultierende_Kraft(welt) self.massenberechnung(zeitschritt) a_gravitation = r/self.m a = a_gravitation + (self.schub/self.m)*(self.richtungsvektor/np.linalg.norm(self.richtungsvektor)) return a # Die naechsten drei Methoden machen genau das gleiche wie die letzen drei Methoden der Himmelskoerper Klasse, # nur das jetzt die Raketenbeschleunigung benutzt wird def berechne_raketengeschwindigkeitsvektor(self,zeitschritt,welt): """ Aus der Beschleunigung wird mit dem uebergebenen Zeitschritt der Geschwindigkeitsvektor berechnet """ a = self.berechne_raketenbeschleunigung(welt,zeitschritt) v = a*zeitschritt return self.v + v #alte und neue Geschwindigkeit addiert def berechne_raketenort(self,zeitschritt,welt): """ Nun wird wieder mithilfe des Zeitschrittes aus der Geschwindigkeit die zurueckgelegte Strecke berechnet und daraus folgt dann der neue Ort """ v = self.berechne_raketengeschwindigkeitsvektor(zeitschritt,welt) self.v = v s = v*zeitschritt #zurueckgelegte Strecke im Zeitintervall ort = self.ort + s return ort def raketenleapfrog(self,zeitschritt,welt): """ Die Methode berechne_ort ist eigentlich schon ausreichend, allerdings gibt es Probleme mit der Genauigkeit. Hier wird das geloest, indem das Leapfrog Verfahren angewendet wird. """ ort_zwischen = self.berechne_raketenort(zeitschritt,welt) + self.v * (zeitschritt/2.0) self.ort = ort_zwischen a_zwischen = self.berechne_raketenbeschleunigung(welt,zeitschritt) self.v = self.v + a_zwischen * zeitschritt self.ort = ort_zwischen + self.v * (zeitschritt/2.0) return self.ort def stufe_starten(self,stufe): """ Eine einfache Methode, in der die Schubkraft aufgrund der brennenden Stufe angepasst wird """ self.momentan_brennende_stufe = stufe self.schub = self.schubkraft self.momentane_gesamtmasse_der_brennenden_stufe = self.leermasse_stufen[stufe-1]*self.fuellmenge/(1-self.fuellmenge) def stufe_abwerfen(self,stufe): """ Hier erfolgt der Massenverlust durch die abgeworfene Stufe """ self.m = self.m - self.leermasse_stufen[stufe-1] def treibstoffverbrauch(self,zeitschritt): """ Eine wichtige Methode, in der der Treibstoffverbrauch ermittelt wird. Faellt die Gesamtmasse der Stufe unter die Leermasse der Stufe, so folgt daraus dass der Treibstoff verbraucht ist und es keinen Schub mehr gibt """ gesamtmasse_der_stufe = self.leermasse_stufen[self.momentan_brennende_stufe-1]*self.fuellmenge/(1-self.fuellmenge) self.momentane_gesamtmasse_der_brennenden_stufe = self.momentane_gesamtmasse_der_brennenden_stufe - (self.brennrate[self.momentan_brennende_stufe-1]*zeitschritt) if self.momentane_gesamtmasse_der_brennenden_stufe <= self.leermasse_stufen[self.momentan_brennende_stufe-1]: self.momentan_brennende_stufe = 0 self.schub = 0 #------------------------------------------------------------------------------------------------------------------------------------ class Welt(object): gravitationskonstante = 6.673 * 10**(-11) #Gravitationskonstante Liste_der_Objekte = [] def objekt_hinzufuegen(self,objekt,ort): """ Fuegt der Welt ein Objekt hinzu """ self.Liste_der_Objekte.append(objekt) objekt.ort = ort def objekt_als_satellit_hinzufuegen(self,objekt,satellit,abstand): """ Fuegt ein Objekt ("satellit" genannt) in die Welt ein, welches sich um ein anderes Objekt ("objekt" genannt) dreht. Es ist dabei in einem stabilen Orbit. Der Abstand ("abstand" genannt) ist dabei der Abstand zwischen den Mittelpunkten. Ansatz: F_g = F_r Nach Aufloesung erhaelt man: v = sqrt(gravitationskonstante*masse_objekt/abstand) """ satellit.v = np.array([0.0,(self.gravitationskonstante*objekt.m/abstand)**0.5,0.0]) ort = objekt.ort + np.array([abstand,0,0]) self.objekt_hinzufuegen(satellit,ort) #------------------------------------------------------------------------------------------------------------------------------------ erde = Himmelskoerper("Erde",5.972E24) # die Erde wird mit korrekter Masse hinzugefuegt rakete = Rakete(np.array([1.0,0.0,0.0]),500000.0,0.8,[40000.0,40000.0],[1000.0,1000.0],"Rakete",2000000.0) # die Rakete erhaelt realistische Werte welt = Welt() # es wird eine Welt definiert, in die Objekte eingefuegt werden koennen welt.objekt_hinzufuegen(erde,np.array([0.0,0.0,0.0])) # die Erde wird in der Mitte eingefuegt welt.objekt_hinzufuegen(rakete,np.array([0.0,3e7,2e7])) # die Rakete wird etwas weiter weg von der Erde eingefuegt raketenbahn = [] # fuer die artist Animation wird eine leere Liste erzeugt, # in die spaeter die einzelnen Positionen der Rakete eingefuegt werden Zeit = 180 # die laenge eines Zeitschritts in Sekunden """ Die folgenden Schleifen legen den Ablauf der Flugbewegungen der Rakete fest. In den ersten 15 Zeitschritten wird der Ort der Rakete durch "raketenleapfrog" bestimmt, allerdings noch ohne Schubkraefte. Danach wird in die Variable "ortr" die aktuelle Position der Rakete geschrieben. Diese wird dann der Liste "raketenbahn" angehangen. Mit "rakete.stufe_starten(1)" wird dann einmalig die erste Stufe der Rakete gestartet. In den naechsten 400 Zeitschritten wird fast das gleiche gemacht wie in der vorherigen Schleife, nur das jetzt die Methode "treibstoffverbrauch" zusaetzlich aufgerufen wird, die den Schub zuruecksetzt sobald der Treibstoff verbraucht ist. Dies ist schon nach wenigen Zeitschritten der Fall. Danach wird die zweite Stufe gestartet und das ganze wird wiederholt. """ for i in xrange(15): rakete.ort = rakete.raketenleapfrog(Zeit,welt) ortr = ax.plot3D(np.array([rakete.ort[0]]),np.array([rakete.ort[1]]),np.array([rakete.ort[2]]),'kv') raketenbahn.append(ortr) rakete.stufe_starten(1) for i in xrange(400): rakete.ort = rakete.raketenleapfrog(Zeit,welt) ortr = ax.plot3D(np.array([rakete.ort[0]]),np.array([rakete.ort[1]]),np.array([rakete.ort[2]]),'kv') raketenbahn.append(ortr) rakete.treibstoffverbrauch(Zeit) rakete.stufe_starten(2) for i in xrange(400): rakete.ort = rakete.raketenleapfrog(Zeit,welt) ortr = ax.plot3D(np.array([rakete.ort[0]]),np.array([rakete.ort[1]]),np.array([rakete.ort[2]]),'kv') raketenbahn.append(ortr) rakete.treibstoffverbrauch(Zeit) #--------------Grafik-------------------------------------------------------------------------- #Erde ----------------------------------------------------------------------------------------- #Hier wird die runde Erde erzeugt, mit einem Radius von 6.000.000 Metern, was annaehernd realistisch ist u = np.linspace(0, 2 * np.pi, 26) v = np.linspace(0, np.pi, 14) x = 6000000 * np.outer(np.cos(u), np.sin(v)) y = 6000000 * np.outer(np.sin(u), np.sin(v)) z = 6000000 * np.outer(np.ones(np.size(u)), np.cos(v)) ax.plot_surface(x, y, z, rstride=1, cstride=1, color='b', shade=0) #Erde Ende ------------------------------------------------------------------------------------ AnimationRakete = animation.ArtistAnimation(fig, raketenbahn, 100) # Aus der Liste mit den Raketenpositionen wird nun eine Animation gemacht # Wird diese Zeile weggelassen werden die einzelnen Positionen gleichzeitig und dauerhaft angezeigt #Timer ---------------------------------------------------------------------------------------- zeit0 = time.time() # die aktuelle Zeit in Sekunden, von der die Differenz berechnet wird def update_title(axes): # eine Methode, die die verstrichene Zeit seit "zeit0" ausgibt timetext.set_text((time.time() - zeit0)*30) fig.canvas.draw() timetext = ax_head.text(0.2,0.04,'0.0') # der erste Wert, 0.0 wird festgelegt ax_head.text(0.05,0.04,'Zeit in min:') ax_head.text(0.4,0.04,'(sehr ungenau)') # leider wahr timer = fig.canvas.new_timer(interval=1000 ) # jede Sekunde wird der Timer aktualisiert timer.add_callback(update_title, ax) # die Aktualisierung greift auf die Methode zu, die oben definiert wurde timer.start() #Timer Ende ----------------------------------------------------------------------------------- grenze = 40000000 # damit die Grafik nicht verzerrt aussieht, werden hier die Grenzen gleich gross festgelegt ax.set_xlim(-grenze,grenze) ax.set_ylim(-grenze,grenze) ax.set_zlim(-grenze,grenze) plt.show() # die Grafik wird angezeigt