std.astro.cr3bp — Librationspunkte und Jacobi-Konstante

Das eingeschränkte Dreikörperproblem: ein masseloser Körper im Feld zweier großer Massen, die einander umkreisen. Diese Unit steht neben dem übrigen std.astro und nicht darin — alles andere rechnet im Zweikörperproblem, und in dem existieren Librationspunkte überhaupt nicht (#1997).

std.astro (Übersicht) · std.astro.orbit · std.astro.state · Bahntypen


1. Das Koordinatensystem

Gerechnet wird dimensionslos im mitrotierenden System: der Abstand der beiden Massen ist 1, ihre Summe ist 1, die Winkelgeschwindigkeit ist 1. Der Ursprung liegt im gemeinsamen Schwerpunkt, die große Masse bei −mu, die kleine bei 1−mu.

Der praktische Vorteil: dieselben Zahlen gelten für Sonne–Erde und Erde–Mond. ADimToMeters rechnet mit dem tatsächlichen Abstand zurück, ADistanceToSecondary liefert direkt den Abstand zur kleineren Masse in Metern.

 
mu ist m2 / (m1 + m2) mit m2 als der KLEINEREN Masse. ACr3bpMu sortiert nicht selbst — wer die Massen vertauscht, bekommt ein mu über 0.5 zurück. Das ist eine ehrliche Antwort auf eine vertauschte Eingabe, kein Fehler; die Punktrechnungen prüfen den Bereich.

2. Die Funktionen

Funktion Rückgabe Bedeutung
ACr3bpMu(m1, m2) f64 Massenparameter aus zwei Massen
AMuSunEarth(), AMuEarthMoon(), AMuSunJupiter() f64 die drei gebräuchlichen Systeme
AHillRadius(a, mu) f64 Hill-Radius der kleineren Masse
ALagrange1(mu)ALagrange5(mu) AVec3 die fünf Punkte, dimensionslos
ALagrange1Gamma(mu)ALagrange3Gamma(mu) f64 Abstand des kollinearen Punktes von der kleinen Masse, dimensionslos
ADimToMeters(v, separation) AVec3 dimensionslos → Meter
ADistanceToSecondary(p, mu, separation) f64 Abstand zur kleineren Masse in Metern
AJacobiConstant(mu, pos, vel) f64 die Erhaltungsgröße des Problems
AEffectivePotential(mu, pos) f64 Potential im mitrotierenden System

3. Wo die Punkte liegen

Sonne–Erde, mu = 3.003·10⁻⁶:

Punkt x y Abstand zur Erde
L1 0.990027 0 1 491 551 km
L2 1.010034 0 1 501 532 km
L3 −1.000001 0 299 195 479 km
L4 0.499997 0.866025 149 597 871 km
L5 0.499997 −0.866025 149 597 871 km

Erde–Mond, mu = 0.0121506:

Punkt x Abstand zum Mond
L1 0.836915 58 019 km
L2 1.155682 64 515 km

Drei Dinge sind an dieser Tabelle bemerkenswert:

  • L2 der Sonne–Erde ist die 1.5-Millionen-Kilometer-Marke, die man aus jeder Meldung zum James-Webb-Teleskop kennt. Hier kommt sie aus der Rechnung: 1 501 532 km.
  • L4 und L5 bilden ein gleichseitiges Dreieck mit den beiden Massen. Die y-Koordinate 0.866025 ist √3/2 — der Abstand zur Erde ist genau eine astronomische Einheit, so weit wie die Erde von der Sonne.
  • L1 und L2 liegen nicht symmetrisch. 1 491 551 gegen 1 501 532 km, knapp 10 000 km Unterschied. Der Hill-Radius der Erde liegt mit 1 496 557 km genau dazwischen.

4. Die Jacobi-Konstante

Im rotierenden System ist die Energie nicht erhalten, wohl aber die Jacobi-Konstante — die einzige bekannte Erhaltungsgröße des eingeschränkten Dreikörperproblems. Sie legt fest, welche Gebiete ein Körper überhaupt erreichen kann.

Ort (Erde–Mond, in Ruhe) Jacobi-Konstante
L1 3.188341
L4 2.987997

Der Wert in L1 ist der Schwellenwert, ab dem sich die Bahngebiete um Erde und Mond zu einem verbinden — unterhalb bleibt ein Körper bei einer der beiden Massen gefangen. Für L4 ergibt sich 3 − mu = 2.9879 — die Formel geht exakt auf.

5. Halo-Bahnen

Die periodischen Bahnen um L1 und L2, auf denen die Teleskope tatsächlich fliegen, sind seit #2001 enthalten — und zwar nicht über Richardsons Reihe dritter Ordnung, sondern über eine Fortsetzung: man löst zuerst eine kleine Bahn, hebt die Auslenkung schrittweise an und nimmt jedes Mal die vorige Lösung als Startwert. Jede Zwischenlösung ist selbst eine periodische Bahn, also ist das Ergebnis nachprüfbar.

Funktion Rückgabe Bedeutung
ACr3bpAccel(mu, pos, vel) AVec3 Beschleunigung im rotierenden System
ACr3bpStep(mu, s, dt) AState ein Integrationsschritt
ACr3bpPropagate(mu, s, dt, schritte) AState mehrere Schritte am Stück
ACr3bpBisEbene(mu, s, dt, maxSchritte) AState integriert bis zum Durchgang durch die xz-Ebene
ACr3bpC2(mu, gamma, istL2) f64 Koeffizient der Linearisierung
ACr3bpLambda(c2), ACr3bpKappa(c2, lambda) f64 Eigenwert und Amplitudenverhältnis
AHaloGuess(mu, gamma, istL2, ax, az, nordSued) AState linearisierter Startwert
AHaloCorrect(mu, s, dt, maxIter) AState Differentialkorrektur auf eine periodische Bahn
AHaloPeriod(mu, s, dt) f64 Umlaufzeit, dimensionslos
AHaloAmplitudeZ(s, separation) f64 Auslenkung aus der Ebene in Metern
AHaloFromAmplitudeZ(mu, gamma, istL2, azZiel, nordSued, dt, schritte) AState Halo-Bahn zu vorgegebener Auslenkung, per Fortsetzung

Was die Korrektur leistet, zeigt der Vergleich am selben Startpunkt: die Linearisierung schätzt vy = 0.003290, die periodische Bahn braucht 0.009846 — Faktor drei. Der linearisierte Wert allein trägt nur bei sehr kleinen Bahnen.

Auslenkung (in gamma) Startpunkt x erreichte Auslenkung Umlaufzeit
0.05 1.008354 75 077 km 3.1021 = 180.3 Tage
0.30 1.007301 450 459 km 3.0847 = 179.3 Tage
0.40 1.006267 600 613 km 3.0609 = 177.9 Tage

Die Umlaufzeit liegt bei rund einem halben Jahr und sinkt mit wachsender Auslenkung leicht — 2.4 Tage zwischen der kleinsten und der größten Bahn. Das ist die Halo-Familie: eine Schar, keine einzelne Bahn.

 
azZiel ist in Einheiten von gamma angegeben, nicht in Metern und nicht dimensionslos in Systemeinheiten. Für Sonne–Erde-L2 ist gamma = 0.01003712, eine Auslenkung von 0.4 entspricht also 600 613 km.

schritte ist die Feinheit der Fortsetzung, nicht die Zahl der Integrationsschritte. Zu grob heißt, dass ein Schritt aus dem Einzugsbereich der Korrektur fällt; die Funktion liefert dann einen Fehlerwert statt einer falschen Bahn. Als Anhalt: 0.3 braucht rund 20 Schritte, 0.4 rund 40, 0.45 bereits 200.
 
Oberhalb von 0.45 gamma bricht die Fortsetzung ab — auch mit 400 Schritten. Die größte erreichte Bahn liegt bei 675 689 km. Das James-Webb-Teleskop fliegt mit rund 800 000 km Auslenkung, und die ist damit noch nicht rechenbar. Gemeldet als #2016.

6. Die Programme

Teil A — die Punkte und die Jacobi-Konstante aus den Abschnitten 3 und 4:

unit main;
import std.astro.cr3bp;
import std.astro.state;
import std.math;
import std.io;
import std.string;

fn f(v: f64, d: int64): pchar { return FloatToStr(v, d); }
fn AU():  f64 { return 149597870700.0; }
fn EM():  f64 { return 384400000.0; }

fn punkt(name: pchar, p: AVec3, mu: f64, sep: f64): void {
    Print(name);
    Print(" | x "); Print(f(p.x, 6));
    Print(" | y "); Print(f(p.y, 6));
    Print(" | Abstand zur kleinen Masse ");
    PrintLn(StrConcat(f(ADistanceToSecondary(p, mu, sep) / 1000.0, 0), " km"));
}

fn main(): int64 {
    var muSE: f64 := AMuSunEarth();
    var muEM: f64 := AMuEarthMoon();
    PrintLn(StrConcat("mu Sonne-Erde:  ", f(muSE, 9)));
    PrintLn(StrConcat("mu Erde-Mond:   ", f(muEM, 7)));
    PrintLn(StrConcat("mu Sonne-Jupiter: ", f(AMuSunJupiter(), 9)));
    PrintLn(StrConcat("aus Massen gerechnet: ", f(ACr3bpMu(1.98847e30, 5.9722e24), 9)));
    PrintLn("");
    PrintLn("-- Sonne-Erde --");
    punkt("L1", ALagrange1(muSE), muSE, AU());
    punkt("L2", ALagrange2(muSE), muSE, AU());
    punkt("L3", ALagrange3(muSE), muSE, AU());
    punkt("L4", ALagrange4(muSE), muSE, AU());
    punkt("L5", ALagrange5(muSE), muSE, AU());
    PrintLn("");
    PrintLn("-- Erde-Mond --");
    punkt("L1", ALagrange1(muEM), muEM, EM());
    punkt("L2", ALagrange2(muEM), muEM, EM());
    PrintLn("");
    PrintLn(StrConcat("Hill-Radius der Erde:  ", StrConcat(f(AHillRadius(AU(), muSE) / 1000.0, 0), " km")));
    PrintLn(StrConcat("Hill-Radius des Mondes: ", StrConcat(f(AHillRadius(EM(), muEM) / 1000.0, 0), " km")));
    PrintLn(StrConcat("gamma L1 (Sonne-Erde): ", f(ALagrange1Gamma(muSE), 8)));
    PrintLn(StrConcat("gamma L2 (Sonne-Erde): ", f(ALagrange2Gamma(muSE), 8)));
    PrintLn("");
    var ruhe: AVec3 := AVec(0.0, 0.0, 0.0);
    PrintLn(StrConcat("Jacobi-Konstante in L4, in Ruhe: ", f(AJacobiConstant(muEM, ALagrange4(muEM), ruhe), 6)));
    PrintLn(StrConcat("Jacobi-Konstante in L1, in Ruhe: ", f(AJacobiConstant(muEM, ALagrange1(muEM), ruhe), 6)));
    PrintLn(StrConcat("effektives Potential in L4:      ", f(AEffectivePotential(muEM, ALagrange4(muEM)), 6)));
    return 0;
}

Ausgabe:

mu Sonne-Erde:  0.000003003
mu Erde-Mond:   0.0121506
mu Sonne-Jupiter: 0.000953875
aus Massen gerechnet: 0.000003003

-- Sonne-Erde --
L1 | x 0.990027 | y 0.000000 | Abstand zur kleinen Masse 1491551 km
L2 | x 1.010034 | y 0.000000 | Abstand zur kleinen Masse 1501532 km
L3 | x -1.000001 | y 0.000000 | Abstand zur kleinen Masse 299195479 km
L4 | x 0.499997 | y 0.866025 | Abstand zur kleinen Masse 149597871 km
L5 | x 0.499997 | y -0.866025 | Abstand zur kleinen Masse 149597871 km

-- Erde-Mond --
L1 | x 0.836915 | y 0.000000 | Abstand zur kleinen Masse 58019 km
L2 | x 1.155682 | y 0.000000 | Abstand zur kleinen Masse 64515 km

Hill-Radius der Erde:  1496557 km
Hill-Radius des Mondes: 61274 km
gamma L1 (Sonne-Erde): 0.00997040
gamma L2 (Sonne-Erde): 0.01003712

Jacobi-Konstante in L4, in Ruhe: 2.987997
Jacobi-Konstante in L1, in Ruhe: 3.188341
effektives Potential in L4:      1.493999

Teil B — die Halo-Bahnen aus Abschnitt 5:

unit main;
import std.astro.cr3bp;
import std.astro.state;
import std.math;
import std.io;
import std.string;

fn f(v: f64, d: int64): pchar { return FloatToStr(v, d); }
fn AU(): f64 { return 149597870700.0; }

fn bahn(name: pchar, az: f64, schritte: int64): void {
    var mu: f64 := AMuSunEarth();
    var g:  f64 := ALagrange2Gamma(mu);
    var s:  AState := AHaloFromAmplitudeZ(mu, g, true, az, 1.0, 0.001, schritte);
    Print(name);
    Print(" | x "); Print(f(s.rx, 6));
    Print(" | Az "); Print(f(AHaloAmplitudeZ(s, AU()) / 1000.0, 0));
    Print(" km | Periode "); Print(f(AHaloPeriod(mu, s, 0.001), 4));
    Print(" = "); Print(f(AHaloPeriod(mu, s, 0.001) * 365.25 / (2.0 * Pi()), 1));
    PrintLn(" Tage");
}

fn main(): int64 {
    var mu: f64 := AMuSunEarth();
    var g:  f64 := ALagrange2Gamma(mu);
    var c2: f64 := ACr3bpC2(mu, g, true);
    var lam: f64 := ACr3bpLambda(c2);
    PrintLn(StrConcat("gamma L2:           ", f(g, 8)));
    PrintLn(StrConcat("c2:                 ", f(c2, 6)));
    PrintLn(StrConcat("lambda:             ", f(lam, 6)));
    PrintLn(StrConcat("kappa:              ", f(ACr3bpKappa(c2, lam), 6)));
    PrintLn("");

    // Der linearisierte Startwert und was die Korrektur daraus macht
    var s0: AState := AHaloGuess(mu, g, true, 0.05, 0.05, 1.0);
    var s1: AState := AHaloCorrect(mu, s0, 0.001, 20);
    PrintLn(StrConcat("Startschaetzung  vy = ", f(s0.vy, 6)));
    PrintLn(StrConcat("nach Korrektur   vy = ", f(s1.vy, 6)));
    PrintLn("");

    PrintLn("Amplitude (in gamma) | Startpunkt | erreichte Auslenkung | Umlaufzeit");
    bahn("0.05 ", 0.05, 5);
    bahn("0.30 ", 0.3,  20);
    bahn("0.40 ", 0.4,  40);
    return 0;
}

Ausgabe:

gamma L2:           0.01003712
c2:                 3.940761
lambda:             2.057073
kappa:              3.187313

Startschaetzung  vy = 0.003290
nach Korrektur   vy = 0.009846

Amplitude (in gamma) | Startpunkt | erreichte Auslenkung | Umlaufzeit
0.05  | x 1.008354 | Az 75077 km | Periode 3.1021 = 180.3 Tage
0.30  | x 1.007301 | Az 450459 km | Periode 3.0847 = 179.3 Tage
0.40  | x 1.006267 | Az 600613 km | Periode 3.0609 = 177.9 Tage