Bis jetzt sind wir: Anton und Philipp.

Idee: Position von Körpern des Sonnensystems vorhersagen.

Ziel:
a) die nächste Mondfinsternis berechnen

alternativ, falls das zu schwer ist
b) berechnen, wann andere Planeten von der Erde aus sichtbar sind
c) berechnen, wann man von wo die ISS sieht

spontane Inspirationen:
[[https://de.wikipedia.org/wiki/Mondfinsternis#Schwierigkeiten_der_Vorausberechnung]]
[[https://de.wikipedia.org/wiki/Mondfinsternis#Schwierigkeiten_der_Vorausberechnung]]
[[http://www.diesonnenseite.de/startseite2/finsternisse_allg/berechnungen_zur_mondfinsternis/pdf/mond.pdf ]]

[[ws1516:himmelskoerper:Protokoll]]


**Grundlagen des Projekts**

In diesem Labor soll versucht werden eine Mondfinsternis so akkurat wie möglich vorherzusagen. Um dies zu Erreichen soll das System Erde-Sonne-Mond mittels einer Python Simulation berechnet werden um so die relative Positionen dieser Himmelskörper zu jedem Zeitpunkt zu ermitteln. Dabei müssen die zugrunde liegenden Differentialgleichungen gelöst werden und daraufhin geprüft werden ob die neu errechneten Positionen die Kriterien für eine Mondfinsternis erfüllen.

**Differentialgleichungen**
Differentialgleichung sind Gleichungen, die nicht einen Wert als Unbekannte haben sondern eine Funktion. So ist in unserem Fall die Veränderung der Geschwindigkeit über die Zeit von der Beschleunigung abhängig. Von der Geschwindigkeit ist nun wieder die Veränderung der Position abhängig. Also müssen wir in jeden Zeitschritt zwei Differentialgleichungen lösen. Da Differentialgleichung in den meisten und auch in unserem Fall nicht analytisch exakt lösbar sind muss hier ein numerisches Verfahren angewendet werden.

Das heißt, dass ein Algorithmus benötigt wird der in jedem Zeitschritt immer den direkt folgenden Wert aus dem Ursprungswert berechnet. Dies kann man durch verschiedene Verfahren bewerkstelligen, wir haben in der eigentlichen Simulation das Runge-Kutta Verfahren 4.Ordnung genutzt. Zur Übung haben wir aber auch des einfachere, aber ungenauere Euler Verfahren implementiert

Im Euler Verfahren wird grundlegend die Ausgangsfunktion genommen und mit dem Zeitschritt multipliziert und das Ergebnis auf den Ausgangswert aufaddiert.

Das Runge-Kutta Verfahren 4.Ordung hingegen splittet die Gleichung zunächst in 4 Differenzenquotienten aufgeteilt und diese anschließen zusammenaddiert. Sie erfordert einen höheren Rechenaufwand bietet bei gleichem Zeitschritt aber ein deutlich genaueres Ergebnis. 

Hier wird die Differentialgleichung x' = x mit den beiden Verfahren im Intervall 0-3,5 und dt = 0,175 mit den beiden Verfahren gelöst. Dies hat als analytische Lösung die bekannt e-Funktion und so lässt sich leicht die Genauigkeit überprüfen.

<code python>
import numpy as np
import matplotlib.pyplot as plt


def euler (f, startwert, intervallende, anzahlschritte, x): # Euler Verfahren
	h = (intervallende - startwert)/anzahlschritte # berechnet den Zeitschritt abhängig von der den Anzahl der Schritte und dem Intervall
	tpoints2 = np.linspace(startwert, intervallende, anzahlschritte) # legt die Werte der x-Achse fest. teilt das Intervall in gleich große Schritte auf
	xpoints2 = [] # Liste für die Funktionswerte
	for t in tpoints2:
		xpoints2.append(x) # erweitert xliste in jedem Schritt um den jeweils nächsten x-Wert
		x += h * f(t,x) # berechnet den neuen Funktionswert nach Euler 
	return tpoints2, xpoints2


def rk4(f, startwert, intervallende, anzahlschritte, x): # Runge-Kutta-Verfahren 4.Ordnung
	h = (intervallende - startwert)/anzahlschritte # genau wie bei Euler
	tpoints = np.linspace(startwert, intervallende, anzahlschritte)
	xpoints = []
	for t in tpoints:
		xpoints.append(x)
		k1 = h*f(t,x) # berechnet die Teilfunktionen nach Runge-Kutta
		k2 = h*f(t+0.5*k1,x+0.5*h)
		k3 = h*f(t+0.5*k2,x+0.5*h)
		k4 = h*f(t+k3,x+h)
		x += (k1+2*k2+2*k3+k4)/6 # fasst die 4 Teilfunktionen zusammen und berechnet neuen Funktionswert
	return tpoints, xpoints
 
def f(t, x): # Funktion x' = x
    return x

vt, vx = rk4(f, 0., 3.5, 20, 1.) # legt Intervall, Anzahl der Schritte und Ausgangswert fest
vx1, vt1 = euler (f, 0., 3.5, 20, 1.)



plt.plot(vx1, vt1)  
plt.show() # zeigt das geplottete Ergebnis für Euler
plt.plot(vt, vx)
plt.show() # zeigt das geplottete Ergebnis für Runge-Kutta 4.Ordnung
</code>

Euler: {{:ws1516:figure_1.png?400|}} 

Hier erreicht die Funktion bei x = 1 den Funktionswert 2,4 was noch deutlich vom korrekten Wert e = 2,71... abweicht.

Runge-Kutta 4.Ordnung: {{:ws1516:figure_2.png?400|}} 

Hier erreicht die Funktion bei x = 1 einen Funktionswert von 2,52 was deutlich näher an dem gesuchten Wert liegt.

**Vpython**

Um eine graphische Ausgabe zu erhalten haben wir mit der Python Erweiterung VPython gearbeitet. In dieser Erweiterung ist eine Vektorklasse implementiert die in einem 3D-Modul Objekten Positionen, Form, Farbe etc übergibt und dieses dann darstellt. Als Übung haben wir zunächst ein kleines 
Programm mit sich statisch bewegenden Kugeln geschrieben.

<code python>
from visual import *
import numpy as np

floor = box (pos=(0,0,0), length=50, height=0.5, width=4, color=color.blue)
ball = sphere (pos=(0,4,0), radius=1, color=color.red) # Ball wird erzeugt
ball.velocity = vector(5,-1,0) # Geschwindigkeit des Balls wird festgelegt
dt = 0.01
wall_left= box (pos= (25,5,0), length=1, height=10, widht=4, color=color.blue) # Begrenzungen werden festgelegt
wall_right= box (pos= (-25,5,0), length=1, height=10, widht=4, color=color.blue)

schnee=[] # List mit diversen kleinen "Schneekugeln"
for i in range (5):
	schnee.append(sphere(pos= ((-25-25)*np.random.rand()+25,10,0) , radius=0.1 , color=color.white)) # Position wird partiell zufaellig erzeugt
	
for flocke in schnee:
	flocke.velocity= vector (0.,-1.,0.) # Geschwindigkeit wird für jede Schneeflocke festgelegt
	flocke.timer=1
		
while 1: # Simulation läuft endlos
	rate (100) # 100 Schleifendurchläufe werden pro Sekunde durchgeführt
	ball.pos = ball.pos + ball.velocity*dt # Position des Balls wird nach Euler geupdated

	if ball.y < ball.radius:
		ball.velocity.y = abs(ball.velocity.y)
	else:
		ball.velocity.y = ball.velocity.y - 9.8*dt # Verhalten des Balls im an den Rändern der erzeugten Box wird festgelegt

	if ball.x > 25-ball.radius:
		ball.velocity.x = -(ball.velocity.x)
	if ball.x < -25+ball.radius:
		ball.velocity.x = -(ball.velocity.x)
	
	if len(schnee)>0 and schnee[-1].y < 9: # Berechnung der Schneeflocken
		for i in range (5):
			schnee.append(sphere(pos= ((-25-25)*np.random.rand()+25,10,0) , radius=0.1 , color=color.white))
			schnee[-1].velocity=vector (0,-1,0)
			schnee[-1].timer=1
	zu_loeschend=[]
	for i,flocke in enumerate(schnee):
		#flocke.velocity= vector (0,-5,0)
		flocke.pos= flocke.pos+flocke.velocity*dt
		flocke.timer=flocke.timer+1
		if flocke.y < flocke.radius+0.18: # legt fest wenn die Scheeflocke den "Boden" erreicht hat
			flocke.velocity=vector(0.,0.,0.)
		else:
		    flocke.velocity= vector (0.,-1.,0.)
		if flocke.timer > 500: # loescht Schneeflocken wenn ihre Anzahl 500 übersteigt
			zu_loeschend.append(i)
			
	schnee_neu=[ f for (i,f) in enumerate(schnee) if i not in zu_loeschend]
</code>

**Simulation des Sonnensystems**

Um das S