KerboScript в примерах и задачах. Часть 3.5. Предсказание положения на орбите.

В предыдущих частях гайда описаны все основные шаги для того, чтобы отправить аппарат с поверхности Кербина к ближайшему спутнику и обеспечить с ним связь. Однако в каких-то случаях всё может пойти немного не по плану, и орбита после трансфера оказывается не такой, какая надо. В этой части рассмотрим, как можно исправить положение.

Постановка задачи

Аппарат находится на экваториальной эллиптической орбите вокруг Кербина, входящей в сферу действия Муны. Необходимо провести коррекцию орбиты на заданном расстоянии от Кербина таким образом, чтобы в сфере действия Муны перицетр оказался на заданной высоте.

На какой высоте проводить коррекцию?

Обычно чем раньше производится коррекция, тем меньше нужно менять орбиту для достижения нужного результата. С другой сторону, на высокоэллиптической орбите даже малые изменения могут требовать большой ΔV, если их проводить вблизи перицентра. Для примера будем считать, что манёвр надо производить на высоте, равной половине расстояния от Кербина до Муны. Задача расчёта оптимальной высоты проведения коррекции оставляется желающим для самостоятельного решения.

Как производится коррекция?

В общем случае, манёвр представляется в виде суммы трёх ортогональных компонент - касательной (prograde), радиальной (radial out) и нормальной (normal). Касательная и радиальная компоненты лежат в плоскости орбиты, нормальная направлена по нормали к плоскости орбиты (NB: хотя в KSP используется левосторонняя система координат, вектор нормали к орбите там направлен согласно правилу правого винта). Поскольку по условию плоскости орбит аппарата и Муны совпадают, то коррекции в нормальном направлении не требуется. Вектор корректирующего импульса имеет только касательную и радиальную составляющие.

Расчёт момента коррекции

В kOS планирование манёвра осуществляется путём задания момента времени, на который нужно поставить узел. Следовательно, требуется по заданной высоте рассчитать, сколько времени потребуется аппарату, чтобы до неё долететь.
Задача нетривиальная, поскольку движение по кеплеровой орбите неравномерное. Оно подчиняется второму закону Кеплера:

• за равные промежутки времени радиус-вектор, соединяющий планету и аппарат, заметает равные площади
  

К счастью, задача описания такого движения в терминах обыкновенных, а не дифференциальных уравнений была решена самим же Кеплером.
Чтобы понять это решение, нужно разобраться, как определяется положение тела на орбите в пространстве и во времени.
В пространстве положение тела можно определить направлением радиус-вектора, проведённого от тела к планете. Это направление задаётся углом, который называется истинной аномалией и обозначается буквой ν. Стандартным считается отсчёт истинной аномалии от направления на перицентр. Положение на орбите в терминах времени также определяется через время, прошедшее от момента прохождения перицентра. Это время переводится в другой угол, называемый средней аномалией M . По определению, средняя аномалия меняется равномерно по времени, т.е.
M(t) = 2π (t - t_Pe) / T,
где t - текущее время, t_Pe - время последнего (или любого другого) прохождения через перицентр, T - орбитальный период, M - средняя аномалия в радианах.
Таким образом, время, требуемое для движения от точки 1 до точки 2 на орбите, равно
Δt = (M₂ - M₁) T / 2π
Для нахождения связи между истинной и средней аномалиями вводится параметр, называемый эксцентрической аномалией E. Геометрический смысл её проще всего понять по картинке:

Рисунок 1. Истинная аномалия ν, эксцентрическая аномалия E и средняя аномалия M.

С помощью эксцентрической аномалии легко записывается параметрическое уравнение эллипса, если начало координат поместить в его центре:
x = acosE, y = bsinE
Длина радиус-вектора, проведённого из фокуса к телу на орбите (SP на рис. 1):
r = a (1 - e · cosE),
где r - длина радиус-вектора, a и e - большая полуось и эксцентриситет эллипса. Идею вывода этого выражения можно посмотреть в Википедии.
Нехитрыми тригонометрическими преобразованиями показывается связь истинной и эксцентрической аномалий:


Средняя же аномалия вычисляется по уравнению Кеплера:
M = E - e · sinE,
где M - средняя аномалия в радианах, E - эксцентрическая аномалия в радианах, e - эксцентриситет.
При помощи этого уравнения легко находим, сколько времени потребуется, чтобы оказаться на высоте, равной половине радиуса орбиты Муны (или любой другой наперёд заданной высоте). Для получения параметров орбиты используются следующие поля структуры ORBIT :

orbit:semimajoraxis // большая полуось
orbit:eccentricity // эксцентриситет
orbit:trueanomaly // истинная аномалия в градусах
orbit:period // орбитальный период

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

function etaheight {
  parameter h is apoapsis.
  if (h > apoapsis) or (h < periapsis) { return 1/0. } // краш, если высота имеет недопустимые значения
  local rmag is body:radius + h. // расстояние от фокуса
  local ta to orbit:trueanomaly.
  local sma to orbit:semimajoraxis.
  local ecc to orbit:eccentricity.
  local eanow to 2 * arctan( sqrt( (1 - ecc) / (1 + ecc) ) * tan(ta) ). // эксцентрическая аномалия в градусах!
  local manow to eanow - constant:radtodeg * ecc * sin(eanow). // радианы нужно перевести в градусы
  local eaheight to arccos( ( 1 - rmag / sma ) / ecc ). // эксцентрическая аномалия на целевой высоте
  // проверки на случай, если первый проход нужной высоты будет при движении вниз
  if -eaheight > eanow { set eaheight to -eaheight. }
  if eaheight < eanow { set eaheight to 360 - eaheight. }
  local maheight to eaheight - constant:radtodeg * ecc * sin(eaheight).
  return orbit:period * (maheight - manow) / 360.
}

Расчёт корректирующего импульса

Для расчёта приращения скорости, которое приведёт перицентр будущей орбиты к желаемому значению, пойдём по пути наименьшего сопротивления, т.е. градиентного спуска.
Градиентом называют направление, вдоль которого функция многих переменных быстрее всего возрастает. В нашем случае эта функция
f(dV_PRG, dV_RAD) = Pe(dV_PRG, dV_RAD) - Pe₀,
где dV_PRG и dV_RAD - приращения скорости в тангенциальном и радиальном направлениях, Pe₀ - желаемая высота перицентра.
Градиент функции f(x, y) определяется как вектор grad f с компонентами (∂f/∂x, ∂f/∂y).
Алгоритм расчёта манёвра такой:

1. 1) ставим начальный манёвр на заданное время с ΔV = 0
   
2. 2) оцениваем частные производные ∂Pe/∂ΔV_PRG и ∂Pe/∂ΔV_RAD; вычисляем градиент
   
3. 3) скорость изменения функции вдоль градиента равна df/d(grad f) = ((∂f/∂x)² + (∂f/∂y)²)^1/2 = |grad f|
   
4. 4) оцениваем необходимую добавку к дельте как dΔV = (Pe₀ - Pe) / (dPe / d(grad Pe) ), раскладываем её по тангенциальной и радиальной компонентам согласно направлению градиента
   
5. 5) обновляем дельту: ΔV -> ΔV + dΔV; вычисляем перицентр с новой дельтой; если он отличается от нужной величины, переходим к п.2 и вычисляем следующую поправку
   

function midcoursecorr {
  parameter pe0.
  parameter pe0tol is 100. // с какой точностью (в метрах) хотим получить перицентр
  local dvrad to 0.
  local dvprg to 0.
  local corrnode to node(time:seconds + etaheight(Mun:apoapsis/2), dvrad, 0, dvprg).
  add corrnode.
  local newpe to corrnode:orbit:nextpatch:periapsis.
  local pe1 to 0.
  local pe2 to 0.
  local gradr to 0.
  local gradp to 0.
  
  until abs(newpe - pe0) < pe0tol {
  // оценка частной производной по dV_rad
    local dvrad1 to dvrad - 0.001.
    set corrnode:radialout to dvrad1.
    set pe1 to corrnode:orbit:nextpatch:periapsis.
    local dvrad1 to dvrad + 0.001.
    set corrnode:radialout to dvrad1.
    set pe2 to corrnode:orbit:nextpatch:periapsis.
    set gradr to (pe2 - pe1) / 0.002.
    
  // оценка частной производной по dV_prg
    local dvprg1 to dvprg - 0.001.
    set corrnode:prograde to dvprg1.
    set pe1 to corrnode:orbit:nextpatch:periapsis.
    local dvprg1 to dvprg + 0.001.
    set corrnode:prograde to dvprg1.
    set pe2 to corrnode:orbit:nextpatch:periapsis.
    set gradp to (pe2 - pe1) / 0.002.
    local dfdgrad to sqrt(gradr^2 + gradp^2).
    set dvrad to dvrad + gradr/dfdgrad * (pe0 - newpe) / dfdgrad.
    set dvprg to dvprg + gradp/dfdgrad * (pe0 - newpe) / dfdgrad.
    
    set corrnode:prograde to dvprg.
    set corrnode:radialout to dvrad.
    set newpe to corrnode:orbit:nextpatch:periapsis.
  }
}

Бонус: расчёт попадания в Муну без опорной орбиты

Умение вычислять средние аномалии в различных точках орбиты можно использовать для расчёта попадания в Муну (или выхода на сближение с другим кораблём) прямо со старта.
Требуются два запуска - один пристрелочный, второй боевой. В первом запуске аппарат сразу выводится на орбиту с апоцентром выше 12 тыс. км (с типичными ракетами из KSP, скорее всего, этот апоцентр будет достигнут ещё до выхода из атмосферы, т.е. это даже не будет стабильной орбитой). После выхода за границу атмосферы нужно отметить:

• а) сколько времени прошло со старта
  
• б) какой угол дуги аппарат прошёл с момента старта
  
• в) сколько времени осталось лететь до высоты орбиты Муны
  
• г) какой угол осталось пройти до достижения этой высоты
  

Первый пункт решается крайне просто - встроенная переменная missiontime хранит время с момента старта (запуска первой ступени).
Третий и четвёртый пункт решаются описанным выше способом - через аномалии на текущей и на будущей высоте.
Для второго пункта нужно посчитать угол между радиус-вектором аппарата и некоторым фиксированным направлением в момент старта и момент выхода из атмосферы. Такое направление, к счастью, имеется, и хранится в системной переменной solarprimevector.
Таким образом, становится известно, какой угол проходит аппарат с момента старта до высоты Муны и сколько времени это занимает. Откуда уже несложно и высчитать угол упреждения для основного запуска:
Δφ₀ = Δφ_ship - Δφ_Mun = Δφ_ship - 360°×Δt/T_Mun
где Δφ₀ - угол упреждения, Δφ_ship - полный угол, проходимый аппаратом со старта до достижения высоты Муны, Δt - время со старта до достижения этой высоты, T_Mun - орбитальный период Муны.
Реализация этого метода оставляется желающим для самостоятельной работы.

Скрипты с коррекцией орбиты для крафта из прошлой части: