Math model

Движение космического аппарата рассчитывается в трехмерной системе координат. Начало отсчета - центр Земли. На полет влияют силы гравитационного притяжения, аэродинамического сопротивления и тяги. Кроме того, рассчитываются перегрузка и аэродинамический нагрев космического аппарата.

Все силы измеряются в (кг * км)/с^2.

Сила гравитационного притяжения (вектор)
F1 = G * mEarth * mShip * rVector / rScalar^3
G - гравитационная постоянная = 6.67408 * 10^(-20) км^3/(кг * с^2)
mEarth - масса Земли = 5.9742 * 10^24 кг
mShip - масса космического аппарата (вместе с массой топлива) (кг)
rVector - координаты относительно центра Земли (км)
rScalar - расстояние до центра Земли (км)

Плотность P (скаляр)
P = p * M * 10^9 / (R * T)
10^9 нужно для перевода из кг/м^3 в кг/км^3
p - давление (Па) = p0 * e^(-M * g * h * 10^6 / (R * T))
10^6 нужно для перевода g и h из км в метры
p0 - нормальное атмосферное давление на уровне моря = 101325 Па
M - молярная масса сухого воздуха = 0.0289644 кг/моль
g - ускорение свободного падения (км/с^2) = G * mEarth / r^2
h - высота над поверхностью Земли (км)
R - универсальная газовая постоянная = 8.31447 Дж/(моль * К)
T - температура (К)

Сила аэродинамического сопротивления (вектор)
F2 = -P * vScalar * vVector * S^2 / 2
vScalar - скаляр скорости (км/с)
vVector - вектор скорости (км/с)
S - длина ребра космического аппарата (км)

Сила тяги (вектор)
F3 = mLevel * impulse * orientationVector / orientationScalar
mLevel - массовый расход топлива (кг/с)
impulse - удельный импульс двигателя (км/с) orientationVector - вектор ориентации аппарата
orientationScalar - длина вектора ориентации (км)

Угловая скорость (вектор)
I * dw / dt = M
w = previuosW + M * size / I
M - момент силы = d * F (кг * км^2 / c^2)
d - координаты центра вращения (км) - здесь центр координат в одной из вершин куба (3 ребра - 3 оси)
F - равнодействующая сила = F1 + F2 + F3
I - момент инерции (кг * км^2)
size - размер кванта (с)
previousW - предыдущий вектор угловой скорости (км/с)

Аэродинамический нагрев (скаляр)
aH = T + v^2 / 2

Вращение космического аппарата
Из вектора угловой скорости создается кватернион (1), с помощью которого и осуществляется поворот вектора ориентации космического аппарата (2).

(1): (X, Y, Z) - вектор угловой скорости
cosX = cos(X / 2)
cosY = cos(Y / 2)
cosZ = cos(Z / 2)
sinX = sin(X / 2)
sinY = sin(Y / 2)
sinZ = sin(Z / 2)

w = cosZ * cosY * cosX - sinZ * sinY * sinX
x = sinZ * sinY * cosX + cosZ * cosY * sinX
y = cosZ * sinY * cosX - sinZ * cosY * sinX
z = sinZ * cosY * cosX + cosZ * sinY * sinX
q = (w, x, y, z) - кватернион (+ еще делается нормализация)

(2): newOrientation = q * orientation * q'
orientation - предыдущий вектор ориентации
q' - обратный к q

Для расчета траектории движения мы делим время полета на кванты, размер которых задается пользователем, и полагаем, что в течение одного кванта скорость аппарата не изменяется. Скорость на каждом кванте рассчитываем с помощью второго закона Ньютона и формулы a = dv/dt.

vVector = v0Vector + (F1 + F2 + F3) * t/m,

где t - размер временного кванта, v0Vector - скорость аппарата в предыдущий квант времени Разделяя скалярную и векторную части получаем формулу:

vVector = v0Vector * (1 - (P * vScalar * S^2)/2 * t/mShip) + orientationVector * mLevel * impulse * t/mShip - rVector * (G * mEarth)/( rScalar^3) * t