¿Altitud del satélite en función del tiempo?

Enfoque 4 (Enfoque exitoso)

En este momento, no puedo continuar con el enfoque analítico, así que encontré un satélite (escombros) a través de http://stuffin.space/?intldes=2012-008N que tenía una órbita similar a la que yo estoy tratando de determinar, escribí un pequeño script de Excel que rastrea los datos de http://www.satview.org/forec.php?sat_id=43273U una vez cada minuto, para que mañana pueda ordenar los datos y obtener la altitud aproximada como una función del tiempo.

El guion es:

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

Después de dos ejecuciones, se traza el siguiente gráfico para la órbita equivalente:

La inspección visual produce una función de órbita clara, después de eliminar los valores atípicos. Se realizan las siguientes comprobaciones de sanidad:

1. la velocidad en el perigeo (la altitud más baja), que debería ser la más alta, lo que significa que el gradiente de la función de altitud de la órbita debería ser el más alto.
2. La velocidad en el apogeo (altitud más alta) debe ser la más baja, lo que significa que el gradiente de la altitud de la órbita debe ser más bajo cerca del perigeo.
3. El gradiente de la altitud de la órbita debe ser igual a 0 (una línea horizontal en el apogeo).
4. El perfil de velocidad se traza por separado y se comprueba la continuidad.
5. El perfil de velocidad se traza por separado y se comprueba visualmente al coincidir con el gradiente de altitud de la órbita.

Las 5 verificaciones se completaron y fueron exitosas, lo que significa que no se encontraron errores en este enfoque.

Dado que la órbita tiene un eje de simetría, solo se necesita la mitad del perfil de velocidad. Como se puede ver, se realizaron dos registros, el registro de la derecha contiene una media órbita simétrica completa, por lo que produce la función de altitud de órbita completa después de una breve manipulación de datos. Estos son los datos de la órbita que representarán la altitud de la órbita en función del tiempo, para fines de estimación de la radiación.

A continuación, los datos se filtraron manualmente (marcando con colores la magnitud de las diferencias entre los puntos de datos de altitud y eliminando las diferencias más grandes hasta que estuve satisfecho con las diferencias), dando como resultado:

Después de filtrar, utilicé el siguiente software/hoja de Excel para hacer un ajuste de mínimos cuadrados: http://www.jkp-ads.com/articles/leastsquares.asp

Se utilizó la siguiente función polinomial para aproximar la función**:

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

donde x representaba el tiempo en minutos e y se calculó y comparó con$y_{actual}$. La diferencia entre las dos y se elevó al cuadrado y se sumó para todos los puntos de datos (tiempos y altitudes) raspados por el primer excel). Las condiciones de contorno (bc) se determinaron inicialmente para el algoritmo evolutivo, seleccionando el valor máximo para el cual se podía obtener la altitud más alta simplemente multiplicando ese coeficiente 1 por su x, y luego multiplicando ese valor por 10. Esto se hizo para garantizar la solución se podía obtener, pero limitaba el espacio de opciones.

Produjo la siguiente función**:

$$\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)$se calcula en$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}$$

La secuencia de comandos de Excel se puede mejorar aún más agregando un controlador de errores cerrando el ienavegador y reiniciando la secuencia de comandos para garantizar un mayor nivel de automatización, con el riesgo de mayores errores de datos.

Discusión del error: como se señaló en los comentarios, el ajuste polinomial de cuarto grado limitará la precisión del ajuste de la órbita. (Esto es visible si observa el inicio de la medición de la órbita real, se curva hacia/alrededor de -t cuando baja, mientras que el polinomio de ajuste solo parece tener 1 radio de curvatura (hacia/alrededor de +t). Esto hizo que el ajuste un poco desafiante, ya que se requirieron varias iteraciones antes de encontrar un ajuste decente. Sin embargo, este grado de precisión se consideró suficiente, por medio de una inspección visual. El procedimiento se puede mejorar mediante un cálculo del error real al cuadrado.

**Creo que un mejor enfoque sería aproximar la altitud de la órbita elíptica en función del tiempo con una sinusoidal.

Enfoque 2 (Enfoque exitoso)

La órbita se describe en la siguiente figura:

Después de una derivación larga, tediosa pero maravillosa, las siguientes 3 ecuaciones son aplicables a la órbita elíptica:

$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}}$con$\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)$= Anomalía media (=Ángulo entre la línea verde, el centro de la elipse, la línea discontinua entre el centro de la elipse y el perigeo)

$E(t)$= Anomalía excéntrica (=Ángulo entre línea roja, centro de elipse, línea entre centro de elipse y perigeo)

$\upsilon(t) = \theta(t)$= Anomalía verdadera (=Ángulo entre la línea negra, el centro de la elipse, la línea entre el centro de la elipse y el perigeo)

$a$= eje semi-mayor

$r(t)$= altitud del satélite

$\epsilon$= excentricidad

Las anomalías se visualizan en la siguiente animación de órbita. Tenga en cuenta que, como se señaló en los comentarios, M(t) es constante y desigual a E(t):

Además, es necesario calcular el período orbital$T$, semieje mayor$a$y excentricidad de la órbita$\epsilon$una vez.

Dado que la altura del perigeo = 320 km y la altura del apogeo 32190 km, el semieje mayor se encuentra usando$r_e = 6371 km$:

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

Ahora la excentricidad se encuentra usando$\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$

El periodo orbital se encuentra con$\mu_{earth} = 398600.4415 \frac{{km}^3}{s^2}$:

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

Entonces, para determinar la órbita en un punto en el tiempo t (por ejemplo, 0.5$s$), se realizan los siguientes pasos de cálculo para cada punto t en el que se desea conocer la altitud del satélite:

Se calcula la anomalía media M(t).$M= \frac{t}{T}2 \pi$

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

Entonces M se sustituye en:

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

Esto deja 1 incógnita: la anomalía excéntrica$E$. Esto se encuentra con un proceso iterativo (simplemente asuma una primera E, luego vea lo que calcula para$M$, compruebe si es igual al cálculo real$M$(probablemente no), entonces baje o aumente su estimación para$E$, y mira si te alejas más de$M$, o acercarse a él, hasta que esté satisfecho con el grado de precisión de$M$, (Lo que implica el grado de precisión de$E$).

Una vez determinada E(t) con suficiente precisión, se sustituye en:

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

Esa es la altitud de la órbita en el tiempo t = 0,5 s. El procedimiento se repite por otro tiempo t volviendo a calcular $M(t).

Además de la altitud, es posible que también desee conocer la ubicación del satélite, en relación con la Tierra, conociendo la verdadera anomalía. Para determinar la verdadera anomalía en el tiempo t, se sustituye la anomalía encontrada iterativamente$E(t)$en la siguiente fórmula:

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

Esta es la documentación del procedimiento explicado en 28:55 a 30:16 de:

.

(No escribí un guión para este {bucle iterativo de determinar/adivinar$E$} (todavía), por lo que no está semiautomatizado, a diferencia de la respuesta del enfoque 4).