Mathesis Wiki

Die graphische Umsetzung

import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
 
fig = plt.figure()
ax = fig.gca(projection='3d')
 
 
class Himmelskoerper(object):
 
	Liste_der_Kraefte = []		#die Liste der Kraefte
	m = 1.0
	v = np.array([0.0,0.0,0.0])
	ort = np.array([0.0,0.0,0.0])
	name = ""
 
	def __repr__(self):
		return self.name
 
	def __init__(self,name,m):
		self.name = name
		self.m = m
 
 
	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):
		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):
		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):
		r = self.berechne_resultierende_Kraft(welt)
		a = r/self.m
		return a
 
	def berechne_geschwindigkeitsvektor(self,zeitschritt,welt):
		a = self.berechne_beschleunigung(welt)
		v = a*zeitschritt
		return self.v + v					#alte und neue Geschwindigkeit addiert
 
	def berechne_ort(self,zeitschritt,welt):
		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):
		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 Planet(Himmelskoerper):
 
	Abstand_zur_Sonne = 0
	min_Temperatur = 0
	max_Temperatur = 0
	Umlaufperiode = 0
	Rotationsperiode = 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 = Planet("Erde",5.972E24)
mond = Planet("Mond",7.35E22)
mond2 = Planet("Mond2",7.35E22)
welt = Welt()
welt.objekt_hinzufuegen(erde,np.array([0.0,0.0,0.0]))
welt.objekt_als_satellit_hinzufuegen(erde,mond,384400000.0)
welt.objekt_als_satellit_hinzufuegen(erde,mond2,384400000.0)
mond2.ort = np.array([0,mond.ort[0],0])
mond2.v = np.array([0,0,mond.v[1]])
abstand_0 = mond.ort-erde.ort
print "Erde: "+ str(erde.ort)
print "Mond: "+ str(mond.ort)
print "Abstand zwischen Mond und Erde: "+str(np.linalg.norm(abstand_0))
 
alter_mond = mond.ort 
mondbahn = []
mondbahn2 = []
for i in xrange(80):
	mond.ort = mond.leapfrog(30000.0,welt)
	mond2.ort = mond2.leapfrog(30000.0,welt)
	ort = plt.plot(np.array([mond.ort[0]]),np.array([mond.ort[1]]),np.array([mond.ort[2]]),'k.')
	mondbahn.append(ort)
	ort2 = plt.plot(np.array([mond2.ort[0]]),np.array([mond2.ort[1]]),np.array([mond2.ort[2]]),'k.')
	mondbahn2.append(ort2)
abstand_n = mond.ort-erde.ort
w = np.arccos(np.dot(abstand_0,abstand_n)/(np.linalg.norm(abstand_0)*np.linalg.norm(abstand_n)))
print "Winkel: "+str(w*360/(2*np.pi))
print "Erde: "+ str(erde.ort)
print "Mond: "+ str(mond.ort)
print "Abstand zwischen Mond und Erde: "+str(np.linalg.norm(abstand_n))
 
 
#--------------Grafik----------------------------------------------------------------------------------------------------------------
 
"""
plt.ion()
plt.plot(np.array([erde.ort[0]]),np.array([erde.ort[1]]),np.array([erde.ort[2]]),'bo')
plt.plot(np.array([mond.ort[0]]),np.array([mond.ort[1]]),np.array([mond.ort[2]]),'wo')
#plt.plot(np.array([alter_mond[0]]),np.array([alter_mond[1]]),np.array([alter_mond[2]]),'wo')
ax.set_xlim(-400000000,400000000)
ax.set_ylim(-400000000,400000000)
ax.set_zlim(-400000000,400000000)
plt.show()
"""
plt.plot(np.array([erde.ort[0]]),np.array([erde.ort[1]]),np.array([erde.ort[2]]),'bo')
Animation1 = animation.ArtistAnimation(fig, mondbahn, 1)
Animation2 = animation.ArtistAnimation(fig, mondbahn2, 1)
ax.set_xlim(-400000000,400000000)
ax.set_ylim(-400000000,400000000)
ax.set_zlim(-400000000,400000000)
 
plt.show()