Ellipse @ PS-Trainer

Ellipse
MathML Die auf dieser Seite verwendeten Symbole orientieren sich vorwiegend an der Wikipedia-Seite zur Ellipse. Allerdings sind die verwendeten Bezeichnungen nicht einmal innerhalb der einzelnen Wiki-Seiten einheitlich.	Deshalb werden zukünftig die aus Wikipedia bezogenen Formel-Bilder von hier verschwinden und durch einheitliche Formeln in der Standard Formel-'Sprache' MathML (W3C, Wikipedia) ersetzt, die auch weniger Speicherplatz brauchen, dafür aber nur mit modernen Browsern angezeigt werden. ♦ Details zu MathML, mit Live-Test
Die → SVG-Grafik (rechts) zeigt die wichtigsten Elemente einer Ellipse mit den meist verwendeten Bezeichnungen. Sie wird von allen modernen Browsern angezeigt. M...Mittelpunkt F1, F2...Brennpunkte a...Große Halbachse b...Kleine Halbachse e = sqrt(a^2 - b^2) Lineare Exzentrizität ε = e/a = sqrt(a^2-b^2)/a Numerische Exzentrizität p = b^2/a Halbparameter Analytische Geometrie: x^2/a^2 + y^2/b^2 = 1 Fläche: A = a * b * Pi Umfang: nächstes ↓ Kapitel	Quelle: Wikipedia

Umfang einer Ellipse

Der Umfang einer Ellipse lässt sich nicht mit einer einfachen Formel berechnen. Eine Möglichkeit ist diese Reihen-Entwicklung:

Quelle: Wikipedia

a ... große Halbachse der Ellipse
b = a * sqrt(1 - ε^2) ... kleine Halbachse der Ellipse
ε = sqrt((a^2-b^2)/a)... numerische Exzentrizität, (in den Funktionen mit e[num] bezeichnet), eine dimensionslose Zahl im Bereich 0...ε...1 welche die Form der Ellipse beschreibt. Der Grenzfall ε=0 ergibt einen Kreis.
• Eine Besonderheit ist die Verwendung von 2 verschachtelten Schleifen: Innen ein Produkt, außen eine Summe.

Rotations-Ellipsoide:

Quelle: Wikipedia

Fast alle größeren Himmelskörper kann man mit guter Näherung als Rotations-Ellipsoide beschreiben:

• Ein Querschnitt durch den Äquator ist kreisförmig, die Berechnung des Äquator-Umfangs daher relativ einfach.

• Ein Querschnitt durch die Pole hat die Form einer Ellipse und ist schwieriger zu berechnen. Man kann dazu u.a. den hier vorgestellten Algorithmus verwenden.

Beispiel-Daten

Als Beispiel werden die Daten der Erde nach dem derzeit meist verwendeten → Erd-Modell WGS84 verwendet:

a = 6378137m [Mittelpunkt - Äquator]
b = a * (1-f) = 6356752.31424518m [MiPu - Pol]

Meist wird nicht b angegeben, sondern die Abflachung

f = 1/298.257223563

Zum Vergleich der Radius der Erde nach dem einfachen Kugel-Modell:

r[E] = 6378000m

Ergebnisse:

u[kugel] = 2 * r[E] * Pi = 40074155.9m
u[eq] = 2 * a * Pi = 40075016.7m
u[pol] = 40007862.9m

Der Umfang der Erde ist am Äquator um ca. 67154m länger als über die Pole.

Kalkulations-Programm

Wenn man ein paar kleine Tricks zulässt, dann ist die Programmierung einfacher als die Formel vermuten lässt:

	A	B	C	D
1	a	6378137
2	b	=B1*(1-B3)
3	f	=1/298,257223563
4	e[num]	=WURZEL(B1^2-B2^2)/B1
5	...
6	Element	Summe	Produkt	dif[m]
7	1	=A7	=2$B$1PI()*B7
8	=(1/2)^2*($B$4^2)/1	=B7-A8	=2$B$1PI()*B8	=C8-C7
9	=(1/23/4)^2($B$4^4)/3	=B8-A9	=2$B$1PI()*B9	=C9-C8
10	=(1/23/45/6)^2*($B$4^6)/5	=B9-A10	=2$B$1PI()*B10	=C10-C9
11	=(1/23/45/67/8)^2($B$4^8)/7	=B10-A11	=2$B$1PI()*B11	=C11-C10

Die Formeln der Zellen B8; C7; D8 kann man wie üblich nach unten ausfüllen. Zur Vereinfachung ist jedoch die Berechnung der einzelnen Elemente in Spalte A nicht allgemein sondern speziell gelöst: Jede Zelle enthält eine andere Formel.
Fortgeschrittene AnwenderInnen können versuchen, auch diese Spalte mit einer allgemeinen Formel zu berechnen.
Das Zwischen-Ergebnis wird in Spalte C berechnet (End-Ergebnis = Umfang der Ellipse in C11). Aus der Differenz der aufeinander folgenden Zwischen-Ergebnisse (Spalte D) erkennt man, dass die Reihen-Entwicklung sehr rasch konvergiert: Die Differenz beträgt - für den Erdumfang - in D8 noch 67km, in D11 nur mehr <1mm

Programmierung mit Basic (Libreoffice-Basic, Visual Basic VBA)

• Die Kreiszahl Pi wird als Konstante c_pi vorgegeben.

• Die Funktion erwartet die beiden Halbachsen a,b der Ellipse als Argumente.

• Im ersten Schritt wird die numerische Exzentrizität ε = e_num berechnet.

• Die Angabe der maximalen Anzahl von Iterationen imax ist nicht notwendig, aber (wie bei allen Reihen-Entwicklungen) eine empfohlene Maßnahme zur Sicherheit gegen endlose Schleifen.

• Initialisierung: Vor Beginn der Schleife werden den Variablen ihre Anfangswerte zugewiesen.

• Die While-Schleife berechnet in jedem Durchlauf 1 Element element der Reihe. Dazu braucht man das in der Formel angegebene Produkt, welches in einer For-Schleife berechnet wird. Wie man aus den Formeln der Kalkulation (Spalte A) erkennt, werden die Zahlen k abwechselnd multipliziert und dividiert. Das wird hier durch die Bool-Variable mul gesteuert, die in jedem Durchlauf der For-Schleife umgekehrt wird.

• Nach der Berechnung eines Elements wird sein Wert von der laufenden Summe summe abgezogen.

• Abbruch-Test: Wenn sich das Ergebnis nicht mehr geänderrt hat (vsumme=summe) oder wenn die maximale Anzahl von Elementen berechnet wurde (hier imax=1000) dann wird doloop=False gesetzt und die while-Schleife abgebrochen.
Der aktuelle Wert von summe wird für den Abbruch-Test des nächsten Durchgangs in vsumme gespeichert.

• Nach dem Ende der While-Schleife wird das Ergebnis berechnet und zurückgegeben. Es ist nicht notwendig, das (Zwischen)-Ergebnis in jedem Durchgang der Schleife zu berechnen, man kann das jedoch zu Debugging-Zwecken durchführen.

Basic-Funktion zur Berechnung des Umfangs einer Ellipse:

Const c_pi = 3.14159265358979

Function ellipse_umfang(a As Double, b As Double) As Double
Dim doloop, mul As Boolean
Dim i, k, imax As Integer
Dim e_num, element, summe, vsumme As Double

e_num = Sqr((a ^ 2 - b ^ 2) / a ^ 2)
imax = 1000

'--- Initialisierung ---

summe = 1
vsumme = 1
i = 1
doloop = True

'--- Schleife ---

While doloop

vsumme = summe
mul = True
element = 1
For k = 1 To 2 * i

If mul Then

element = element * k

Else

element = element / k

End If
mul = Not mul

Next k
element = element ^ 2 * (e_num ^ (2 * i)) / (2 * i - 1)
summe = summe - element
i = i + 1
If (vsumme=summe Or i>imax) Then doloop = False

Wend
ellipse_umfang = 2 * a * c_pi * summe

End Function

Anwendung der Basic-Funktion in einer Kalkulation:

=ellipse_umfang(B1;B2)

Planeten-Bewegung (2. Kepler-Gesetz)

• Das 1. Kepler-Gesetz beschreibt die Bahnen der Planeten als Ellipsen, in deren einem Brennpunkt die Sonne ● steht.

• Das 2. Kepler-Gesetz beschreibt die Bewegung eines Planeten auf seiner Bahn: Sie verläuft ungleichmäßig: Am Sonnen-nächsten Punkt (Perihel) ist die Geschwindigkeit am größten, am Sonnen-fernsten Punkt (Aphel) am kleinsten.

• Das Demo-Beispiel (rechts) zeigt eine → SVG-Grafik mit → Javascript-Animation. Ein Demo-Planet umkreist die Sonne im Erd-Abstand (1 AU), jedoch auf einer deutlich elliptischen Bahn e[num]=0.5, damit die Änderung der Geschwindigkeit erkennbar ist.
Klick auf ● zeigt Bahnpunkte, die in gleichem zeitlichen Abstand liegen.
Klick auf ┼ zeigt Achsen und Winkel an.
Klick auf ▼ ▲ ändert die Geschwindigkeit der Animation.
Die Zeit wird in Tagen angegeben.

• Das Beispiel ist eine starke Vereinfachung:
Die Dimensionen von Sonne, Bahn und Planet sind stark verändert.
Die Erde erreicht ihren Sonnen-nächsten Punkt nicht mit Neujahr, sondern ca. am 81. Tag

Bahn-Daten einiger Körper des Sonnensystems (ohne Gewähr)
Achsen in AE: 1 AE = 149.6 Mio km
a...Große Halbachse (semi-major axis)
b...Kleine Halbachse (semi-minor axis)
e[num]...Numerische Exzentrizität (0...e[num]...1)
T[d]...(Siderische) Umlaufzeit in Erd-Tagen (1 d = 86400 s)

Die Daten stammen aus unterschiedlichen Quellen, u.a. Wikipedia. Man findet meist die Angaben für a,e,T und berechnet daraus alle anderen Daten, z.B. diese Entfernungen:

b = a * sqrt(1 - e[num]^2)
e[lin] = a * e[num]
Aphel = a + e[lin]
Perihel = a - e[lin]

Für einen Demo-Körper kann man die Umlaufzeit aus a nach dem 3.Kepler-Gesetz berechnen:

C[kep] = T[erde]^2 / a[erde]^3 = 133412
T[demo] = sqrt(C[kep] * a[demo]^3)

Alle Werte wurden zwecks besserer Übersicht gerundet.

Körper	a	b	e[num]	T[d]
Merkur	0.387	0.379	0.2056	87.969
Venus	0.723	0.723	0.0068	224.701
Erde	1.000	1.000	0.0167	365.256
Mars	1.524	1.517	0.0935	686.980
Jupiter	5.203	5.197	0.0484	4331.936
Saturn	9.582	9.568	0.0542	10759.346
Uranus	19.201	19.180	0.0472	30685.522
Neptun	30.047	30.045	0.0113	60190.536
Pluto	39.482	38.240	0.2488	90466.606
Halley	17.834	4.544	0.9670	27511.082
Demo ↑	1.000	0.866	0.5000	365.256

• Die Halbachsen a,b der Bahnen sind in der Einheit AU (Astronomical Units) angegeben, die mit der großen Halbachse der Erde (rot) definiert ist. Die beiden Achsen von Venus und Erde sind nicht gleich groß, die Differenzen sind jedoch kleiner als 0.001 AU.
• Die numerische Exzentrizität e[num] ist dimensionslos.
• Die Umlaufzeiten sind in Erd-Tagen angegeben.
• Man kann man die Funktionen ohne Änderung auch mit anderen Einheiten verwenden.

Berechnung der XY-Koordinaten

eines umlaufenden Körpers als Funktion der Zeit (ohne Gewähr).
Vor Beginn definiert man zweckmäßig die Werte einiger globalen Konstanten:

c_pi = 3.14159265358979 (→ Kreiszahl Pi)
c_pi_2 = 1.5707963267949 (rechter Winkel 90°)
c_3_pi_2 = 4.71238898038469 (rechter Winkel 270°)

Man definiert einige globale Variable, deren Werte (Bahn-Daten) für einen bestimmten Körper konstant sind, und setzt Vorgabe-Daten bzw. berechnete Daten ein:

a = ... (große Halbachse der Bahn)
b = ... (kleine Halbachse der Bahn)
T = ... (Umlaufzeit)

Der Algorithmus wurde für die Programmiersprache → Basic (LibreOffice-Basic, Visual Basic, VBA) in 3 Teile (einzelne Funktionen) zerlegt:
• Berechnung des Bahnwinkels w aus a,b,T,t (Zeit)
• Berechnung der X-Koordinate aus a,b,w
• Berechnung der Y-Koordinate aus a,b,x

In modernen Programmiersprachen kann man mehrere Ergebnisse als Array zurückgeben und in diesem Fall auf die Koordinaten (x,y) beschränken.
Unabhängig davon kann man auf alle mathematischen Funktionen verzichten und die Funktionen selbst mit Reihen-Entwicklungen berechnen.

● Alle Beispiele wurden mit einem Demo-Körper (letzte Zeile der Tabellen-Daten) berechnet: Diese Bahn wurde absichtlich mit e[num]=0.5 wesentlich stärker elliptisch angenommen als reale Planeten-Bahnen, damit man den Effekt besser erkennt.

Berechnung des Bahnwinkels

Diese Rechnung wird als Reihen-Entwicklung ausgeführt:

► Zunächst werden die Werte der Argumente geprüft und bei Bedarf geändert.

► Man berechnet als 1.Näherung den Bahnwinkel in einer Kreisbahn:

w[0] = w[k] = 2 * c_pi * t / T

Beispiel für t=10 Tage:

w[k] = 2 * 3.1416 * 10 / 365.256 = 0.31476

Winkel werden in der Informatik immer im Boganmaß verwendet, es entspricht

w[k] = 0.31476 → 0.31487/Pi*180 = 18.034° (im Gradmaß)

► Die 2.Näherung wird aus der ersten berechnet:

w[1] = w[k] + e[num] * sin(w[0])

Beispiel:

w[1] = 0.31476 + 0.5 * sin(0.31476) = 0.46955

► Alle weiteren Näherungen werden jeweils aus der vorhergehenden berechnet:

w[n+1] = w[k] + e[num] * sin(w[n])

► Die Reihe wird abgebrochen, wenn man eine gewünschte → Genauigkeit erreicht hat - spätestens dann, wenn sich das Ergebnis nicht mehr ändert (und daher auf ca. 15 Stellen genau ist). Das ist bei diesem Algorithmus nach ca. 30-50 Schritten der Fall:

Beispiel:

w[49] = w[50] = 0.595 (→ 34,09°)

Basic-Funktion zur Berechnung des Bahnwinkels:

Const c_pi = 3.14159265358979
Const c_pi_2 = 1.5707963267949
Const c_3_pi_2 = 4.71238898038469

Function kepler_w( _

Optional e_num As Double = 0, _
Optional T_Orbit As Double = 1, _
Optional t_calc As Double = 0, _
Optional nmax As Integer = -1) As Double

Dim doloop As Boolean
Dim w, wk, wv As Double
Dim n As Integer

' Argumente

If (e_num < 0) Then e_num = 0
If (e_num > 1) Then e_num = 0.9999999999
If (T_Orbit <= 0) Then T_Orbit = 1
While t_calc < 0

t_calc = t_calc + T_Orbit

Wend
While t_calc > T_Orbit

t_calc = t_calc - T_Orbit

Wend
If (nmax < 3) Then nmax = 3
If (nmax > 100) Then nmax = 100

' Anfangswert (Kreis)

wk = 2 * c_pi * t_calc / T_Orbit
wv = wk
doloop = True
iloop = 0

' Schleife

While doloop

wv = w
w = wk + e_num * Sin(wv)
n = n + 1

' Test

If (w=wv Or n >=nmax) Then doloop = False

Wend

' Ergebnis

kepler_w = w

End Function

Berechnung der X-Koordinate

x = sqrt( 1 / (1/a^2 + tan^2(w)/b^2) )

Für Winkel im Bereich 90°...w...270° muss man das Vorzeichen korrigieren. Diese Winkel werden als globale Konstanten (Absatz oben) vorgegeben.

Beispiel:

w[10d] = 0.595 nach 10 Tagen
x = 0.92675 AU (Astronom.Units)

Basic-Funktion zur Berechnung der X-Koordinate:

Function kepler_x( _

a As Double, _
b As Double, _
w As Double) As Double

Dim x As Double

x = Sqr(1 / ((1 / a ^ 2) + ((Tan(w)) ^ 2 / b ^ 2)))
If (w > c_pi_2 And w < c_3_pi_2) Then x = -x
kepler_x = x

End Function

Berechnung der Y-Koordinate

Dazu gibt es mehrere Alternativen:

y = x * tan(w)
y = sqrt( b^2 * (1 - x^2/a^2) )

In der 2.Variante muss man das Vorzeichen für Winkel im Bereich >180° korrigieren. Nur dazu wird auch der Bahn-Winkel w als Argument übergeben.

Beispiel:

y[10d] = 0.3253 AU

Basic-Funktion zur Berechnung der Y-Koordinate:

Function kepler_yy(x As Double, w As Double) As Double

kepler_yy = x * Tan(w)

End Function

Alternative Berechnung der Y-Koordinate:

Function kepler_y( _

a As Double, b As Double, _
x As Double, _
w As Double) As Double

Dim y As Double

y = Sqr(b ^ 2 * (1 - x ^ 2 / a ^ 2))
If (w > c_pi) Then y = -y
kepler_y = y

End Function

Anwendung in einem Kalkulations-Programm:

Berechnung der XY-Koordinaten als Funktion der Zeit. Als Längen-Einheit wird 1 AU = 149.6 Mio km verwendet, als Zeit-Einheit Tage.

	A	B	C	D
1	a	1
2	b	=B1*WURZEL(1-B3^2)
3	e[num]	0,5
4	T	=WURZEL(133412*B1^3)
5
6	t	w	x	y
7	0	=kepler_w($B$3;$B$4;A7)	=kepler_x($B$1;$B$2;B7)	=kepler_y(C7;B7)
8	=A7+10	=kepler_w($B$3;$B$4;A8)	=kepler_x($B$1;$B$2;B8)	=kepler_y(C8;B8)
...
43	=A42+1	=kepler_w($B$3;$B$4;A42)	=kepler_x($B$1;$B$2;B42)	=kepler_y(C42;B42)

• Sie können die Formel von Zelle A8 sowie der Zellen B7:D7 nach unten bis Zeile 43 ausfüllen.
• Wenn sie aus dem Bereich C7:D43 ein XY-Diagramm erzeugen, dann sollte es so ähnlich aussehen wie die Demo-Grafik am Beginn dieses Kapitels. Die Abstände der berechneten Punkte wurden mit 10 Tagen absichtlich groß gewählt, damit man die Geschwindigkeit an verschiedenen Punkten der Bahn ohne Hilfsmittel erkennen kann: Je weiter die Bahn-Punkte auseinander liegen, desto größer die Geschwindigkeit.
• Die Sonne befindet sich am Brennpunkt mit den Koordinaten x[sol]=e[num]; y[sol]=0

Javascript

JS-Versionen der vorgestellten Programme finden sie im Live-Beispiel, welches am ↑ Beginn dieses Kapitels gezeigt wird.
Mit Rechtsklick in die Grafik können sie den Quelltext anzeigen.
• SVG-Grafik wird von allen modernen Browsern seit vielen Jahren problemlos angezeigt.

Javascript funktioniert mit jedem gängigen Browser, arbeitet sehr schnell und lässt sich zusammen mit allen Mitgliedern der → XML-Familie verwenden, z.B. mit → XHTML (diese Webseite) oder → SVG (Demo-Grafik).

Algorithmen	Ausgewählte IT-Rezepte, Reihen-Entwicklung (Iteration)
Ellipse	Allgemeine Angaben
Umfang	Berechnung des Umfangs einer Ellipse mit Reihen-Entwicklung
Planeten-Bahn	Berechnung der XY-Koordinaten von Planeten als Funktion der Zeit (2.Kepler-Gesetz)

Ellipse

Berechnung einiger Aufgaben mit Reihen-Entwicklung

Ellipse

MathML

Umfang einer Ellipse

Rotations-Ellipsoide:

Kalkulations-Programm

Programmierung mit Basic (Libreoffice-Basic, Visual Basic VBA)

Planeten-Bewegung (2. Kepler-Gesetz)

Berechnung der XY-Koordinaten

Berechnung des Bahnwinkels

Berechnung der X-Koordinate

Berechnung der Y-Koordinate

Anwendung in einem Kalkulations-Programm:

Javascript