Brouwer-Lyddane-Mean-Semi-Major-Axis-Bias

Ich ergänze meinen Kommentar :

Ich bin mit der Software nicht vertraut, aber wenn Sie den Satelliten mit Ausnahme der Zeit (12 Stunden Unterschied) zum gleichen Zeitpunkt starten, startet der Satellit in beiden Fällen mit unterschiedlichen Gravitationspotentialen. Ich vermute, wenn Sie es 100 Mal mit gleichmäßigen Abständen zwischen Ihren Startepochen über 24 Stunden ausführen würden, würden Sie sehen, dass sich ihre "Höhenverzerrungen" nach oben und unten verteilen würden, je nachdem, ob Sie in einem Gravitationsgipfel oder -tal begonnen haben. Es sollte nicht erwartet werden, dass der Mittelwert gleich dem momentanen Startwert sein sollte, Sie beginnen zufällig.

Um es einfach zu halten, verwende ich nur die$J_2$Multipolkomponente des Gravitationspotentials der Erde zusätzlich zum Monopolterm$GM_E$. Ich habe ein kurzes Python-Skript geschrieben, das Zahlen und Gleichungen aus Wikipedia verwendet (Links werden im Skript angezeigt).

Ich habe eine um 90 Grad geneigte Umlaufbahn in einer Höhe von etwa 901 Kilometern konstruiert und für 14 Umlaufbahnen oder etwa 1 Tag propagiert. (Ich habe diese Höhe gewählt, damit ich die Neigung auf 99 Grad einstellen und bestätigen konnte, dass sie sonnensynchron war, indem ich 91 Tage lang propagierte und beobachtete, wie sich der aufsteigende Knoten von der x-Achse zur y-Achse bewegte, wodurch überprüft wurde, dass der Gradient von das$J_2$potenzielle Laufzeit ist in Ordnung.)

Ich habe eine große Halbachse und die Umlaufgeschwindigkeit berechnet$a_0$und$v_0$für eine Kreisbahn im Monopolfeld$GM_E/r$und diese für Startbedingungen verwendet.

Die erste Handlung zeigt$x(t), y(t), z(t)$und$r(t)-a_0$für 24 Stunden (14 Umläufe), beginnend am Äquator mit phase = 0. Die Berechnung wurde für Grad wiederholt phase = 30, 60, 90, so dass der letzte über dem Pol begann, ebenfalls in der gleichen Höhe von 901 Kilometern.

Im zweiten Plot zeige ich die vier Plots von$r(t)-a_0$die auf der gleichen Höhe beginnen, aber an unterschiedlichen Orten über der Erde. Sie können sofort sehen, dass die "Vorspannung", von der Sie sprechen, vom Startort und dem dort vorhandenen besonderen Gravitationspotential abhängt.

Beginnend über dem Äquator mit phase = 0Grad, wo man tiefer in der Schwerkraft der Erde sitzt, ist die Geschwindigkeit für eine kreisförmige Umlaufbahn etwas niedrig, und so befindet sich der Satellit tatsächlich in der Apoapsis und fällt eine halbe Periode später auf eine etwas niedrigere Höhe ab.

Beginnt man jedoch über dem Pol mit phase = 90Graden, die höher in der Schwerkraft der Erde sind, ist die Geschwindigkeit für eine kreisförmige Umlaufbahn etwas hoch, und daher ist dies jetzt Periapsis, und die Höhe steigt eine halbe Periode später an.

Zusammenfassung: Durch die Wahl eines leichter verständlichen Gravitationsfeldes ist das Verhalten der „Bias“ hier systematischer. In Ihrer Frage verwenden Sie ein sehr komplexes, uneben aussehendes Gravitationsfeld, das sich ständig dreht. Indem Sie am selben Ort im Weltraum, aber in verschiedenen Epochen, begannen, drehten Sie tatsächlich die Erde und bewegten unterschiedliche Gravitationspotentiale zum Startort Ihres Satelliten .

Ich sage jetzt voraus , dass, wenn Sie das Multipolfeld in JGM-2 ausschalten und ein kugelsymmetrisches Potential verwenden, Ihre Berechnungen eine flache Höhe, keine Variationen und keine Verzerrungen zeigen werden.

Schlussgedanke: Realistische Bahnen sind keine perfekten Kegelschnitte, und daher repräsentieren sie und ihre Keplerschen Elemente keine realistischen Bahnen . Sie sind nur Annäherungen an die Realität und stimmen daher nicht, obwohl sie nah dran sind.

Keplerische Elemente wurden verwendet, als Menschen mit Federn schrieben, indem sie Licht aus der Verbrennung von Tierfett verwendeten (wenn sie nicht damit beschäftigt waren, selbst auf dem Scheiterhaufen verbrannt zu werden). Sie sind ein gemischter Segen im 21. Jahrhundert, wo alles so viel mehr Ziffern hat.

def deriv(X, t):

    rr, vel = X.reshape(2, -1)

    acc0, acc2 = accs(rr)

    return np.hstack([vel, acc0 + acc2])


def accs(rr):

    x,   y,   z   = rr
    xsq, ysq, zsq = rr**2

    rsq = (rr**2).sum()
    rm3 = rsq**-1.5
    rm7 = rsq**-3.5

    acc0  = -GM_earth * rr * rm3

    # https://en.wikipedia.org/wiki/Geopotential_model#The_deviations_of_Earth.27s_gravitational_field_from_that_of_a_homogeneous_sphere
    acc2x = x * rm7 * (6*zsq - 1.5*(xsq + ysq))
    acc2y = y * rm7 * (6*zsq - 1.5*(xsq + ysq))
    acc2z = z * rm7 * (3*zsq - 4.5*(xsq + ysq))

    acc2  = J2_earth * np.hstack((acc2x, acc2y, acc2z))

    return acc0, acc2


import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

halfpi, pi, twopi = [f*np.pi for f in [0.5, 1, 2]]
degs, rads        = 180./pi, pi/180.

GM_earth = 3.986004418E+14 #  m^3/s^2  https://en.wikipedia.org/wiki/Standard_gravitational_parameter
J2_earth = 1.7555E+25      #  m^5/s^2  https://en.wikipedia.org/wiki/Geopotential_model

a0  = (901 + 6378) * 1E+03 # meters (roughly) for a n=14 sun-synchronous orbit

v0  = np.sqrt(GM_earth/a0)  # m/s  (vis-viva)
T   = twopi * a0 / v0       # sec

print T/60, ' minutes'

phases = [rads*d for d in [0, 30, 60, 90]]

X0s = []
for phase in phases:
    sp, cp = np.sin(phase), np.cos(phase)

    X0 = np.hstack(( [ a0*cp, 0, a0*sp],   # initial positions
                     [-v0*sp, 0, v0*cp] )) # initial velocitys
    X0s.append(X0)


n, days = 14, 1.  # orbits/day, days

time = np.linspace(0, n*days*T, int(n*days*100)+1)  

answers = []
for X0 in X0s:

    answer, info = ODEint(deriv, X0, time, full_output=True)

    answers.append(answer)

if 1 == 1:

    names  = 'x(t)', 'y(t)', 'z(t)', 'r(t)-a0'
    answer = answers[1]   # just choose one for the plot
    posn   = answer.T[:3]
    rdiff  = np.sqrt((posn**2).sum(axis=0)) - a0
    things = [thing for thing in posn] + [rdiff]

    plt.figure()
    for i, (thing, name) in enumerate(zip(things, names)):
        plt.subplot(4, 1, i+1)
        plt.plot(time/(24.*3600), thing)   # x-axis in days, not seconds
        plt.title(name, fontsize=16)
        plt.xlim(0, days+0.01)  # a little more than n days
    plt.show()

if 1 == 1:

    plt.figure()

    for answer in answers:

        posn   = answer.T[:3]
        rdif   = np.sqrt((posn**2).sum(axis=0)) - a0

        plt.plot(time/(24.*3600), rdif)   # x-axis in days, not seconds

    plt.xlim(0, days+0.1)  # a little more than n days
    plt.text(1.01, -6000, u' 0\u00B0', fontsize=15 )
    plt.text(1.01,  2000, u'30\u00B0', fontsize=15 )
    plt.text(1.01, 16000, u'60\u00B0', fontsize=15 )
    plt.text(1.01, 22000, u'90\u00B0', fontsize=15 )
    plt.title('r(t)-a0', fontsize=16)
    plt.show()

Unten: eine schnelle Überprüfung, die 91 Tage lang mit einer Neigung von 99 Grad läuft, um zu bestätigen, dass die Umlaufbahn ungefähr sonnensynchron ist, um sicherzustellen, dass es zumindest keine groben Fehler bei der Art gibt, wie ich die Gleichungen für die Steigung der eingetippt habe$J_2$Komponente oder die numerischen Werte. Das obige Skript lässt keine Neigung zu, aber Sie würden es durch Mischen von einbeziehen$v_y$und$v_z$Werte, wenn X0berechnet wird.