После того, как мы собрали небольшое приложение для решения достаточно абстрактной задачи, перейдем к такой теме, как основы небесной механики. С ней связаны задачи исследования движения спутников, выполнение орбитальных маневров, оценка области видимости ИСЗ с поверхности Земли, исследование движения спутников в точках Лагранжа. А заодно посмотрим на реализацию новых вычислительных алгоритмов и способов визуализации.
Начнем с математической модели.
В нашей задаче основными являются гравитационные силы, вычисляемые по закону Всемирного тяготения. Принцип прост - все тела притягивают друг друга с силой, прямо пропорциональной их массам и обратно пропорциональной квадрату расстояния между ними. Соответственно, каждый объект мы опишем через его массу, скорость и координату.
Данными нашей модели будет массив из объектов, а методы модели будут определять взаимное расположение этих объектов, находить действующие на них силы, вычислять ускорения и на основе кинематических параметров строить орбиты.
Потом мы дополним нашу модель различными дополнительными факторами - атмосферой, тягой, давлением излучения Солнца, рассмотрим движение в поле тяготения сразу нескольких тел, но пока - заложим базис.
Ключевой метод объекта - расчет сил, которые действуют на него под влиянием другого объекта, и которые он сам на него оказывает. Метод accelerations принимает списки кинематических параметров двух объектов и массу второго объекта, с которым происходит взаимодействие.
Сначала мы определяем разность координат, затем квадрат расстояния, затем - ускорение, действующее на основное тело ( мы сразу же делим силу притяжения на массу основного тела) и готовимся к тому, чтобы определить проекции силы тяготения на оси нашей системы координат.
Стоит обратить внимание на параметр lightweight. Мы специально ввели его, чтобы уменьшить объем вычислений модели для расчета орбит малых объектов, притяжением которых можно пренебречь. Например, так можно исследовать движение спутника в "троянской точке" системы "Земля - Луна". Если этот параметр является истинным, то мы не вычисляем силу, которую наш объект оказывает на другие тела.
Мы можем вычислять ускорения и знаем скорости каждого объекта. Пора приступить к оценке ускорений и скоростей всех объектов в модели. Для этого будем последовательно проходиться по всем объектам и вычислять силы, которые он оказывает (и которые оказывают на него) остальные объекты. При таком подходе мы заметим несколько нюансов, благодаря которым можно сократить в 2 раза объем вычислений:
- Вычислив силы для первого тела, мы будем знать, с какой силой он будет действовать на 2, 3 и так далее объекты (и они на него). Это значит, что для 2 объекта (и всех после него) мы уже не должны считать взаимодействие с первым объектом. Распространим это правило и получим, что если у нас есть всего M объектов, то при расчете сил, действующих на N объект мы можем пропустить первые N-1 объектов.
- Нам не нужно проводить вычисления для последнего объекта - все силы, действующие на него и оказываемые им, вычисляются при прогонке предыдущих объектов.
Таким образом, расчет состоит из основного прогона от первого до предпоследнего тела (i = 0...M-1), внутри которого для i-го мы совершаем прогон от i+1 до M тела. Мы немного уменьшили размерность задачи с O(2) до 0(2)/2. По сравнению с нашим "наивным" подходом есть и более эффективные методы, основанные на группировке нескольких тел (метод Барнеса - Хата), но мы пока сосредоточимся на физическом, а не алгоритмическом аспекте задачи.
Теперь у нас есть производные кинематических параметров, и можно приступать к численному интегрированию уравнений движения. Воспользуемся проверенным методом Рунге-Кутта IV порядка точности с фиксированным шагом. Схема метода выглядит следующим образом:
K0 = df(x, t)/dt
K1 = df(x + K0 * dT/2, t)/dt
K2 = df(x + K1 * dT/2, t)/dt
K3 = df(x + K2 * dT, t)/dt
f(x, t + dT) = f(x) + (dT/6) * (K0 + 2K1 + 2K2 + K3)
Численное интегрирование происходит в методе integrate. Получив на вход шаг интегрирования и количество шагов, метод запускает схему Рунге-Кутта и на каждом шаге сохраняет результаты расчетов в массиве result.
Расчетная схема выглядит достаточно громоздко, но это тоже продиктовано требованиями быстродействия. Попробуем заменить прямое перечисление элементом массива лаконичным map, в начале вычислений запустим счетчик производительности методом performance.now(). Этот метод дает нам текущую временную метку с точностью до микросекунд.
Сохраним время запуска метода в переменной timeStamp1, а перед выдачей результата - найдем разность во времени (в миллисекундах) между запуском и окончанием интегрирования сначала для for-цикла, а затем для map.
В прямом цикле на вычисления (всего 5200 шагов) уходит ~ 80 миллисекунд, в случае map ~ 115 миллисекунд.
Разумеется, в ответственных за интерфейс участках программы нет нужды заниматься подобной дичью и "развертывать" массив, а циклы for повсюду, где это возможно, заменяются на map, reduce, forEach или filter.
Физический движок развернут, и можно загрузить в него первые данные. Начнем с тестовой задачи на вращение искусственного спутника Земли вокруг нашего глобуса. Ее легко проверить аналитическим расчетом
На загрузку документа создается модельная система, запитанная данными Земли и тестового ИСЗ (высота орбиты - 100 км, масса - 100 кг, скорость - 7848 м/с). Для проверки по формуле V1 = Sqrt(G0 * m / R) определим для высоты в 100 км первую космическую скорость - 7848.26 м/с, период обращения - 5180 секунд.
Вывод результатов модельного эксперимента в консоли браузера показывает, что спустя 5180 секунд наш спутник будет на расстоянии всего 3 км от точки, в которой он начал движения.
Если бы мы уменьшили шаг интегрирования или выводили результаты расчета с меньшей скважностью, то добились бы точного совпадения, но видно, что разность координат между начальной и конечной спустя один виток - меньше расстояния, которое спутник пролетает за один шаг интегрирования.
Мы смогли рассчитать самый основной случай - движение спутника в центральном поле тяготения!
