Как рассчитать угол траектории полета γ по вектору состояния?

Question

Как рассчитать угол траектории полета γ по вектору состояния?

ооо

Отличный ответ @Julio описывает угол траектории полета и объясняет, что это угол между тангенциальным направлением (перпендикулярным радиальному вектору к центральному телу) и текущим вектором скорости.

Сначала я попытался получить угол из этого выражения, но это явно неправильно, так как $\arccos$ является четной функцией, и угол может идти от $-\pi/2$ к $\pi/2$ :

арккос (\frac{р \cdot в}{| р | | в |}) - \frac{π}{2} (неправильно!)

$\arccos\left(\frac{\mathbf{r \centerdot v}}{|\mathbf{r}| \ |\mathbf{v}|} \right) - \frac{\pi}{2} \ \ \ \text{ (incorrect!)}$

Я интегрировал орбиты для GM ( $\mu$ ) и СМА ( $a$ ) единицы и стартовые расстояния от 0,2 до 1,8. Это делает период всегда $2 \pi$ . Когда я рисую результат своей функции, я получаю слишком много покачиваний.

Какое выражение я могу использовать, чтобы получить правильную гамму угла траектории полета, начиная с векторов состояния?

Пересмотренный python для ошибочной части был бы оценен, но, конечно, не обязателен для ответа.

def deriv(X, t):
    x, v = X.reshape(2, -1)
    acc = -x * ((x**2).sum())**-1.5
    return np.hstack((v, acc))

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]

T    = twopi
time = np.linspace(0, twopi, 201)

a       = 1.0
rstarts = 0.2 * np.arange(1, 10)
vstarts = np.sqrt(2./rstarts - 1./a)  # from vis-viva equation

answers = []
for r, v in zip(rstarts, vstarts):
    X0 = np.array([r, 0, 0, v])
    answer, info = ODEint(deriv, X0, time, full_output= True)
    answers.append(answer.T)

gammas = []
for a in answers:
    xx, vv = a.reshape(2, 2, -1)
    dotted = ((xx*vv)**2).sum(axis=0)
    rabs, vabs = [np.sqrt((thing**2).sum(axis=0)) for thing in (xx, vv)]
    gamma = np.arccos(dotted/(rabs*vabs)) - halfpi
    gammas.append(gamma)

if True:

    plt.figure()
    plt.subplot(4, 1, 1)
    for x, y, vx, vy in answers:
        plt.plot(x, y)
        plt.plot(x[:1], y[:1], '.k')
    plt.plot([0], [0], 'ok')
    plt.title('y vs x')

    plt.subplot(4, 1, 2)
    for x, y, vx, vy in answers:
        plt.plot(time, x, '-b')
        plt.plot(time, y, '--r')
    plt.title('x (blue) y (red, dashed)')
    plt.xlim(0, twopi)

    plt.subplot(4, 1, 3)
    for x, y, vx, vy in answers:
        plt.plot(time, vx, '-b')
        plt.plot(time, vy, '--r')
    plt.title('vx (blue) vy (red), dashed')
    plt.xlim(0, twopi)

    plt.subplot(4, 1, 4)
    for gamma in gammas:
        plt.plot(time, gamma)
    plt.title('gamma?')
    plt.xlim(0, twopi)

    plt.show()

пользователь20636

должен ли этот вопрос быть TLDRed, чтобы указать, что это была ошибка кодирования, поскольку он все еще, кажется, спрашивает, что не так с формулой

Ответы (2)

Как рассчитать угол траектории полета γ по вектору состояния?

должен ли этот вопрос быть TLDRed, чтобы указать, что это была ошибка кодирования, поскольку он все еще, кажется, спрашивает, что не так с формулой

Том Спилкер · Answer 1

Это проблема, которая преследует группы людей, очень хорошо осведомленных об орбитальной динамике, но которые учились по разным учебникам: есть два разных определения «угла траектории полета»!!

В дополнение к $\gamma$ , угол между тангенциальным направлением и вектором скорости, есть $\beta$ , угол между радиальным направлением и вектором скорости. Люди часто говорят «угол траектории полета», не говоря, какое определение они используют . Сбивает с толку! (Я только что заметил, что диаграмма в ответе Хулио также показывает $\beta$ )

Если вы работаете с $\beta$ вместо того $\gamma$ , $\beta$ дан кем-то

\begin{matrix} (1) & арккос (\frac{р \cdot в}{| р | | в |}) \end{matrix}

$\arccos\left(\frac{\mathbf{r \centerdot v}}{|\mathbf{r}| \ |\mathbf{v}|} \right) \tag{1}$

который идет от 0 ("прямо вверх") до $\pi$ ("прямо вниз"). С использованием $\gamma$ , "прямо вверх" это $\pi/2$ и "прямо вниз" $-\pi/2$ , так что преобразование $\beta$ к $\gamma$ ты просто вычитаешь $\beta$ из $\pi/2$ :

\begin{matrix} (2) & γ знак равно π / 2 - арккос (\frac{р \cdot в}{| р | | в |}) \end{matrix}

$\gamma = \pi/2 - \arccos\left(\frac{\mathbf{r \centerdot v}}{|\mathbf{r}| \ |\mathbf{v}|} \right) \tag{2}$

Это эквивалентно

\begin{matrix} (3) & γ знак равно арксин (\frac{р \cdot в}{| р | | в |}) \end{matrix}

$\gamma = \arcsin\left(\frac{\mathbf{r \centerdot v}}{|\mathbf{r}| \ |\mathbf{v}|} \right) \tag{3}$

Я не знаком с языком, который вы использовали для своих расчетов и графиков, поэтому я не смотрел ваш алгоритм, чтобы понять, почему «слишком много покачиваний».

Спасибо! Я добавил теги (числа) к уравнениям. Скажете ли вы, что покачивания слишком много, или такое поведение покачивания на самом деле разумно? Поскольку ваш $\beta$ (уравнение 1) совпадает с моим ошибочным $\gamma$ за исключением смещения в половину числа пи, тогда колебания на моем графике должны быть такими же, как и на правильном графике вашего $\beta$ (уравнение 1).
Мне кажется слишком много люфтов. Я проверю это позже.
@uhoh, На самом деле, мое уравнение 1 - это просто отрицательная часть вашего уравнения. Что-то еще не так. Конечно, вы знаете, что что- то не так, потому что все задуманное $\gamma$ s отрицательны или равны нулю, чего не может быть, кроме внутренней спирали. Для кеплеровской эксцентрической орбиты $\gamma$ должен пересекать ноль ровно дважды, в перицентре и апоапсисе, и быть монотонным между экстремумами как для короткого (экстремум через перицентр до другого экстремума), так и для длинного (экстремум от апоцентра до другого экстремума) отрезков. Я посмотрю, смогу ли я привести пример того, что $\gamma$ кривая должна выглядеть.
@uhoh, О, да: действительно слишком много покачиваний. Вот что значит "монотонность". Каким-то образом $\gamma$ расчет сбился. Убедитесь, что аргумент arccos соответствует вашим намерениям.
Хорошо спасибо. Я ожидаю, что кто-то сможет понять, что на самом деле не так. Это может оказаться я, я вернусь к этому позже сегодня.
Упс, я должен был сказать выше: «Мое уравнение 2 — это просто отрицательное выражение вашего». Я должен выйти из системы и лечь спать!
Спасибо за всю твою помощь! Я нашел ошибку , а также лучше понял $\beta$ и $\gamma$ .
Я процитировал вас здесь , посмотрите и убедитесь, что я сделал это правильно, или у вас есть что добавить. Спасибо!
Теперь я официально запутался! Аварийный словарный запас! Как может касательная скорость эллиптической кеплеровской орбиты не касаться орбиты?
@uhoh «тангенциальный» к сфере с центром в первичной точке, а не к орбите. Лично я бы предпочел сказать «поперечная скорость», но мой первый профессор орбитальной динамики в Стэнфорде использовал «тангенциальную».

ооо · Answer 2

Я нашел ошибку в скрипте, это было из-за моего "доморощенного" точечного продукта. У меня был лишний квадрат:

dotted = ((xx*vv)**2).sum(axis=0)     # WRONG

dotted = (xx*vv).sum(axis=0)          # Correct

Итак, используя это плюс отличные разъяснения @TomSpilker, я использовал следующие два метода для расчета гаммы:

Способ 1:

\begin{matrix} (3) & γ_{1} знак равно арксин (\frac{р \cdot в}{| р | | в |}) \end{matrix}

$\gamma_1 = \arcsin\left(\frac{\mathbf{r \centerdot v}}{|\mathbf{r}| \ |\mathbf{v}|} \right) \tag{3}$

Способ 2:

Альтернативный метод грубой силы для двойной проверки:

θ_{р} знак равно арктический 2 (у, Икс)

$\theta_r = \arctan2(y, x)$

θ_{в} знак равно арктический 2 (в у, Икс)

$\theta_v = \arctan2(vy, x)$

θ_{т а н Дж} знак равно θ_{р} + \frac{π}{2}

$\theta_{tanj} = \theta_r + \frac{\pi}{2}$

γ_{2} знак равно θ_{т а н Дж} - θ_{в}

$\gamma_2 = \theta_{tanj} - \theta_v$

γ_{2 м о д} знак равно мод (γ_{2} + π, 2 π) - π

$\gamma_{2mod} = \mod(\gamma_2+ \pi, 2\pi) - \pi$

Операция по модулю действительно нужна только в компьютерной программе, поскольку каждая тета получается из отдельной операции arctan2:

gammas_1, gammas_2 = [], []
for a in answers:

    xx, vv = a.reshape(2, 2, -1)

    dotted = (xx*vv).sum(axis=0)
    rabs, vabs = [np.sqrt((thing**2).sum(axis=0)) for thing in (xx, vv)]
    gamma_1 = np.arcsin(dotted/(rabs*vabs))   # Per Tom Spilker's answer Eq. 3

    theta_r = np.arctan2(xx[1], xx[0])
    theta_v = np.arctan2(vv[1], vv[0])
    theta_tanj = theta_r + halfpi

    gamma_2 = theta_tanj - theta_v
    gamma_2 = np.mod(gamma_2 + pi, twopi) - pi

    gammas_1.append(gamma_1)
    gammas_2.append(gamma_2)

plt.figure()
plt.subplot(2, 1, 1)
for gamma_1 in gammas_1:
    plt.plot(time, gamma_1)
plt.title('gammas_1', fontsize=16)

plt.subplot(2, 1, 2)
for gamma_2 in gammas_2:
    plt.plot(time, gamma_2)
plt.title('gammas_2', fontsize=16)

Действительно новый $\gamma$ сюжет то, что я ожидал. Ура! Хорошее расследование.

Как рассчитать угол траектории полета γ по вектору состояния?

ооо

пользователь20636

Ответы (2)

Том Спилкер

ооо

Том Спилкер

Том Спилкер

Том Спилкер

ооо

Том Спилкер

ооо

ооо

ооо

Том Спилкер

ооо

ооо

Том Спилкер

Почему мое математическое решение vis-viva оказалось таким близким, несмотря на то, что оно было неправильным? При каких условиях это было бы хорошим приближением?

Сколько delta-v я использовал здесь? Каково «официальное» уравнение для дельта-v от параметрической тяги?

Как работает пример пакета poliastro python «Going to Mars with Python»? Что он делает на самом деле?

Как лучше всего представить Матрицу переходов состояний и как использовать ее для поиска периодических гало-орбит?

Системы координат для векторов состояния

Ограниченная задача трех тел: необходимо создать орбиту вокруг двух массивных тел.

Что именно означает универсальная переменная x и z?

Существуют ли еще точки Лагранжа, если на третье тело от первого оказывается значительное радиационное давление?

Почему угол изгиба гиперболической траектории дает разные результаты?

Библиотека с открытым исходным кодом для расчета изображения из сигналов распределенного набора радиотарелок?