Алгоритм симуляция движения космического аппарата по спиральной траектории под действием малого ускорения

от автора

Рассмотрим алгоритм симуляция движения космического аппарата (КА) по спиральной траектории (СТ) под действием малого ускорения (МУ).

Алгоритм подразумевает замену непрерывного малого ускорения на серию мгновенных импульсов прикладываемых к КА через равные промежутки времени, при этом промежутки времени подбираются автоматически таким образом чтобы при разгоне КА никогда не покидал «условный» I квадранта координатной плоскости, что наиболее точно соответствует движению с непрерывным ускорением.

Алгоритм будет рассмотрен аналитически без привязки к какому либо компьютерному языку или конкретной планете для возможности его последующей адаптации к любому компьютерному языку и/или планете по выбору.

Алгоритм был проверен автором на MASM64 с применением  Extended precision (расширенной точностью).

  • Параметры нулевой точки

Зная экваториальный радиус (Re) планеты и высоту (Halt) положения КА над экватором, находим радиус перицентра (rp), по формуле:

r_p = R_e + H_{ alt }

1) Зная гравитационный параметр (mu) находим скорость (V0) относительно (ЦГ) планеты, по формуле:

V_0 = \sqrt { \mu / r_p}

2) Находим время (Tp) полного оборота (ПО) по формуле:

T_p = \frac{2 \pi r_p}{V_0}

3) Зная мощность (N0) двигателя и скорость (Vrm) истечения реактивной массы, находим расход реактивной массы (mr) в секунду, по формуле:

m_{r} =\frac{2 N_0}{V_{rm}^2}

4) Находим начальную продолжительность (dTI) импульса двигателя, определяя ее как 1/4 периода ПО по формуле:

\Delta T_I = T_p / 4

5) Принимаем начальное время (dT0) равным нулю, по формуле:

\Delta T_0 = 0.0

  • Точка рестарта

6) Находим расход реактивной массы (mrdt) за один импульс продолжительностью dTI, по формуле:

m_{ r \Delta T_I } = m_{ rt } \cdot \Delta T_I

7) Находим массу (mc) после импульса, по формуле:

m_c = m_0 - m_{ r \Delta T_I }

8) Находим дельту скорости (dV0) после импульса, по формуле:

\Delta V_0 = \ln ( m_0 / m_c ) \cdot V_{rm}

9) Находим скорость (V) после импульса, по формуле:

V = V_0 + \Delta V

10) Находим эксцентриситет (e) орбиты, по формуле:

e = \frac { V^2 \cdot r_p } { \mu } - 1.0

11) Находим большую полуось (a), по формуле:

a = \frac { r_p } { 1.0 - e }

  • Параметры N точки траектории до импульса

12) Находим средние движение (n), по формуле:

n = \sqrt { \mu / a^3 }

13) Находим среднею аномалию (M), по формуле:

M = n ( \Delta T_I + \Delta T_{ n-1 } )

14) Принимаем приблизительную среднюю аномалию Mapprox(0) равной среднею аномалию (M), по формуле:

M_{ approx (0) } = M

15) Находим приблизительную эксцентрическую аномалию Eapprox(n), по формуле:

E_{ approx (n) } = M_{ approx (n-1) } - \frac { \sin (M_{ approx (n-1) } ) \cdot e} { 1.0 - \cos ( M_{ approx (n-1) } ) \cdot e }

16) Находим приблизительную среднюю аномалию Mapprox(n), по формуле:

M_{ approx (n) } = E_{ approx (n) } - \sin ( E_{ approx (n) } )

17) Если условие указанное ниже верно, переходим к пункту № 15:

M \neq M_{ approx (n) }

18) Принимаем эксцентрическую аномалию (E) равной Eapprox(n), по формуле:

E = E_{ approx (n) }

19) Если условие указанное ниже верно, переходим к пункту № 21 :

E < \frac{ \pi } {2}

20) Принимаем новое значение dTI по формуле указанной ниже и переходим к пункту № 6:

\Delta T_I = \Delta T_I / 2

21) Находим катет эксцентрической аномалии (X), по формуле:

X = ( \cos (E) - e) \cdot a

22) Находим фокальный параметр (p), по формуле:

p = a \cdot (1.0 - e^2)

23) Находим катет эксцентрической аномалии (Y), по формуле:

Y = \sin (E) \cdot a \cdot \sqrt { (1.0 - e^2) }

24) Находим радиус-вектор (r), по формуле:

r = \sqrt { X^2 + Y^2}

25) Находим синус фета (sin(theta)), по формуле:

\sin (\theta) = Y / r

26) Находим радиальную скорость (Vrad), про формуле:

V_{rad} = \sqrt { \frac { e^2 \cdot \sin^2 ( \theta ) \cdot \mu } {p} }

27) Находим поперечную скорость (Vnor), по формуле:

V_{nor} = \sqrt {p / r^2}

28) Находим полную скорость (Vn), по формуле:

V_n = \sqrt { V_{rad}^2 + V_{nor}^2}

29) Находим косинус фи (cos(phi)), по формуле:

\cos ( \phi ) = V_{rad} / V_n

30) Находим синус фи (sin(phi)), по формуле:

\sin ( \phi ) = V_{nor} / V_n

  • Параметры N точки после импульса

31) Находим массу (mc) после импульс, по формуле:

m_c = m_n = m_{ (n-1) } - m_{ r \Delta T_I }

32) Находим дельту скорости (dV) после импульса, по формуле:

\Delta V = \ln ( m_{n-1} / m_n ) \cdot V_{rm}

33) Находим скорость (Ve) выхода из гравитационного поля планеты, по формуле:

V_e = \sqrt{ 2 \cdot V_n }

34) Если условие указанное ниже верно, выходим из алгоритма:

V_{n+1} = V_n + \Delta V_n

36) Находим радиальную скорость (Vrad), по формуле:

V_{rad} = \cos (\phi) \cdot V_{n+1}

37) Находим фокальный параметр (p), по формуле:

p = \frac { (V_{n+1} \cdot \sin ( \phi ) \cdot r )^2 } { \mu }

38) Находим эксцентриситет (е), по формуле:

e = \sqrt { \frac {V_{rad}^2 \cdot p} { \mu } + \bigg( \frac { p - r } { r } \bigg)^2 }

39) Находим косинус фета cos(theta), по формуле:

\cos (\theta) = \frac {p-r} {r \cdot e}

40) Находим катет эксцентрической аномалии (Y), по формуле:

Y = r \sqrt { 1 - \cos^2 ( \theta ) }

41) Находим большую полуось (а), по формуле:

a = \frac { p }{ 1.0 - e^2 }

42) Находим катет эксцентрической аномалии (X), по формуле:

X = \cos ( \theta ) \cdot r + e \cdot a

43) Находим синус фета sin(theta), по формуле:

\sin ( \theta ) = Y / r

44) Находим эксцентрическую аномалию E, по формуле:

E = \arcsin (E)

45) Находим среднею аномалию (Mn), по формуле:

M_n = E - \sin (E)

46) Находим средние движение (n), по формуле:

n = \sqrt { \frac { \mu } { a^3 } }

47) Находим среднею аномалию (Mn+1), по формуле:

M_{n+1} = M_n / n + \Delta T_I

48) Переход к пункту № 14

ссылка на оригинал статьи https://habr.com/ru/post/540840/


Комментарии

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *