Синтаксически эта программа проще предыдущих: ни векторов, ни циклов, ни своих функций кроме двух. Зато в ней есть три вещи, о которых стоит поговорить отдельно.
if (R < 1e-4f) R = 1e-4f;
Запись 1e-4 означает $1 \cdot 10^{-4}$, то есть 0.0001. Буква e читается как «умножить на десять в степени».
| Запись | Значение |
|---|---|
1e6 |
1 000 000 |
2.5e-3 |
0.0025 |
6.67e-11f |
гравитационная постоянная, float |
Научная запись делает порядок величины явным и облегчает сравнение малых и больших чисел. Ошибка в показателе степени всё ещё возможна, поэтому числовые константы следует проверять.
Суффикс f на конце по-прежнему нужен – без него число будет double.
В программе нет типа «вектор на плоскости», хотя задача про векторы целиком. Вместо этого каждая векторная величина разложена на пару переменных:
float x, y; // радиус-вектор
float v_x, v_y; // скорость
float a_x, a_y; // ускорение
Это самый простой способ, и для двух измерений он вполне рабочий. Плата за простоту – каждую операцию приходится писать дважды, по разу на компоненту, и легко ошибиться, скопировав строку и забыв заменить x на y. Такая опечатка – одна из самых частых в численных расчётах, и ищется она тяжело: программа работает, картинка выглядит правдоподобно, но физика неверна.
Как делают, когда компонент становится много: заводят структуру
struct Vec2 { float x, y; };
и учат её складываться и умножаться на число (это называется перегрузкой операторов). Тогда шаг метода пишется одной строкой на вектор. Для одной планеты это лишняя работа, а вот в задаче о молекулярной динамике, где частиц сотни, такая структура уже есть.
Про структуры – в статье «Свои структуры: struct».
Формула $a = \alpha / R^2$ при $R = 0$ содержит деление на ноль. В реализациях с арифметикой IEEE 754 результатом обычно становится inf или nan; полагаться на это как на обработку ошибки нельзя.
Дальше inf умножается на координату, получается inf или nan, и через несколько кадров все переменные состояния становятся нечисловыми. На экране это выглядит как исчезнувшая планета: рисовать nan нельзя, и график просто перестаёт что-либо показывать.
В программе стоит простейшая защита – ограничитель снизу на R. Приём называется клампингом (от clamp – зажать) и встречается постоянно:
if (R < R_MIN) R = R_MIN; // вручную
R = std::max(R, R_MIN); // то же самое, короче
R = std::clamp(R, R_MIN, R_MAX); // с обеих сторон, C++17
Первые две строки ограничивают R только снизу. std::clamp ограничивает значение и снизу, и сверху; если верхнее ограничение не требуется, эквивалентом первой строки является std::max. Последние две записи требуют #include <algorithm>.
Клампинг – приём языка, а не решение задачи: почему он здесь только затычка и чем центр поля лечат по-настоящему, разобрано в статье о численном методе.