04.02.2021 08:28
K-ILYA-V 8 1400 Источник

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

Алгоритм подразумевает замену непрерывного малого ускорения на серию мгновенных импульсов прикладываемых к КА через равные промежутки времени, при этом промежутки времени подбираются автоматически таким образом чтобы при разгоне КА никогда не покидал "условный" 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) Если условие указанное ниже верно, выходим из алгоритма:

\Delta V > V_e

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

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