Satellitenhöhe als Funktion der Zeit?

Ansatz 4 (Erfolgreicher Ansatz)

Da ich zum jetzigen Zeitpunkt den analytischen Ansatz nicht weiter verfolgen kann, habe ich über http://stuffin.space/?intldes=2012-008N einen Satelliten (Debris) gefunden , der eine ähnliche Umlaufbahn hatte wie ich versuche herauszufinden, habe ein kleines Excel-Skript geschrieben, das die Daten von http://www.satview.org/forec.php?sat_id=43273U jede Minute durchsucht, damit ich morgen die Daten sortieren und die ungefähre Höhe erhalten kann als Funktion der Zeit.

Das Skript lautet:

Sub absorb_data_from_vdb()
Application.Wait (Now + TimeValue("0:00:10")) 'H:MM:SS
'source: https://plus.google.com/105053415343007690451/posts/1pjy2PFVwN5
Dim ie As Object
Dim objelement As Object
Dim c As Integer
Dim i
Dim strCaptcha As String
Dim strImageSource As String
Dim img As HTMLhtmlElement 'tools =>reference Microsoft Html Object Library
Set ie = CreateObject("InternetExplorer.Application")

Dim click_changed_additions  As Integer
Dim electric_bullet_count As Integer
Dim gas_bullet_count As Integer


With ie
    .Visible = True
    .navigate "http://www.satview.org/forec.php?sat_id=43273U"
    'wait until first page loads
    Do Until .readyState = 4
        DoEvents
    Loop
    Application.Wait (Now + TimeValue("0:00:02")) 'H:MM:SS
    whole_day = 2
    Do While whole_day < 10
        .Refresh
        Do Until .readyState = 4
            DoEvents
        Loop
        Application.Wait (Now + TimeValue("0:00:02")) 'H:MM:SS

        Set elements0e = ie.document.getElementsByClassName("link_mudar") 'Get innertext to get value
        If elements0e.Length <> 0 Then
            'MsgBox (elements0e.item(2).innerText)
            elements0e.item(2).Focus
            elements0e.item(2).Click
            'MsgBox ("clicked link")
        End If

        If ThisWorkbook.Worksheets("Sheet1").Range("A1").Value - 2 < 2 Then
            ThisWorkbook.Worksheets("Sheet1").Range("A1").Value = 2
        End If
        For Iteration = ThisWorkbook.Worksheets("Sheet1").Range("A1").Value - 2 To ThisWorkbook.Worksheets("Sheet1").Range("A1").Value - 2 + 100
            Set elements0e = ie.document.getElementsByClassName("texto_track2") 'Get innertext to get value
            If elements0e.Length <> 0 Then
                ThisWorkbook.Worksheets("Sheet1").Range("A" & 2 + Iteration).Value = elements0e.item(0).innerHTML
                ThisWorkbook.Worksheets("Sheet1").Range("F" & 2 + Iteration).Value = Now()
                'MsgBox ("found from")
            End If
        Next Iteration
        For Iteration = ThisWorkbook.Worksheets("Sheet1").Range("A1").Value - 2 To ThisWorkbook.Worksheets("Sheet1").Range("A1").Value - 2 + 100
            Set elements0e = ie.document.getElementsByClassName("texto_track2") 'Get innertext to get value
            If elements0e.Length <> 0 Then
                ThisWorkbook.Worksheets("Sheet1").Range("G" & 2 + Iteration).Value = elements0e.item(1).innerHTML
                'MsgBox ("found from")
            End If
            If ThisWorkbook.Worksheets("Sheet1").Range("A1").Value < Iteration + 2 Then
                ThisWorkbook.Worksheets("Sheet1").Range("A1").Value = Iteration + 2
            End If

        Next Iteration

        proceed_boolean = False
        Do While proceed_boolean = False
        If Left(Right(Now(), 5), 2) <> Left(Right(ThisWorkbook.Worksheets("Sheet1").Range("F" & 2 + Iteration - 1).Value, 5), 2) Then
            'MsgBox ("A minute has passed")
            'MsgBox (Left(Right(Now(), 5), 2))
            'MsgBox (Left(Right(ThisWorkbook.Worksheets("Sheet1").Range("F" & 2 + Iteration - 1).Value, 5), 2))
            proceed_boolean = True
            Application.Wait Now + #12:00:02 AM# 'This waits for 5 seconds
        End If
        Loop


    Loop
End With
MsgBox ("Done")

End Sub

Nach zwei Läufen wird die folgende Grafik für die äquivalente Umlaufbahn gezeichnet:

Die visuelle Inspektion ergibt eine klare Bahnfunktion, nachdem die Ausreißer entfernt wurden. Folgende Plausibilitätsprüfungen werden durchgeführt:

1. die Geschwindigkeit am Perigäum (niedrigste Höhe), die am höchsten sein sollte, was bedeutet, dass der Gradient der Bahnhöhenfunktion am höchsten sein sollte.
2. Die Geschwindigkeit am Apogäum (höchste Höhe) sollte am niedrigsten sein, was bedeutet, dass der Gradient der Bahnhöhe nahe dem Perigäum am niedrigsten sein sollte.
3. Der Gradient der Umlaufbahnhöhe sollte gleich 0 sein (eine horizontale Linie am Apogäum).
4. Das Geschwindigkeitsprofil wird separat gezeichnet und auf Kontinuität geprüft.
5. Das Geschwindigkeitsprofil wird separat aufgezeichnet und visuell auf Übereinstimmung mit dem Bahnhöhengradienten überprüft.

Alle 5 Überprüfungen sind abgeschlossen und erfolgreich, was bedeutet, dass bei diesem Ansatz keine Fehler gefunden wurden.

Da die Umlaufbahn eine Symmetrieachse hat, wird nur das halbe Geschwindigkeitsprofil benötigt. Wie man sieht, wurden zwei Aufnahmen gemacht, die rechte Aufnahme enthält eine komplette halbsymmetrische Umlaufbahn, ergibt also nach kurzer Datenmanipulation die komplette Bahnhöhenfunktion. Dies sind die Umlaufbahndaten, die die Umlaufbahnhöhe als Funktion der Zeit für die Zwecke der Strahlungsschätzung darstellen.

Als nächstes wurden die Daten manuell gefiltert (indem die Größe der Unterschiede zwischen den Höhendatenpunkten farblich markiert und die größten Unterschiede entfernt wurden, bis ich mit den Unterschieden zufrieden war).

Nach dem Filtern habe ich die folgende Software/Excel-Tabelle verwendet, um die kleinsten Quadrate anzupassen: http://www.jkp-ads.com/articles/leastsquares.asp

Zur Annäherung an die Funktion** wurde folgende Polynomfunktion verwendet:

$y=ax^4+bx^3+c(x-f)^2+dx+e$

wobei x die Zeit in Minuten darstellt und y berechnet und verglichen wird$y_{actual}$. Die Differenz zwischen den beiden y's wurde quadriert und für alle Datenpunkte (Zeiten und Höhen) summiert, die vom ersten Excel geschabt wurden). Die Randbedingungen (bc) wurden zunächst für den Evolutionären Algorithmus bestimmt, indem der Maximalwert, für den die höchste Höhe erhalten werden konnte, nur durch diesen 1 Koeffizienten multipliziert mit seinem x ausgewählt und dieser Wert dann mit 10 multipliziert wurde. Dies wurde getan, um sicherzustellen Die Lösung war verfügbar, begrenzte jedoch den Optionsraum.

Es ergab die folgende Funktion**:

$$\begin{equation} \begin{split} const_a = 4.52765E-06\\ const_b = -0.002389915\\ const_c = -0.074298946\\ const_d = 252.5709003\\ const_e = 260.8620904\\ const_f = -108.1781644 \end{split} \end{equation}$$

$h(t)$eingerechnet wird$km$

$$\begin{equation} h=Const_a*xValues^4+Const_b*xValues^3+Const_c*(xValues-Const_f)^2+Const_d*xValues+Const_e \label{altitude_orbit_as_func_time} \end{equation}$$

Das Excel-Skript kann weiter verbessert werden, indem ein Errorhandler hinzugefügt wird, der den ieBrowser schließt und das Skript neu startet, um einen höheren Automatisierungsgrad zu gewährleisten, wobei das Risiko größerer Datenfehler besteht.

Fehlerdiskussion: Wie in den Kommentaren darauf hingewiesen wurde, wird die Polynomanpassung 4. Grades die Genauigkeit der Umlaufbahnanpassung einschränken. (Dies ist sichtbar, wenn Sie sich den Beginn der tatsächlichen Umlaufbahnmessung ansehen, es krümmt sich in Richtung / um -t, wenn es nach unten geht, während das passende Polynom nur einen Krümmungsradius zu haben scheint (in Richtung / um + t). Dies hat die Anpassung gemacht ein wenig herausfordernd, da mehrere Iterationen erforderlich waren, bevor eine anständige Anpassung gefunden wurde. Dieser Grad an Genauigkeit wurde jedoch mittels visueller Inspektion als ausreichend angesehen. Das Verfahren kann durch eine Berechnung des tatsächlichen Fehlerquadrats verbessert werden.

** Ich denke, ein besserer Ansatz wäre tatsächlich, die elliptische Umlaufbahnhöhe als Funktion der Zeit mit einer Sinuskurve anzunähern.

Ansatz 2 (Erfolgreicher Ansatz)

Die Umlaufbahn ist in der folgenden Abbildung beschrieben:

Nach einer langen, mühsamen, aber wunderbaren Herleitung sind die folgenden 3 Gleichungen auf die elliptische Umlaufbahn anwendbar:

$M= \frac{t}{T}2 \pi$

$M = E-\epsilon sin{E}$

$r = a(1-\epsilon cos {E})$

$\upsilon = \theta = 2 arctan (\sqrt{\frac{1+\epsilon}{1-\epsilon}} tan{\frac{E}{2}})$

$T=2\pi\sqrt{\frac{a^3}{\mu}}$mit$\mu_{earth} = 398600.4415 \frac{{km}^3}{s^2}$

$r = \frac{a(1-e^2)}{1+e cos{\theta}} = a(1-\epsilon cos{E})$

$M(t)$= Mittlere Anomalie (=Winkel zwischen grüner Linie, Ellipsenmittelpunkt, gestrichelte Linie zwischen Ellipsenmittelpunkt und Perigäum)

$E(t)$= Exzentrische Anomalie (=Winkel zwischen roter Linie, Ellipsenmittelpunkt, Linie zwischen Ellipsenmittelpunkt und Perigäum)

$\upsilon(t) = \theta(t)$= Wahre Anomalie (=Winkel zwischen schwarzer Linie, Ellipsenmittelpunkt, Linie zwischen Ellipsenmittelpunkt und Perigäum)

$a$= große Halbachse

$r(t)$= Höhe des Satelliten

$\epsilon$= Exzentrizität

Die Anomalien werden in der folgenden Orbit-Animation visualisiert. Beachten Sie, dass, wie in den Kommentaren erwähnt, M(t) konstant und ungleich E(t) ist:

Außerdem muss die Umlaufzeit berechnet werden$T$, große Halbachse$a$und Bahnexzentrizität$\epsilon$einmal.

Da die Perigäumshöhe = 320 km und die Apogäumshöhe 32190 km beträgt, wird die große Halbachse durch Verwendung gefunden$r_e = 6371 km$:

$a = \frac{320+r_e+32190+r_e}{2}= 25971.5 km$

Jetzt wird die Exzentrizität mit ermittelt$\theta = 0 rad$:

$r = \frac{a(1-e^2)}{1+e cos{\theta}}$

$r = \frac{a(1-e^2)}{1+e} = a(1-e)$

$e =1-\frac{r}{a} = 1-\frac{r_e+320}{a} =1-\frac{6371+320}{25971.5}=0.74237144562308$

Die Umlaufzeit wird mit gefunden$\mu_{earth} = 398600.4415 \frac{{km}^3}{s^2}$:

$T=2\pi\sqrt{\frac{a^3}{\mu}}=41685.4 s$

Um also die Umlaufbahn zu einem Zeitpunkt t (zB 0,5$s$) werden für jeden Punkt t, an dem wir die Satellitenhöhe wissen wollen, folgende Rechenschritte durchgeführt:

Die mittlere Anomalie M(t) wird berechnet.$M= \frac{t}{T}2 \pi$

$M= \frac{0.5}{41685.4}2 \pi=0.0000753643.. rad$

Dann wird M eingesetzt in:

$M = E-\epsilon sin{E}$

Damit bleibt eine Unbekannte: die exzentrische Anomalie$E$. Dies wird mit einem iterativen Prozess gefunden (Nehmen Sie einfach ein erstes E an und sehen Sie dann, wofür Sie berechnen$M$, überprüfen Sie, ob es dem tatsächlich berechneten entspricht$M$(wahrscheinlich nicht), also senken oder erhöhen Sie Ihre Schätzung für$E$, und sehen Sie, ob Sie sich weiter davon entfernen$M$, oder nähern Sie sich ihm, bis Sie mit dem Genauigkeitsgrad von zufrieden sind$M$, (Was den Grad der Genauigkeit von impliziert$E$).

Sobald E(t) mit ausreichender Genauigkeit bestimmt ist, wird es ersetzt durch:

$r = a(1-\epsilon cos{E})$finden$r{t}$.

Das ist die Bahnhöhe zum Zeitpunkt t = 0,5 s. Der Vorgang wird für eine weitere Zeit t wiederholt, indem $M(t) erneut berechnet wird.

Zusätzlich zur Höhe möchte man vielleicht auch die Position des Satelliten relativ zur Erde wissen, indem man die wahre Anomalie kennt. Um die wahre Anomalie zum Zeitpunkt t zu bestimmen, ersetzt man die iterativ gefundene$E(t)$in folgender Formel:

$\upsilon = \theta = 2 arctan (\sqrt{\frac{1+\epsilon}{1-\epsilon}} tan{\frac{E}{2}})$

Dies ist die Dokumentation des in 28:55 bis 30:16 erläuterten Verfahrens von:

.

(Ich habe kein Skript für diese {iterative Schleife des Bestimmens/Ratens geschrieben$E$} (noch) also nicht halbautomatisiert, im Gegensatz zur Antwort von Ansatz 4.)