Спутник на дисплее

Квант 1989, №10 с. 60-64 Тег 5ZN9

Спутник на дисплее

Н. Н. НОСКОВ, С. К. СТОЛБОУШКИН

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

Р. Фейнман

Все мы когда-то строили модели кораблей или макеты домиков и городов. А не хотелось ли вам сделать действующую модель Вселенной или, хотя бы, Солнечной системы? Такую модель построить вполне возможно, и очень удобным инструментом для ее создания является компьютер.

Мы начнем с простого — с моделью движение спутника Земли. Сначала займемся теоретической постановкой задачи, т. е. установим физические законы, определяющие движение спутника, а затем составим алгоритм, по которым можно писать программу для ЭВМ.

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

При моделировании неизбежны упрощения. Поэтому для начала мы предположим, что в космосе находится только два тела: Земля и спутник. Чтобы построить модель их взаимодействия, нужно знать физические законы, определяющие их движение. Эти законы два: второй закон Ньютона и закон всемирного тяготения, открытий им же. Кроме того, нужно задать исходные данные. С этого и начнем.

На спутник действует единственная сила — сила прятяжения Земли. Ее можно вычислить по формуле закона всемирного тяготения

$F = G \frac{M_1 m_c}{R},$

где $G$ — gravitationная постоянная, она равна $6,67 \cdot 10^{-11} \text{ m}^3 / (\text{кг} \cdot \text{с}^2)$ ;

$M_3$ — масса Земли, равная $5,98 \times 10^{24} \text{ кг}$ ;

$m_c$ — масса спутника.

По второму закону Ньютонаускорение $\vec{a}$ спутника равно: $\vec{a} = \vec{F} / m_c$ . Далее, зная координаты $\vec{S} \vec{\theta} = (x\theta; y\theta)$ и скорость $\vec{v} \vec{\theta}$ спутника в некоторый момент времени и применяемый кinemатрику

$\vec{v} = v\theta + \vec{a} \Delta t, \quad \vec{S} = S\theta + \vec{v} \Delta t,$

(1)

казаловы бы, можно определить скорость и положение спутника через время $\Delta t$ . Но это не так! Сила, действующая на спутник, да и скорость $\vec{v}$ за время $\Delta t$ изменится. Что же делать?

Выход есть. Предположим, что промежуток времени $\Delta t$ , через которого нужно вычислить положение и скорость спутника, мал. Тогда за это время расстояние от спутника до Земли изменится незначительно и изменением силы $\vec{F}$ и скорости $\vec{v}$ за этот промежуток времени можно будет пребежать. Нужно только разумно выбрать величины $\Delta t$ . Результатом вычисления будет новые значения координат и скорости спутника. Для них можно повторить эти же вычисления при том же $\Delta t$ и получить очередные величины координат и скорости. Последовательно применя такую схему вычислений, мы получим координаты спутника в последовательные моменты времени $T1, T2, T3, ...$ , а соединив эти точки отрезками прямых, мы получим ломанную линию, которая приближен-

но представляет траекторию движения спутника (рис. 1). Такие ломаны называются ломанными Эйлера, а сам метод вычислений также носит имя Леонарда Эйлера — великого российского ученого XVIII века.

Метод Эйлера применяется не только для нахождения траекторий небесных тел. Он используется для решения широкого класса уравнений, называемых дифференциальными уравнениями. Наше исходное уравнение (второй закон Ньютона) — тоже дифференциальное. Причем метод Эйлера применяется для решения уравнений, которые точными методами не решаются.

Теперь, зная метод Эйлера, можно составить подробную схему вычисления траектории спутника на околоземной орбите.

Поскольку орбита спутника расположена в одной плоскости с центром Земли, установим систему координат так, чтобы плоскость орбиты совпада с плоскостью $xy$ . В таком случае координата $z$ всегда равна нулю и ее можно нерасматривать. Начало системы координат помещим в центр Земли (см. рис. 2). Пусть расстояние от центра Земли до спутника равно $R$ , а его скорость относительно Земли — $\vec{v}$ . Для того, чтобы однозначно определить положение спутника, одной величины $R$ недостаточно. Нужно либо к $R$ добавить угол $\alpha$ , под которым виден спутник из начала координат, либо вместе пары $R$ , $\alpha$ задать пару координат $x, y$ .

Если известно расстояние $R$ и угол $\alpha$ , то координаты $x, y$ можно определить по формулам $x=R \cos \alpha, y=R \sin \alpha$ .

Скорость спутника и сила пряжения Земли — величины векторы, разложим их на компоненты по осям $x, y$ . Из подобия соответствующих треугольников на рисунке 2 видно, что

$\frac{F_x}{F} = -\frac{x}{R}, \quad \frac{F_y}{F} = -\frac{y}{R},$

отсюда

$F_x = -|F|x/R, \quad F_y = -|F|y/R.$

Если на спутник действует только пряжение Земли, то второй закон Ньютона можно записать в виде

$m_c a_x = -G \frac{M_3 \cdot m_c}{R^2} \frac{x}{R}, \quad m_c a_y = -G \frac{M_3 \cdot m_c}{R^2} \frac{y}{R}.$

Сократив на $m_c$ и подставив значение ускорения в формулы kinematики (1), получим формулы для вычисления новых значений координат и компонент скорости через время $\Delta t$ :

$x^\mathrm{H} = x + v_x \Delta t, \quad y^\mathrm{H} = y + v_y \Delta t, \quad v_x^\mathrm{H} = v_x - G \cdot M_3 x \Delta t/R^3, \quad v_y^\mathrm{H} = v_y - G \cdot M_3 y \Delta t/R^3.$

Формулы (2) представляют собой готовую схему алгоритма для ЭВМ. Реализация этого алгоритма мы и заявимся. Но сначала выясним, как задаать шаг $\Delta t$ . Существуют специальные методы расчета шага с тем, чтобы через дительное время погрешность вычисления координат не превышала наперед заданного значения. Мы поступим проще — подберем значение шага опитным путем. Для этого нужно выбрать произвольное $\Delta t$ и найти решение для этого значения. Затем уменьшить шаг в два раза и снова найти траекторию. Если решение для нового $\Delta t$ не отличимо от предыдущего, то будем считать, что выбор шага удовлетворит.

Теперь можно приступать к реализации модели на компьютере.

Составление программы

У нас есть алгоритм вычисления траектории спутника, т. е. набор формул (2), позволяющих нахождение и скорости спутника в любой момент времени. Но если применять

эти формулы в СИ, то все величины будут иметь больше значения, соответствующие космическим масштабам. Поэтому выгодновести свою систему единиц, в котором произведение $GM_3$ равно единице. Это получится, если принять за единицу длины радиус Земли, а за единицу скорости — первый космический скорость $v_{1k} = \sqrt{GM_3/R_3}$ . В этом случае время будет изменять в единицах $R_3/v_{1k}$ . Соотношение между единицами новой системы измерений и СИ для Земли получится следующие: единица длины — $6,38 \cdot 10^6$ м, скорости — $7,8 \cdot 10^3$ м/c, времени — $0,82 \cdot 10^3$ с, а формулы (2) принимают вид
$x^H = x + v_x \Delta t,$
$y^H = y + v_y \Delta t,$
$v_x^H = v_x - x \Delta t/R^3,$
$v_y^H = v_y - y \Delta t/R^3.$
Формулы стали более пригодными для вычисления и, главное, универсальными. Они «работают» не только для Земли, но и для любого космического тела, естественно, в системе единиц, в которой единица длины — радиус этого тела, а единица скорости — первый космическая скорость для этого тела.

Посленее, что нам осталось сделать перед написанием программы, — это уяснить, что такое графический экран и как на нем расположить нашу модель. Дисплей в графическом режиме высокого разрешения представляет собой прямоугольное поле точек размером $m \times n$ . Ось $x$ на нем всегда расположена слева направо, а ось $y$ может расположаться как сверху вниз, так и сизу вверх. Давайте считать, что ось $y$ идет сизу вверх и что $m=n=256$ . Чтобы вывести на экран точку, нужно задать ее координаты и цвет. Чтобы нарисовать линю, следует указать координаты начальной и конечной точки. Для изображения окружности следует указать координаты центра и радиус. Операторы для этих действия будем обозначать так:
точка (X, Y), цвет;
линия (X1, Y1) — (X2, Y2), цвет;
круг (X, Y), R, цвет.

Однако в качестве экранных координат нельзя просто взять физические координаты. Как вы помните, за единицу мы взяли радиус Земли, а в центр Земли помещили начало координат. Тепер нам надо помещить центр Земли в центр экран (точку с координатами ( $X_\pi$ , $Y_\pi$ ), например (128, 128)) и все уведомить в $k$ раз. Получаются следующие формулы: если физические координаты точки (спутника) $x, y$ , то координаты на дисплее

$XD = X_\pi + k*x, \quad YD = Y_\pi + k*y.$

Величину $k$ мы определим позже.

Тепер можно переодить непосредственно к составлению программы. Сначала построим общую схему, а затем будем ее детализировать. В группе приближении программы для движения спутника Земли можно записать в следующем виде:

начало
ввод начальных данных;
подготовка графического дисплея;
нарисовать Землю и спутник;
повторять
вычислить новые значения координат;
нарисовать отрезок траектории и спутник в новом положении;

конец повторяка конец.

Первые две инструкции не вызывают проблем. Пусть спутник в начальной момент времени ТФ расположен на оси $x$ , тогда начальное координату $y$ надо положить равной нулю и оставаться расстояние $x$ в единицах физической системе координат, т. е. в единицах земного радиуса.

Начальную скорость тоже несложно в виде пары чисел VX, VY. Однако оказывается, что удобнее ввести значение модели скорости и ее направление (угол с осью $x$ ), а затем преобразовать в компоненты прямой системы координат по формулам

$VX = V*\cos (y\theta*\pi/180),$
$VY = V*\sin (y\theta*\pi/180).$

Выражения в скобках нужны для перевода угла из градусов в радианы,

потому что во всех языках трионометрические функции вычисляются только в радианах.

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

Тепер изобразим Землю и спутник. Для того чтобы рисовать, требуется определить масштабный коэффициент $k$ . Мы предлагаем задать его по формуле $k=100/x$ . Тогда изображение Земли и спутника (с экранами координатами XD, YD) реализуется операторами

круг (X $_\pi$ , Y $_\pi$ ), $k$ , голубой;
XD←X $_\pi$ + $k^*x$ ; YD←Y $_\pi$ ;

точка (XD, YD), красный.

Масштабирование в таком случае получается автоматически. Если вы решили запустить спутник с поверхности Земли, т. е. $x=1$ , то изображение Земли займет почти весь экран, его радиус будет сто экранных единиц. Если же спутник помещить на расстоянии $x=10$ , то Земля займет кругок радиусом 10 точек.

Осталось определить величины шага по времени T $\theta=\Delta t$ . Экспериментально определено, что значение в интервале 0.1—0.2 выполняет удовлетворительно.

Все подготовительные действия законены, приступим к моделированию полета спутника на основании формул (3).

Повторять
XS←XD; YS←YD;
X←X+VXT $\theta$ ;
Y←Y+VYT $\theta$ ;
XD←X $_\pi$ + $k^*X$ ;
YD←Y $_\pi$ + $k^*Y$ ;

Линия (XS, YS) — (XD, YD), синий;
Точка (XD, YD), красный;

$R\leftarrow\sqrt{X^2+Y^2}$ ;
VX←VX-XT $\theta/R^3$ ;
VY←VY-YT $\theta/R^3$ ;

Конец повторя.

Как вы, наверное, договорались, XS и YS — старые экраны координаты спутника. Причем мы не просто проводим линию или ставим точку, а делаем и то, и другое. В результате по экрану будет дивигаться красная точка, оставляя за собой синий след. При вывод линии можно указать только координату конечной точки, если начальную точку компьютер помнит сам.

Программа практического конечна, осталось решить, как должно повторять цикл. Можно задать число повторений заранее, но не исключено, что работа с моделью надоест нам раныше, чем кончится цикл или, наоборот, цикл закончится в самый интерсный момент. Поэтому удобнее было бы перевать цикл по нашему желанию. Например, нажатим какой-нибудь клавиши. И еще одно условие окончания цикла надо добавить. Дело в том, что если R<1, ничего не помещает программе продолжать работу. А это означает, что спутник временя в Землю и продолжает движение в ее недрах. Поэтому ведем так называемую бульевскую переменную УПАЛ, которая принимает только два значения: истина и ложь. И будем вычислять ее так: УПАЛ←(R<1). Условие повторя в этом случае:

ПОВТОРЯТЬ ПОКА (НЕ УПАЛ) И (КЛАВИША ≠“S”)

Цикл будет продолжаться до тех пор, пока не нарушится хотя одно из условий.

Мы составили программу на гипотетическом языке. Ее легко перевести на язык любого компьютера.

Тепер поговорим о том, что вы можете самостоятельно изменить и дополнить в программе. Первая часть изменений касается естественно оформления. Если спутник упал, то можно в точке падения изобразить взрыв. Да и сам спутник можно рисовать не точкой, а в виде фигурки (владекошему искусством спрайтов это

не составит труда). Можно через 10 или 20 циклов выводить на экран текущее значения скорости и расстояния до Земли. Это первое, что приходит в голову. Надемся, вы придумеете что-нибудь еще.

Помимо художственных дополнительный, можно изменить и саму модель. В этом случае простora для фантазии еще большее. Если в формулах для VX и VY заменить минусы на плюсы, то закон всемирного тяготения превращится в закон всемирного отталкивания. Тем самым модель спутника превращится в модель взаимодействия одноменного заряженных частиц. А можно содействовать несуществующую Вселенную. Сила в законе тяготения обратно пропорциональна квадрату расстояния. А что получится, если она будет обратно пропорциональна кубу или первой степени расстояния? Нам несложно содействовать и эту ситуацию — достаточно в формулах заменить $R^3$ на $R^4$ или на $R^2$ .

А вот следующее изменение потребуют бóльших трудов. Мы создали модель спутника Земли, но не можем им управлять. Кроме того, в космосе существуют и другие тела, на которые можно летать. Интересно было бы на «Союз» с «Фобоса» высадить «Попрыгунчик» на Фобос. Этот мы займемся в одном из следующих номеров журнала.

Список читателей:
приславших правильные решения
(Начало см. на c. 38)

48, 55, 58, 59; П. Молодов (Ломоносов) 55; С. Мурин (Врест) 59; В. Мытько (Ленинград) 48, 52, 55, 59, 60; А. Насибов (Шекинский р-н АэССР) 52; А. Насиров (Шекинский р-н АэССР) 52; Н. Немировская (Киев) 48, 52, 58; А. Николаев (Киев) 48, 49; И. Новиков (Калининград) 58; Р. Новоселов (Москва) 55, 58, 59, 61; Д. Омецинский (Киев) 48, 49, 51, 58—61; В. Онищук (Минск) 58, 59; Г. Оренер (Киев) 58; А. Орловский (Киев) 48, 49, 51, 58—61; Д. Островский (Днепропетровск) 60; А. Павлоцук (Киев) 49, 52, 58; Ю. Пайко (Врест) 58; П. Пензин (Куйбышев) 58—60; М. Перепечко (Москва) 55; С. Петухов (Горький) 52; И. Пилосин (Алма-Ата) 52, 58; Н. Плахтий (Киев) 58; Т. Подаляко (Ленинград) 58, 61; С. Подрятов (Киев) 49, 52; Э. Помарский (Ташкент) 49, 52; В. Полушук (Канев) 51, 61; С. Полышин (Харьков) 52, 58—61; С. Пономаренко (Новосибирск) 49; С. Погорочин (Ижевск) 52, 55, 58, 60; Е. Призант (Одеcca) 51, 55, 58, 59, 61; Е. Проскурина (Киев) 58; В. Пузанов (Донецк) 51; И. Птайкин (Москва) 55; А. Ревский (Москва) 61; А. Распереза (Врест) 48, 49, 52, 55, 58, 61; И. Рассадин (Минск) 52; Р. Ролинский (Киев) 52; А. Руевский (Ивано-Франковск) 52; Р. Рузев (Ханкинский р-н Хорезмской обл.) 52; И. Рябова (Харьков) 52, 59; А. Савин (Горловка) 58, 59, 61; А. Савченко (с. Мошны

Черкасской обл.) 58, 59, 61; Н. Самойлова (Киев) 58; А. Самуенко (Семипалатинск) 48, 49, 55; С. Севченко (Минск) 52; А. Скабелин (Баранович) 48, 51, 55, 58—60; Т. Скрипник (Киев) 48; В. Солдатов (п. Краснореченский Пrimорского кр.) 55, 59; Ю. Спектор (Киев) 48, 51, 52, 58, 59, 61; Ю. Старостин (Якутск) 55; М. Субботин (Стайый Оскол) 49; А. Сувханов (Чарджоу) 59; Н. Сулейманов (Шекинский р-н АэССР) 52; И. Сысовев (Семипалатинск) 55, 59, 61; С. Табальцев (Москва) 58; В. Тамошина (Вильнюс) 48, 49, 51, 52, 59; А. Тенекеджи (Воронеж) 58, 59; С. Тимашов (Алма-Ата) 48, 49, 55; М. Титов (Киев) 48, 55, 58—60; С. Тихонов (Воронеж) 48, 51; Р. Тодоров (Варна, НРВ) 52; А. Тутов (Рига) 51; Д. Украинский (Москва) 58; А. Усинский (с. Птичья Ровенской обл.) 48, 49, 52, 60, 61; Ю. Фенок (Киев) 58; Л. Френклах (Воронеж) 58, 60; А. Фридлянд (Саратов) 48, 49, 55; Г. Хачатран (Арташатский р-н АрмССР) 48; Р. Христок (Киев) 48, 51; Е. Чашечкина (п. Черноголовка Москвской обл.) 48, 49, 52, 55, 58—60; И. Чемолосов (Москва) 60, 61; Л. Чернышеев (Москва) 49, 59, 61; Д. Чокин (Алма-Ата) 48, 49, 52, 55, 58, 61; В. Чуев (Стайый Оскол) 52; Е. Шааров (Грозный) 48, 59; А. Шазинян (Ростов) 59; Т. Шарафутдино (Усть-Каменогорск) 48, 49; Е. Швец (Черновцы) 49, 52, 60; С. Шинкевич (Верезники) 48, 49, 51, 52, 58—60; С. Шубский (Рест) 48; И. Шуляк (Киев) 52, 55, 58—60; К. Шурунов (Куйбышев) 58, 59; Н. Шербина (Артемовск) 52; М. Энтит (Тула) 48, 49, 51, 52, 55, 58—61; В. Юрасов (Волгодонск) 55; И. Ягольцер (Черновцы) 52, 55, 58, 59.