Высота спутника как функция времени?

Подход 4 (Успешный подход)

На данный момент времени я не могу дальше использовать аналитический подход, поэтому я нашел спутник (обломки) через http://stuffin.space/?intldes=2012-008N , который имел орбиту, похожую на ту, которую я я пытаюсь определить, написал небольшой сценарий Excel, который сканирует данные с http://www.satview.org/forec.php?sat_id=43273U каждую минуту, чтобы завтра я мог отсортировать данные и получить приблизительное значение высоты. как функция времени.

Сценарий:

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

После двух прогонов для эквивалентной орбиты строится следующий график:

Визуальный осмотр дает четкую функцию орбиты после удаления выбросов. Выполняются следующие проверки работоспособности:

1. скорость в перигее (наименьшая высота), которая должна быть самой высокой, что означает, что градиент функции высоты орбиты должен быть самым высоким.
2. Скорость в апогее (наибольшая высота) должна быть самой низкой, что означает, что градиент высоты орбиты должен быть самым низким вблизи перигея.
3. Градиент высоты орбиты должен быть равен 0 (горизонтальная линия в апогее).
4. Профиль скорости строится отдельно и проверяется на непрерывность.
5. Профиль скорости строится отдельно и визуально проверяется на совпадение с градиентом высоты орбиты.

Все 5 проверок завершены и успешны, что означает, что в этом подходе не было обнаружено ошибок.

Поскольку орбита имеет ось симметрии, требуется только половина профиля скорости. Как видно, были сделаны две записи, правая запись содержит полную полусимметричную орбиту, поэтому после непродолжительной обработки данных она дает полную функцию высоты орбиты. Это данные об орбите, которые будут представлять высоту орбиты как функцию времени для целей оценки излучения.

Затем данные были отфильтрованы вручную (помечая цветом величину различий между точками данных о высоте и удаляя самые большие различия, пока я не был удовлетворен различиями). Получилось:

После фильтрации я использовал следующее программное обеспечение/таблицу Excel, чтобы подогнать метод наименьших квадратов: http://www.jkp-ads.com/articles/leastsquares.asp

Для приближения к функции** использовалась следующая полиномиальная функция:

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

где x представляет время в минутах, а y вычисляется и сравнивается с$y_{actual}$. Разница между двумя y была возведена в квадрат и суммирована для всех точек данных (время и высота), считанных первым Excel). Граничные условия (bc) были первоначально определены для эволюционного алгоритма путем выбора максимального значения, для которого наибольшая высота может быть получена только с помощью этого 1 коэффициента, умноженного на его x, а затем умножения этого значения на 10. Это было сделано для обеспечения решение было доступно, но ограничивало пространство вариантов.

Получилась следующая функция**:

$$\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)$вычисляется в$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}$$

Сценарий Excel можно дополнительно улучшить, добавив обработчик ошибок, закрывающий ieбраузер и перезапускающий сценарий, чтобы обеспечить более высокий уровень автоматизации с риском большего количества ошибок данных.

Обсуждение ошибки: как было указано в комментариях, полиномиальная подгонка 4-й степени ограничит точность подгонки орбиты. (Это видно, если вы посмотрите на начало фактического измерения орбиты, оно изгибается в направлении/около -t при движении вниз, тогда как подгоночный полином имеет только 1 радиус изгиба (в направлении/вокруг +t). Это сделало подгонку немного сложно, поскольку потребовалось несколько итераций, прежде чем была найдена подходящая подгонка.Однако эта степень точности была сочтена достаточной с помощью визуального осмотра.Процедура может быть улучшена путем вычисления фактического квадрата ошибки.

** Я думаю, что лучшим подходом действительно было бы аппроксимировать высоту эллиптической орбиты как функцию времени синусоидой.

Подход 2 (Успешный подход)

Орбита описана на следующем рисунке:

После долгого, утомительного, но замечательного вывода следующие 3 уравнения применимы к эллиптической орбите:

$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}}$с$\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)$= Средняя аномалия (= Угол между зеленой линией, центром эллипса, пунктирной линией между центром эллипса и перигеем)

$E(t)$= Эксцентрическая аномалия (= Угол между красной линией, центром эллипса, линией между центром эллипса и перигеем)

$\upsilon(t) = \theta(t)$= Истинная аномалия (= Угол между черной линией, центром эллипса, линией между центром эллипса и перигеем)

$a$= большая полуось

$r(t)$= высота спутника

$\epsilon$= эксцентриситет

Аномалии визуализируются на следующей анимации орбиты. Обратите внимание, что, как было указано в комментариях, M(t) постоянна и не равна E(t):

Кроме того, необходимо вычислить орбитальный период$T$, большая полуось$a$и эксцентриситет орбиты$\epsilon$однажды.

Учитывая, что высота перигея = 320 км, а высота апогея 32190 км, большая полуось находится с помощью$r_e = 6371 km$:

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

Теперь эксцентриситет находится с помощью$\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$

Орбитальный период находится с$\mu_{earth} = 398600.4415 \frac{{km}^3}{s^2}$:

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

Таким образом, чтобы определить орбиту в момент времени t (например, 0,5$s$), для каждой точки t, в которой мы хотим узнать высоту спутника, выполняются следующие шаги вычисления:

Вычисляется средняя аномалия M(t).$M= \frac{t}{T}2 \pi$

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

Затем М заменяется на:

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

Остается одно неизвестное: эксцентрическая аномалия.$E$. Это поиск с помощью итеративного процесса (просто предположим, что сначала E, а затем посмотрим, что вы вычислите для$M$, проверьте, равно ли оно фактическому вычисляемому$M$(вероятно, нет), поэтому уменьшите или увеличьте свою оценку для$E$, и посмотрите, не отклонитесь ли вы дальше от$M$, или приближайтесь к нему, пока не будете удовлетворены степенью точности$M$, (От чего зависит степень точности$E$).

Как только E(t) определено с достаточной степенью точности, оно подставляется в:

$r = a(1-\epsilon cos{E})$найти$r{t}$.

Это высота орбиты в момент времени t = 0,5 с. Процедура повторяется еще раз t путем повторного вычисления $M(t).

Помимо высоты, можно также захотеть узнать местоположение спутника относительно Земли, зная истинную аномалию. Чтобы определить истинную аномалию в момент времени t, подставляют найденное итеративно$E(t)$в следующей формуле:

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

Это документация по процедуре, описанной с 28:55 по 30:16:

.

(Я не писал скрипт для этого {итеративного цикла определения/угадывания$E$} (пока), поэтому он не полуавтоматический, в отличие от ответа подхода 4.)