Ещё один путь к Земле

Представляем вашему вниманию оригинальную реализацию программ серии "Путь к земле" в виде веб-страниц со скриптами на сайте Станислава Володарского.

"На этом сайте есть всё необходимое, чтобы повторить перелёт «Кон-Тики» с Луны на Землю и решить многие другие задачи космонавигации."

Кроме программ-симуляторов, лунолетчикам может потребоваться книга "Задачи для юного космонавта".

Undefined

Комментарии

Похоже, что сетевой Лунолет не очень совместим с самым популярным кораблем современности - iPad-ом :)
У меня не получилось вводить команды

---------------------------
Истина где-то рядом
aldebaran.ru/author/samurov_vitaliyi/kniga_dozvonitsya_do_devyi/

Работает прекрасно на iPad'е. Просто управление не столь очевидное, советую инструкцию прочитать.

Если кратко, то вводишь исходный комплект данных в левом столбике, потом два дефиса, потом команды.

Да, работает. Ну что ж, RTFM, как говорится

---------------------------
Истина где-то рядом
aldebaran.ru/author/samurov_vitaliyi/kniga_dozvonitsya_do_devyi/

Тексты выложены в хорошем текстовом виде — с индексами и греческими буквами, в Юникоде. Где требуется, подобран кегль, курсив, жирный шрифт. На вскидку веб-страничку сложно отличить от ксерокса, т.к. вёрстка идеальная, с картинками из ТМ.

Переработаны и выложены все программы цикла, хотя и без исходного кода. Программы серии «ОС» объединены в одну, механизм стыковки как-то доработан. Вместе с «Лунолётом-3» выложена программа «Лунолёт-3.1» с лучшей точностью вычислений. Опять же, без исходного текста сразу не сказать, какие внесены изменения. Даны ссылки на ключевые места в Сети, включая мои переводы эпопеи на QBasic.

В-общем, работа Станислава Володарского — профессиональная. Дизайн минималистичный, но всё необходимое есть. Хотя сайт может быть улучшен, он сейчас — лучшее воплощение «Пути к Земле» в Интернете.

Поля веб-страничек я сделал бы побольше. Несколько неудобней то, что при работе программы не видно ни текста, ни инструкции. Впрочем, оконная система большинства операционных систем позволяет получить на экране и то, и другое. Также, играя за компьютером, их можно вывести на планшетку, и наоборот.

Для развития художественной части можно посоветовать озвучку — первую попытку в XXI веке восстановить старую магнитофонную запись Александра Перепёлкина. В принципе, OS X Lion самостоятельно зачитывает русский текст по Option-ESC, а бесплатный CoolReader 2 позволяет под Windows записывать озвучку в файл. Профессиональные аудиокниги также сейчас делаются на потоке. В идеале художественная часть может стать туториалом в игре, проигрываясь в виде видео или комикса, позволяя в нужный момент пользователю вводить рекомендуемые команды.

Неплохо, если программная часть сможет отображать траекторию полёта в графическом виде, координатные сетки Oxt или Oxy. Это можно сделать как в отдельной клетке таблицы, так и на фоне журнала. Когда я совершал свои первые полёты, миллимитровка существенно помогала осознавать происходящее. Кстати, она где-то была упомянута в КЭИ, хотя сходу точного места не нашёл.

Дальнейшее развитие программ может отражать следующие заключительные слова из «Техники-молодежи» 1986 №5, которые мы легко можем превратить в пророческие: «Персональный компьютер позволяет также вести игру в реальном масштабе времени. Это, конечно, приближает ее к действительности. Пока администрация КЭИ такими программами не располагает и надеется здесь на помощь читателей.»

Тексты:
Тексты сканировал (ещё до меня) Сергей Павлович Хлынин. На его сайте http://epizodsspace.no-ip.org/ большое количество материалов по космонавтике. Когда мне понадобились рассказы из ТМ, я списался с Сергеем Павловичем и, с его разрешения, переверстал их. Они по-прежнему размещены на его сайте.

Исходный код:
Он есть, хотя и спрятан. Странички с программами написаны на JavaScript и скачиваются на клиентскую машину, где и работают. Вот прямые ссылки на код:
http://slavav.ru/way_to_earth/programs/Lunolet3.js,
http://slavav.ru/way_to_earth/programs/Lunolet3.1.js. Несколько замечаний:
1. Интерпретатора кода БЗ-34 у меня нет. Код для калькулятора переписывался в виде программы на JavaScript, которая имитировала его исполнение в браузере (отличия - другая точность и все операции со сверхчислами реализовывались "вручную" по описанию из журналов). Затем он проверялся (сверялись результаты на калькуляторе МК-52 и компьютере). Затем код рефакторился в "нормальный" код на JavaScript. Затем я заменял посадочный/стыковочный блок (если программа зациклится на калькуляторе, то пользователь вмешается; если на компьютере, то данные будут потеряны). Затем еще одна сверка и всё. Варианты решений из журнала можно прогнать на компьютере с теми же результатами (это была первая цель проекта).
2. Математическая модель третьего лунолёта точная, но надо выбирать маленький шаг по времени для получения высокой точности. Говоря технически - используется метод Эйлера. Он имеет первый порядок точности. Выбор шага интегрирования (времени манёвра) ложится на пользователя. Чем больше шаг - тем больше ошибка. Самое неприятное в этом для меня то, что нельзя один манёвр разбить на две "половинки" с тем же реактивным ускорением. Точность возрастёт и лунолёт окажется в другом месте. При посадке "на пределе" это неприятно.
Лунолёт-3.1 использует метод интегрирования пятого порядка с автоматическим контролем точности. То есть, точность решения гарантируется. Но модель та же: оригинальный третий лунолёт с уменьшением времени манёвра будет "сходится" к решению Лунолёта-3.1. Вторая цель проекта - получение физически достоверного решения навигационных задач.
Ещё в Лунолёте-3.1 задаётся не расход топлива, а реактивное ускорение. Надоело в каждом манёвре подбирать расход под ускорение. Всё-таки в навигации ускорения и скорости важнее расхода топлива.

Художественная часть:
Я не художник, и пока не планировал доработок по этой части (разве что сделать зеленое на чёрном и добавить звездно-лунные фоны :).
Графика для траекторий есть в планах (кстати, посоветуйте технологию - Flash, HTML5, генерация картинок на серверной стороне, Java-апплет?. Хочется, что бы траектории были доступны на максимальном количестве платформ).
С тьюториалом действительно надо что-то сделать. Возможно, для Лунолёта-1 будет всплывающее окошко, которое "проведёт" пользователя по рассказу "Истинная правда".

Реальный масштаб времени:
Есть идеи, как сделать управление в браузере мышью в динамике. Это подойдёт для первого и второго лунолётов и атмосферы-2. Вообще взлёты, посадки и стыковки делать в динамике прикольно и занятно. В длительных манёврах нужно дополнительное навигационной оборудование - траектометры всякие, которые помогут пилоту быстро выбирать углы, ускорения, скорости. Тут для меня многое не ясно.

Рад увидеть ваш ответ, Станислав. Сайт "Эпизоды космонавтики" знаю давно, достойный. Спасибо за исходный код. Впервые увидел хорошую математику, реализованную на JavaScript. Впрочем, адреса этих .js-файлов можно легко достать из Page Source на HTML. Меня просто смутило, что html'ки страниц отдаются php-скриптом.

Пятый порядок, да ещё явное численное интегрирование это сильно. Улучшения третьего "Лунолёта", опубликованные здесь уважаемым basvic, использовали третий порядок:
http://pmk.arbinada.com/node/214#comment-877
http://pmk.arbinada.com/node/320

Автоматический контроль точности при интегрировании был у проф. Дьяконова для советских ПМК и развит мной для ЭКВМ здесь:
http://pmk.arbinada.com/node/369
http://pmk.arbinada.com/node/408

Использование реактивного ускорения при задании манёвра -- вещь интересная. Неясно только, могут ли так двигатели работать. Это ведь серьёзное изменение в конструкции впрыскивателя топлива, который был рассчитан на постоянный расход. Любопытно, как реальные манёвренные двигатели работают.

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

JavaScript сейчас идеальная технология и ваше решение работает везде -- даже на MacOS X и iOS. Дальше идёт HTML5, который поддерживается Firefox и iOS. Adobe в её битве со стандаром постепенно идёт на уступки со своим собственническим и тяжеловесным флэшем. Сложно придумать аргументы, зачем кому-то в России надо вставать на её сторону:
http://blog.store-apple.ru/flash-ios/

В реальном времени, думаю, можно уже сейчас провести эксперимент (судя по некоторым опубликованным фразам, такова и была задумка Михаила Пухова):

  • встроить тикающий на экране таймер, отсчитывающий время манёвра
  • после окончания счёта таймера игровое время останавливается, ожидание следующей команды
  • возможность быстро пропускать время, по нажатию соответствующей кнопки
  • возможность подготовить данные для нового манёвра, во время предыдущего -- ввести в экранную форму
  • нажать-утопить кнопку "следующий манёвр" -- который будет введён сразу после окончания действия таймера текущего манёвра

Интересно будет впервые пройти "Истинную правду" в реальном времени...

Пятый порядок:
Я использовал метод Dormand-Prince. Он рассчитывает два решения в конце шага и по разнице выдаёт рекомендацию по уменьшению/увеличению шага интегрирования. Таких методов много, отличаются они только набором коэффициентов.

К Дьяконову отношусь с уважением, но с тех пор техника интегрирования ушла вперёд. Точно не знаю, но думаю что Дорманд-Принс на МК-152 даст лучшие результаты.

Постоянное реактивное ускорение:
Это было самое сомнительное изменение в Лунолёте 3.1. Мотивация такая: темп расхода топлива устанавливается пилотом и поддерживается насосами. Несложно в эту цепочку встроить управление, которое будет снижать расход в процессе манёвра. В конце концов, пилоту не запрещается разбивать манёвры на секундные интервалы и так моделировать переменную производительность насоса для получения максимального ускорения на всём протяжении манёвра. Далее к управлению приделываем акселерометр - и готова обратная связь для поддержания постоянного ускорения. Мне пришло в голову, что именно так мог быть устроен лунолёт Михаила Пухова - отсюда и постоянство ускорения при манёвре.

Или немного по другому: Можно считать, что "журнальная" математика неточна, так как меняется ускорение при постоянном расходе. А можно считать, что ускорение постоянно, но надо точнее считать расход (самое обидное, что последнее изменение можно было вставить во все калькуляторные лунолёты - формула Циолковского очень просто считается). Оба изменения имеют право на существование.

Я понял, что ваш лунолёт летает на постоянном расходе. Мне нужен ваш прогон. Я изменю математику в 3.1 (это не сложно, так как никаких аналитических формул (кроме формулы Циолковского) там нет и переменное ускорение ляжет в ту же схему без затруднений) и проверю совпадение результатов.

Интересно, как пилоты лунных модулей при посадке меняли тягу? Плавно или ступенчато?

Спасибо за Дорманд-Принс, было интересно ознакомиться. До этого я использовал таблицы коэффициентов лишь в методе Гаусса-Кронрода. Но здесь, насколько я понял, разбиение отрезка интегрирования идёт на равные части. Ещё и адаптивный метод. И товарищ под рукой, который это уже реализовал. :-) Пока не уверен, что этим в ближайшее время займусь, но ведь ЭКВМ не только у меня есть.

Лунолёт и его математику развивал другой участник форума -- basvic. Надеюсь, он здесь объявится и даст необходимый прогон. Сам я прогнать не могу, т.к. для ЭКВМ этой программы нет. Только для разных Касио.

Если тяга постоянна, уравнения решаются проще. Но ваш подход с численным интегрированием позволяет сделать "Лунолёт" с постоянным расходом без аналитического решения. Действительно, интересно сравнить результаты.

Кстати, вы учитываете изменение центробежного и кориолисова ускорений с высотой, во время манёвра -- или ограничились переменной массой? Нам именно это не дало разработать точную аналитическую модель:
http://pmk.arbinada.com/node/214#comment-883

> ускорение постоянно, но надо точнее считать расход (самое обидное, что последнее
> изменение можно было вставить во все калькуляторные лунолёты - формула
> Циолковского очень просто считается)

Ого! Оказывается, в XXI веке мы можем совершенствовать ПМКшные "Лунолёты"! Постоянную тягу, кстати, на КА действительно делают. Возможно, именно из-за удобства расчёта их движения. В динамическом "Лунолёте" можно сделать переключатель -- постоянная тяга или расход.

> Интересно, как пилоты лунных модулей при посадке меняли тягу? Плавно или ступенчато?

Ага! И есть много людей, которым не лениво задавать НАСА всё более точные вопросы... :-)))

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

Ускорения:
... пересчитываются на каждом шаге интегрирования. Благодаря этому и переменному шагу я могу летать целыми витками. Или имитировать ионный двигатель (малая тяга, большое время). Аналитическое решение, к сожалению, невозможно, если у вас есть реактивное ускорение (даже если он постоянное).

Я-то надеялся, что решение можно свести к интегралу. Который, если не получается взять аналитически, брать численными методами.

Получается же, что надо просто подбирать как можно более мелкий шаг и менять на нём массу, ускорения и т.д. И точного решения на 12-14 знаков -- пусть не аналитического, а численного нет? Или же ваш Дорманд-Принс учитывает все изменяющиеся факторы, включая разные гравитационные и инерционные ускорения при отдалении от планеты?

Летать витками интересно. При больших вычислительных ошибках не факт, что возникнет такая вещь, как стабильная орбита.

Дорманд-Принс родственник Рунге-Кутта и работает также. Раличаюся только таблицы коэффициентов.
В каждой точке расчёта (в пространстве и времени) все ускорения вычисляются в зависимости от текущих координат и скоростей (это верно и внутри шага экстраполяции).
В отличие от Рунге-Кутта конечная точка вычисляется по двум линейным комбинациям пятого и шестого порядков. Разница в решениях даёт ошибку. Если ошибка превосходит ожидаемую, то решение выбрасывается и шаг уменьшается (формула для пересчёта шага учитывает пятый порядок ошибки). Если нет, то вычисляется ожидаемый шаг для следующей итерации - благодаря этому шаг может увеличивается в процессе интегрирования.
В процессе манёвра учтываются изменения всех характеристик кроме расхода топлива, который непосредственно на динамику не влияет, только через реактивное ускорение, которое просто вычислять как функцию времени или считать постоянным.
Работа с ошибками позволяет дать гарантии на точность. При ошибке 10^-9 в секунду (м/с и м) через 11 суток полёта решение будет отличаться от истинного не более чем на один миллиметр. И это гарантируется, лишь бы формулы для ускорений были верны (ещё надо бы прибавить ошибки округления которые накапливаются, но можно показать, что они ещё меньше).
Формулы я проверял путём сличения решений в полярной и прямоугольной сетке координат (у меня есть Лунолёт-3.2 решающий ту же задачу в прямоугольных координатах). На прямоугольной решётке формулы гораздо проще, сделать там ошибку малореально. Результаты совпали с точностью до погрешности.
Стабильная орбита возникает. Вот вариант: выход на орбиту; три витка половинками, чтобы показать периселений и апоселений; Затем 5 витков; затем 50 витков.
Единственное чего я ещё не проверял это сохранение энергии и момента количества движения. Проверка за мной.

Про расход здесь написал (и стёр) глупость. Надо будет посидеть подольше над формулами. Значит, изменение массы корабля, особенно сложное при постоянном реактивном ускорении, не влияет ни на гравитационное, ни на инерциальные ускорения. Следовательно, не оказывает влияние ни на скорость, ни на координаты корабля. Если же расход постоянен, то реактивное ускорение можно решить аналитически и независимо от основной системы. Я правильно вас понял?

А как в Дорманд-Принсе вычисляется ошибка если переменных несколько? Насколько я понял, ошибка это число, а не вектор.

Сила тяготения пропорциональна массе корабля. Следовательно его гравитационное ускорение не зависит от его массы. Инерциальные ускорения (центробежное и кориолисово) - это фикция для компенсации непрямоугольности системы координат. Их нет, например, если вы считаете все в прямоугольной решётке (это называется инерциальная система). Так что вы поняли правильно.
Ошибка берется как максимум ошибки по всем координатам и скоростям. Дорманд-Принс возвращает мне вектор, я из него делаю скаляр (длину). Её и контролирую, и, как следствие, контролирую ошибки по координатам и скоростям.

То есть чисто математическое решение. Складываем метры и метры в секунду. :( Относительная погрешность здесь вряд ли поможет, т.к. точность определения координат должна зависеть от высоты лишь при посадке и взлёте.

Кстати, каков ваш подход к урезанию количества десятичных знаков координат при межпланетных перелётах? Если он имеется, конечно. Даже 384.400.000 метров от Луны жрут целых девять знаков. В результате точность измерения координат значительно хуже, чем скорости и ускорения. Про расстояние до Солнца вообще молчу...

Математическое решение:
Можно оценивать ошибки покомпонентно. Если хотя бы одна координата или скорость дают большую ошибку, то пересчитываем шаг. Такой контроль даёт те же результаты, что и вычисление длины вектора ошибки (в действительности проверка длины вектора строже контролирует ошибки).

Точность:
Размеры солнечной системы около 40 ае, 1ае = 150 * 10^6 км. Размер солнечной системы в миллиметрах будет 40 * 150 * 10^6 * 10^3 * 10^3 = 6000 * 10^12 = 6 * 10^15. То есть мне нужно 16 знаков мантиссы, чтобы садится на все планеты с точностью до одного миллиметра. Арифметика компьютера работает с точностью между 16 и 17 знаками. На первый взгляд достаточно.
На второй взгляд - нет. Во-первых полезут ошибки округления (надо добавить, скажем, три знака), во-вторых интегрирую я с контролем абсолютной ошибки 10^-9 метра. Чтобы это работало на краю солнечной системы, нужно ещё шесть знаков. Итого 25 знаков. Такой разрядности под рукой нет.
Решений два: первое - двойная точность (32 знака), написанная вручную. Это медленно и не нужно. Второе - перемещать начало системы координат вместе с КА (только не непрерывно, а шагами определённого размера). Это то что надо - расстояния до всех близких предметов считаются максимально точно. Расстояния до далёких предметов точно можно не считать - ошибки проявятся только в гравитационных ускорениях, но они будут очень малы, фактически незаметны.
Я таких штук ещё не пробовал, но думаю, что получится.
Здесь важна ещё такая штука: в навигации точные координаты и скорости КА не известны. Их определяют с большим трудом с помощью наблюдений с поверхности небесных тел. Последнее время GPS почти решил задачу (шаттлы определяются по GPS, которая для этого отчасти и создавалсь). Да только при первом подлёте к Плутону никакого GPS там не будет.
Знание точных координат и скоростей КА - не очень реалистичный элемент "лунолётов".

> Если хотя бы одна координата или скорость дают большую ошибку, то пересчитываем шаг.

Считаю, что достаточно оценивать не отдельные компоненты векторов, а модуль ошибки координаты и модуль ошибки скорости. В этом случае физический смысл ясен.

> Решений два: первое - двойная точность (32 знака), написанная вручную. Это медленно и не нужно.
> Второе - перемещать начало системы координат вместе с КА (только не непрерывно, а шагами определённого размера).

По большому счёту, это одно и тоже решение. Хранить каждую координату корабля в двух регистрах и учитывать связанную с этим математику. Я так на МК-152 раскидывал 14-значное число из стека по двум 12-разрядным регистрам. Если же не полностью полагаться на встроенную арифметику, а расписать самому операции, можно точность до 20-24 разрядов увеличить.

Тут ещё вопрос, как рассчитывать движение планет и спутников. Если это 32 знакоместа прямоугольной сетки, им просто можно разрешить вращаться по Ньютону и это будет учитывать взаимное влияние и друг на друга, и на планетолёт. Пару лет назад я видел старшеклассника, реализовавшего этот проект на Паскале. :-) В противном случае нам придётся разбивать полёт на множество мелких этапов, как в ТМ, и использовать "простые формулы", которые так любил Лунный Коршун. Кто и как в этом случае будет учитывать взаимное расположение небесных тел, непонятно.

Есть, конечно, вариант сделать планеты и спутники отдельными объектами, расписав для них формулы движения, принятые у астрономов. При этом планетолёт будет единственным телом, летающим среди них по Ньютону.

> Знание точных координат и скоростей КА - не очень реалистичный элемент "лунолётов".

Блок, который отвечает за физику, это одно. Блок ввода-вывода, который отвечает за реалистичность кабины самого корабля, это другое.

> Считаю, что достаточно оценивать не отдельные компоненты векторов, а модуль ошибки координаты и модуль ошибки скорости. В этом случае физический смысл ясен.
Согласен. Но не всё диктуется физическим смыслом. Мне надо выдать три знака после запятой через 10 суток полёта. Это чистая математика.

> По большому счёту, это одно и тоже решение.
Вы правы. Трудности могут возникнуть возникнут при симуляции больших промежутков времени (большой объём вычислений).

> Тут ещё вопрос, как рассчитывать движение планет и спутников.
Нет, только не по Ньютону! 10 тел - это 45 пар взаимодействий (спутники и Луну я опустил). А ведь есть готовые модели, которые всё равно будут нужны, чтобы стартовать симуляцию.

> Есть, конечно, вариант сделать планеты и спутники отдельными объектами, расписав для них формулы движения, принятые у астрономов. При этом планетолёт будет единственным телом, летающим среди них по Ньютону.
Я бы так и делал.

> Мне надо выдать три знака после запятой через 10 суток полёта. Это чистая математика.

Не совсем чистая. Если мы задаём координаты в километрах, три знака после запятой это одно. В метрах, другое. В миллиметрах, третье. Точно также задаётся и необходимая точность для скорости, по-разному для км/ч и м/с. При задании точности важно учитывать размерность.

> Заманчиво реализовать его (или вообще табличные методы с контролируемой точностью) для ЭКВМ.

Можно вас попросить написать эту систему дифуров, решение которой нужно для «Лунолётов» — так, чтобы мне было проще отвлечься от физики и погрузиться в стоящие на полочках талмуды численных методов?

Пока я не вижу принципиальных сложностей, мешающих начать решать эту систему на МК-161. Опыт успешного решения табличным методом с контролируемой точностью есть — хотя, как вы правильно заметили, более лёгкой задачи. Да и сдавал я все эти Рунге-Кутты в МГУ, кажись на Паскале. От вашей поддержки, правда, тоже не откажусь. Давно это было, в 90-е. :-)

Есть две системы координат, в которых можно решать задачу: полярные и прямоугольные. "Калькуляторные" лунолёты летали в полярных координатах: точность у них выше (точность расчётов можно обменять на скорость расчётов) но и формулы сложнее. С какой начать - выбирайте сами.

Прямоугольные координаты (не вращаются, ноль в центре масс планеты):

gc - гравитационная постоянная небесного тела (квадрат радиуса умножить на ускорение свободного падения на поверхности)

dt - шаг по времени

x, y - координаты корабля в осях
vx, vy - проекции скорости корабля на оси
rax, ray - проекции реактивного ускорения на оси

r2 = x^2 + y^2
r = sqrt(r2) // расстояние от корабля до центра планеты

cs = x / r // косинус направления планета-корабль в осях
sn = y / r // синус направления планета-корабль

gax = -cs * (gc / r2);
gay = -sn * (gc / r2);

ax = ray + gax;
ay = ray + gay;

vx_new = vx + ax * dt
vy_new = vy + ay * dt

x_new = x + vx * dt
y_new = y + vy * dt

При программировании важно помнить что rax и ray - функции времени (при постоянном расходе модуль ускорения растёт) и координат корабля x, y (корабль сохраняет свой угол тяги в полярных координатах, а в прямоугольной системе этот угол меняется во время манёвра).

Чтобы учесть меняющийся угол, дополним нашу модель:

a - модуль реактивного ускорения (функция времени)
da - направление реактивного ускорения в полярных координатах (константа, 0 - от центра планеты)

at = a * sin(da) // проекция ускорения на местную горизонталь (трансверсальное ускорение)
ar = a * cos(da) // проекция ускорения на местную вертикаль (радиальное ускорение)

cs = x / r // косинус направления планета-корабль в осях
sn = y / r // синус направления планета-корабль

rax = ar * cs - at * sn // перевод ускорений в прямоугольные координаты
ray = ar * sn + at * cs

Есть две системы координат, в которых можно решать задачу: полярные и прямоугольные. "Калькуляторные" лунолёты летали в полярных координатах: точность у них выше (точность расчётов можно обменять на скорость расчётов) но и формулы сложнее. С какой начать - выбирайте сами.

Полярные координаты (не вращаются, ноль в центре масс планеты):

gc - гравитационная постоянная небесного тела (квадрат радиуса умножить на ускорение свободного падения на поверхности)

dt - шаг по времени

x - горизонтальная координата (угол места в радианах)
y - радиус (расстояние до центра планеты)
vx - проекция скорости на местную горизонталь (трансверсальная скорость)
vy - проекция скорости корабля на местную вертикаль (радиальная скорость)
rax - трансверсальное реактивное ускорение
ray - радиальное реактивное ускорение

трансверсальное ускорение складывается из реактивного и кориолисова
радиальное ускорение складывается из реактивного, центробежного и гравитационного

ax = rax - vy * vx / y
ay = ray + vx * vx / y - gc * y^2

vx_new = vx + ax * dt
vy_new = vy + ay * dt

x_new = x + vx / y * dt // преобразование в радианы!
y_new = y + vy * dt

И опять rax и ray - функции времени. Но не функции координат:

a - модуль реактивного ускорения (функция времени)
da - направление реактивного ускорения в полярных координатах (константа, 0 - от центра планеты)

rax = a * sin(da) // проекция ускорения на местную горизонталь (трансверсальное ускорение)
ray = a * cos(da) // проекция ускорения на местную вертикаль (радиальное ускорение)

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

Формулы понятны и близки. Но в таком примере это обычное линейное приближение. Пока не очень вижу крупную картину, саму интригу. :-) Можно составить, сто раз уже делал, такой вычислительный блок для одного шага по dt:

1. вычисляем реактивное ускорение rax, ray
2. вычисляем ускорение ax, ay
3. вычисляем новую скорость vx, vy
4. вычисляем новые полярные координаты x, y

Можно даже этот блок в цикл взять, с постоянным dt или изменяя его по какому-либо закону (здесь этот Дорманд-Принс зарылся?).

Стоит этот блок кодировать или решение системы требует чего-то более хитрого? Также, требуется ли сохранять старые vx, vy, x, y — или это не обязательно и можно всё записывать в старые регистры? По математике получается странно, т.к. vx ≠ dx/dt, ведь x это радианы, а vx это м/с.

И вспомнил ещё один свой вечный вопрос. :-) Это da относительно чего — на звезду, планета (система координат) вращается? Если da задаётся относительно вертикали или векторов скорости/ускорения, проекции реактивного ускорения очень даже будут зависеть от координат.

Ваше представление верно. Если вы реализуете цикл, то получете интегрирование системы диффуров прямым методом Эйлера. При уменьшении шага он будет сходится к точному решению. Все численные методы решения используют эту идею: из точки пространства отступаем вдоль производной, чтобы найти новую точку пространства. Чем изощрённее метод, тем больше раз он выполняетс этот элементарный шаг.

Берём Рунге-Кутта 4 порядка .
Первый прогноз соответствует методу Эйлера.
Второй прогноз считается для средней точки на отрезке интегрирования (пользуемся первым прогнозом).
Третий прогноз также для средней точки, но учитывает результаты второго прогноза.
Четвёртый прогноз считается для конца шага интегрирования.

Каждый прогноз вычисляет смещение по координата и скоростям.

Затем считаем взвешенную сумму и "передвигаем" корабль.
Шаг сделан. Сложная схема гарантирует, что если решением задачи является многочлен 4-ой степени, то задача будет будет решена точно.
Кроме того, если шаг уменьшить в два раза, и проделать эти операции два раза, то точность вырастет примерно в 16 раз (2^4).

Дорманд-Принс считает решение сразу по двум схемам. Разница между двумя решениями помогает оценить "расстояние" до точного решения. В принципе похожую схему можно сделать из Рунге-Кутты: считаем с шагом один раз, считаем с двумя половинными шагами. Оцениваем разницу в решениях. Корректируем шаг в зависимости от разницы. Недостаток тут будет только один - две трети расчётов выбрасываются на ветер (рассчёт с половинными шагами используется только для контроля ошибки).
Дорманд-Принс для оценки ошибки делает только одну лишнюю экстраполяцию. Сравним эффективность:

Рунге-Кутта с половинками: 12 прогнозов из них 4 идут в результат.
Дорманд-Принс: 6 прогнозов из них 5 идут в результат.

Кроме того у ДП порядок выше - он будет давать меньшую ошибку на больших интервалах.

Блок вам кодировать стоит - это элеменарный кирпич любой схемы (кроме того правильность удобно контролировать рассчётами по Эйлеру).
Старые координаты и скорости надо сохранять. Сами методы требуют 4-6 копий. И если вы не удовлетворены точностью, то придётся выбросить шаг и начать новый (меньший).

vx = y * dx/dt - с размерностью всё в порядке. Но надо понимать, что наша система криволинейная. Поэтому такие странности.

da - угол откладываемый от местной вертикали (управление взято из старого "лунолёта"). Мы его считаем (в реальности для этого нужен микропроцессор с обратной связью и системой ориентации).

Сам блок написал, дело нехитрое:
http://the-hacker.ru/2011/lun-103.mkl
http://the-hacker.ru/2011/lun-103.mkp
http://the-hacker.ru/2011/lun-103.htm

Файлы в UTF-8, в исходнике перевод строки 0x0a. Пока не оптимизировал. Файлы должны читаться MK.EXE, ВК-6 и другими программами, а также загружаться в железный МК-161. С углом da разберёмся позже. Пусть пока будет, как в ТМ.

Взять этот блок в цикл несложно, но вопрос с шагом dt (постоянство или закон изменения которого меня по-прежнему волнуют) и количеством повторений. Также нужна формула a(t) и какие-нибудь начальные данные, для тестов.

Насколько я понял, под «прогнозом» вы имеете в виду не однократный подсчёт всех этих величин моим блоком, как функции dt, а уже конечный результат выполнения написанного блока в цикле, причём N = T/dt раз. Предлагаю пока сосредоточиться на делании одного такого прогноза (чтобы код мог посчитать любой прогноз из 6 требуемых), а потом уже написать подбор шага на основании 6 прогнозов. Промежуточные значения (x, y, vx, vy) в каждом из циклов нам ведь не нужно запоминать, только конечные — когда прогноз уже посчитан?

Чем отличается подсчёт разных прогнозов в Дорманд-Принсе? Только значениями N и dt? Или придётся делать два варианта — для прогноза по Эйлеру и Рунге-Кутту?

Для простоты можно начать с Эйлера. Рецепт выбора шага по времени такой:

e - ошибка в секунду
error - оцениваемая ошибка в процессе расчётов
DT - длительность манёвра
dt - длительность шага внутри манёвра

Мы хотим чтобы error полёт.

Ускорение пока будем считать постоянным.

Манёвр DT как правило будет разбит на много маленьких шагов dt. "Прогноз" - это шаг на dt или его часть. Метод Эйлера делает единственный прогноз на шаг.

Дорманд-Принс делает шесть прогнозов на dt со сложными коэффициентами. Давайте пока отложим сложные методы и отработаем контроль ошибок на Эйлере. Затем я бы попробовал Рунге-Кутта (потребуется ещё уточнить формулы для расчета error и new_dt. А там и до Дорманд-Принса рукой подать будет.

Сразу адаптивный метод предлагаете. :-) С другого конца хотел взяться, ну да ладно. Можно и по-вашему сделать, за 2-3 этапа. Не трудно. Пока задача мало отличается от интегрирования функции. Для надёжности перед кодированием поищу в книжечках метод Эйлера для системы дифуров, а потом перейдём к более быстрым и точным.

Но я бы взялся менять не dt, а m = DT/dt — тогда решение всегда будет заканчиваться точно на границе DT. Насколько я понимаю, мы последовательно вычисляем F(t), F(t+dt), F(t+2dt)… F(t+m*dt)=F(t+DT). При этом каждый последующий набор переменных F(t+k*dt) рекурсивно опирается на предыдущий результат F(t+(k-1)*dt), то есть погрешности накапливаются и повторно использовать значения F(t+k*dt) для чётных значений k, уже вычисленные при предыдущем значении m, мы к сожалению не сможем. За ошибку вы предлагаете считать удвоенный модуль разницы между F(t+DT), посчитанных для m и 2*m. Я правильно изложил?

new_dt вычислять заранее, для следующего манёвра, считаю бесполезным. Т.к. манёвры разные и независимые. Можно каждый раз брать m=2 и удваивать, пока не достигнем нужной точности. Почему вы предлагаете взять более грубое решение, непонятно. Когда я считал интегралы, применял формулы линейной экстраполяции на основе двух наиболее точных значений: I ≈ (a*Im+b*I2m)/(a+b), что практически даром даёт почти такую же точность, как расчёт следующего решения I4m. Можно даже сохранять прошлые значения и брать квадратную/кубическую экстраполяции, но само по себе I2m должно быть неплохим приближением.

Большей точности, если это не связано с дополнительными вычислениями, я не боюсь — да и ерунда это по сравнению с заниженным значением new_dt, которое может возникнуть в вашем варианте.

Да, сразу адаптивный. И вот почему: Верхний уровень алгоритма писать сложнее. И его будет удобнее отладить на простом Эйлере. Сам метод интегрирования на одном шаге можно воспринимать как чёрный ящик: на вход вы ему подаёте точку и dt. Он возвращает новую точку и оценку ошибки. По ошибке вы принимаете решение и вычисляете новый шаг.

про DT/m вы изложили правильно. Возможно построить адаптивный метод таким способом. Например удваивать m пока нас не удовлетворит накопленная ошибка в конце DT. Собственно это всё что нужно для адаптивного метода (так решал задачу интегрирования Коши). Ну это вы и сами понимаете.

Но времени вы потратите много: постоянный шаг внутри манёвра плохо соответствует меняющимся условиям: манёвром может быть оборот по вытянутой эллиптической орбите. У апогея dt увеличивается, у перигея уменьшается. Тоже верно под тягой: маленькие скорости - большой шаг, скорость растёт - шаг уменьшается.

new_dt вычисляется для следующего шага внутри того же манёвра. Новый манёвр делает первую попытку при dt = DT. По ней, скорее всего неудачной, вычисляется new_dt подходящий под условия.

Про грубое решение: я беру самое грубое решение, которое укладывается в заданную точность. Большей точностью пользователь воспользоваться не сможет (пока мы ему её не покажем, а мы её не показываем). Но на большую точность тратится время пользователя - ценный ресурс, а котором программисты часто забывают :).

Заниженное значение new_dt продержится только один шаг. В конце шага метод оценит ошибку и увеличит шаг не дожидаясь окончания манёвра.

Простите, так за время одного манёвра DT шаг dt всё-таки меняется?? Мы вычисляем интегралы с dt и dt/2 не последовательно (для всего манёвра DT), а одновременно — корректируя dt на ходу. При этом, в рамках одного манёвра длительностью DT, заданного человеком, могут быть использованы совершенно разные значения dt — и полученные с помощью них данные войдут в конечный результат, если все локальные ошибки в пределах нормы.

Я правильно вас понял?

Когда писал похожую программу, шаг выбирал в зависимости от ускорения, действующего в данный момент на тело. Чем оно больше, тем шаг меньше. Использовал метод Рунге-Кутты 4 порядка, он очень простой в реализации: http://spiff.rit.edu/richmond/nbody/OrbitRungeKutta4.pdf

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

Файл оказался на удивление полезным, хотя и пестрит отпечатками. Действительно, такой подход точнее ПМКшного. К сожалению, меня уже развратили адаптивными методами, позволяющими оценивать точность приближения — а не считать ошибки «случайными возмущениями со стороны неучтённых масс». :-)

Адаптивность не зависит от метода экстраполяции. Тот метод, что вы напишете, нормально отработает и с Эйлером и с Рунге-Куттой и с Дормандом-Принсом. Всякий раз будет получено точное решение, только времени на это уйдёт разное количество.

Вы же сами писали, что при адаптивном Рунге-Кутте приходится делать лишние вычисления.

Я ожидаю, что в Дорманд-Принсе при вычислениях p1, pp1 и pp2 можно будет повторно использовать часть значений, уже вычисленных моим блоком — а не гонять его всё время заново. Там синусы вычисляются, между прочим.

ДП считает p1 и pp2 одновременно. Но на принципиальную сторону метода это не влияет. Да, РК будет считать медленнее ДП, но зато много быстрее Э. И РК проще реализовать чем ДП - там таблица меньше.

Я веду вас путём программиста - отлаживаем внешний блок на медленном простом методе (Э), пробуем более продвинутый метод (РК), и наконец штурмуем вершину - табличные методы (ДП). Заодно мы убедимся в работоспособности отдельных блоков, будет легче искать ошибки, и сравним производительность всех методов.

Синусы можно вынести из циклов расчётов - они все не меняются внутри манёвра. Но это мы забегаем вперёд. Сперва работающий алгоритм, затем алгоритмическая оптимизация (выбор метода). Затем низкоуровневая (вынос синусов из циклов и т.п.).

Уже загрузил, попробую собрать.
На вашей программе я обязательно буду летать.

Если делать JavaScript страничку, то можно обойтись без куки. Даже интерактивность сохранится.

Отзыв будет. Надо получить на моём и вашем "лунолётах" одинаковые результаты. Этого, думаю, мы сумеем добиться.

Планируете ли вы сделать арифметику и функции эмулятора побитово совпадающими с "железным" калькулятором?

Эмулятор собрал. Завтра попробую запускать программу.

Закодировал Эйлера — и даже вшил в программу исходные данные, соответствующие вашему примеру отсюда:
http://pmk.arbinada.com/node/804#comment-4102

Если у вас запускается MKL2MKP, можете эти исходные данные менять по своему усмотрению и компилировать программу lun-103.mkl заново. Брать на прежнем месте:
http://the-hacker.ru/2011/lun-103.mkl
http://the-hacker.ru/2011/lun-103.mkp
http://the-hacker.ru/2011/lun-103.htm

Правда, я пока сам этот монстр не запускал и не проверял — боюсь, да и спать уже надо. Так что ошибок, наверное, куча. Но уже есть, где их исправлять. :-))

Код совершенно не оптимальный и перестраховочный, т.к. пока не знаю, что придётся менять и переставлять при переходе на более тяжёлые алгоритмы. Зато хорошо комментирован. Пришлось использовать почти 50 регистров и пять страниц, хотя этот код лишь один манёвр способен рассчитать! То ли ещё будет… :-)

Приведу пример: DT = 10c

p0 - начальная точка
вычисляем позицию p1 для начальной точки p0, t = 0 и dt = DT = 10:
p1 = f(p0, 0, 10)

Вычисляем туже позицию в два шага:
pp1 = f(p0, 0, 5)
pp2 = f(pp1, 5, 5)

по pp2 - p1 оцениваем ошибку. Предположим, что она оказалась в пять раз больше требуемой. Тогда вычисляем новый шаг:

dt = dt / 5 (dt = 2)

и снова

p1 = f(p0, 0, 2)
pp1 = f(p0, 0, 1)
pp2 = f(pp1, 1, 1)

снова проверка ошибки. Предположим, что она подходит, и даже в два раза меньше ожидаемой. тогда делаем следующий шаг:

dt = dt * 2 (dt = 4)

p2 = f(p1, 2, 4)
pp1 = f(p1, 2, 2)
pp2 = f(pp1, 4, 2)

и так далее

Да, вы поняли меня правильно:
:loop
p1 = f(p0, t, dt)
pp1 = f(p0, t, dt/2)
pp2 = f(pp1, t + dt/2, dt/2)
er = 2 * |pp2 - p1|
if (er

Большое спасибо. Да, такое закодировать просто.

1. Какое ограничение снизу для dt предлагаете? Сверху это DT-t, насколько я понимаю.

2. Также вот этот фрагмент ведёт себя как-то странно (а с размерностями там вообще ужас):
> er = 2 * |pp2 - p1|
> if (er dt = 0.9 * er / e

Предположим, мы получили очень большую ошибку er. По логике, надо проделать тот же шаг с меньшим dt.

Правильно, что условие не выполняется и время t не увеличится. Ведь er > e*dt. Следовательно, dt er/e. Но присваивание запишет в dt величину явно меньшую. Возможно ли, что правильно написать dt = 0.9 * e / er ?

3. Также повторю ещё один свой вопрос — какова аналитическая запись a(t) ? Без кодирования этой формулы мой вычислительный блок бесполезен. Простите, что сосредотачиваюсь на математике и программировании МК-161, а не физической модели.

под DT я понимал длительность манёвра.

1. если t0 - время начала манёвра, t - текущее время, то dt

Модуль ускорения для модели постоянного расхода (постоянной тяги):

// время с начала манёвра
dt = t - t0

// полная масса корабля = сухая + топливо - то что израсходовано
fm = M + m - (dt / DT) * mass

// мнгновенное реактивное ускорение
ra = C * m / (DT * fm);

Спасибо. Количество переменных растёт, как на дрожжах. :-) Может ли быть такое, что последняя формула на самом деле такова (насколько я понял, C это скорость истечения продуктов сгорания):

ra = C * mass / (DT * fm)

Хотелось бы услышать подтверждение. Не дожидаясь его, я применил это изменение и выложил расчётный блок по старым адресам:
http://the-hacker.ru/2011/lun-103.mkl
http://the-hacker.ru/2011/lun-103.mkp
http://the-hacker.ru/2011/lun-103.htm

Вы правы

Извините за назойливость, но поскольку уже есть точные формулы для реактивного ускорения, скорости и координаты, почему бы их не использовать http://pmk.arbinada.com/node/214#comment-883

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Наверное потому, что я ещё в 2008 году не до конца понял, как совмещать точное аналитическое решение реактивного движения с нерешённой аналитически задачей гравитационных и инерциальных сил — а также почему такой аналитически-численный гибрид будет работать и в чём его преимущество. Если теорию Дорманд-Принса я ещё немного понимаю и доверяю ей, с комбинированными методами ранее не сталкивался вообще.

Но я помню о вашем решении формулы Циолковского. Когда закончим данный проект с точным численным решением (о котором мечталось с детства), с удовольствием пройду у вас курс обучения на тему комбинированных аналитически-численных методов или изучу указанную вами литературу на эту тему.

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

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

В третьем "лунолёте" криволинейная система координат. Ваше решение хорошо подойдёт для прямоугольной системы. Но, оказывается, что криволинейная система численно лучше (даёт в два раза меньшие ошибки, проверено на двух "точных" "лунолётах"). Математики недаром её ввели в те времена когда стоимость расчётов была очень велика (расчёты могли проводится только в головах очень умных людей).

Так что неясно где истина. Окончательное (и убедительное для всех участников процесса) решение можно получить только сделав тесты. Выиграет метод, который:

1. При уменьшении некоторого параметра сходится к точному решению (параметром может быть шаг интегрирования или оценка ошибки сверху).
2. Обгоняет другие методы - при одинаковых отклонениях от истинного решения считает быстрее.

Первый пункт даёт корректность. Второй - эффективность. Если хотите, можно согласовать условия более формально и вступить в борьбу в поисках истины.

В настоящий момент у меня не так много свободного времени, чтобы вступать в соревнования. Но я подумаю. Может вечерок-другой на это потрачу.

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Вот эту линейность системы я и ставил под сомнение, т.к. потенциал гравитационного поля явно нелинейный. Про силы инерции же вообще не говорю. Уважая ваше решение, я просто не знаю, как его в данной задаче применять.

Другими словами, когда призрачный безгравитационный, но точный реактивный «Лунолёт» идёт с тягой под углом 45°, перемещаясь на 10, 20 и 30 метров, гравитационное поле изменяется нелинейно. Благодаря этому сила тяжести, действующая на реальный «Лунолёт» будет на каждом из этих этапов в 10, 20 и 30 метров отличаться от силы тяжести, действующей на ваш точный реактивный. Вычислить по вашей формуле, например, реактивную составляющую сдвига ∆r несложно. Но добавлять её под квадрат радиуса при учёте гравитации? Без веских аргументов на это у меня руки не поднимаются. Сдвиг на ∆r неправильных (без учёта других сил, но с учётом тяги) метров по радиальной оси приведёт к неправильному и нелинейному вычислению силы тяжести — причём явно отличному от правильного классического сдвига на dr. У меня нет причин верить, что отличающийся результат будет более точным!

Для меня вообще не очевидно, почему наложение на гравитационный потенциал точного безгравитационного реактивного решения, при котором Лунолёт прилетит в совершенно другую точку с отличной гравитацией, будет правильным. Также я не знаю, какой числовой фокус может сделать такое наложение правильным. Отсутствие теории численных вычислений, в которых некоторые слагаемые в системе дифуров решены аналитически, делает меня неспособным к кодированию этой задачи, кроме как сомнительным методом повезло / не угадал.

Это не обязательно означает, что подобную теорию создать нельзя или что я не приветствую работу в этом направлении. Возможно, вы нашли хорошую тему для чьей-либо научной работы, даже докторской. До комментария Станислава я даже тешил себя мыслью, что эта теория уже кем-то разработана и её можно взять/скачать и изучить. Просто на данном этапе можно гарантированно написать программу, использующую уже проверенные методы решения дифуров — и одно это существенно улучшит ситуацию. Потом, при необходимости, будем её улучшать дальше. Возможно, опираясь на уже сделанный шаг, пробовать что-то, что другие ещё не пробовали.

Точное решение не учитывает меняющихся кориолисова, центробежного и гравитационного ускорений. Совместить точное решение с приближёнными методами будет сложно. Для приближённых методов доказаны теоремы об их точности, которые дают уверенность, что при правильно выбранном шаге ошибка будет меньше некоего порога. Как ни странно, хотя гибридный метод кажется лучше, для него такие теоремы ещё не доказаны. Для этого есть две причины: 1. гибридный метод может оказаться хуже, 2. просто не умеем доказывать/не брались за это дело.
Положение таково, что интуитивной убеждённости в том что один метод лучше другого мало.

>Точное решение не учитывает меняющихся кориолисова, центробежного и гравитационного ускорений.

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

Переменные:
exy — желаемая ошибка координат на секунду полёта
ev — желаемая ошибка скорости на секунду полёта
dt_min — минимальный шаг по времени
R — радиус планеты
DT — время манёвра
∆t — шаг по времени
t — время от начала манёвра
p0 — вектор координаты и вектор скорости корабля
p1 — вектор координаты xy1 и вектор скорости v1 корабля на следующем шаге
pp1 — вектор координаты и вектор скорости корабля на половинном шаге
pp2 — вектор координаты xy2 и вектор скорости v2 корабля через два половинных шага
erxy — оценка ошибки координат корабля
erv — оценка ошибки скорости корабля

Алгоритм:
1. Начальная инициализация p0, DT, exy, ev, dt_min, gc и т.д. (вручную?)
2. Начальная инициализация t = 0, ∆t = DT

3. Вычисление p1 = f(p0, t, ∆t) // шаг
4. Вычисление pp1 = f(p0, t, ∆t/2) // первые пол-шага
5. Вычисление pp2 = f(pp1, t+∆t/2, ∆t/2) // вторые пол-шага
6. Вычисление erxy = 2 * |xy2-xy1| // оцениваем ошибку координат
7. Вычисление erv = 2 * |v2-v1| // оцениваем ошибку скорости
8. Если ∆t ≤ dt_min, перейти к пункту 10 // такой маленький шаг объявим точным
9. Если ∆t < max(erxy/exy, erv/ev), перейти к пункту 14 // если ошибка er грубее, чем e*∆t, уточним

10. Если |xy1| ≤ R, перейти к пункту 18 // врезались!
11. Вычислить p0 = p1 // полетели!
12. Вычислить t = t + ∆t
13. Если t + dt_min/2 > DT, остановиться // прилетели!

14. Вычислить ∆t = 0.9 * min (exy/erxy, ev/erv) // коррекция шага
15. Вычислить ∆t = max (∆t, dt_min) // ограничиваем шаг снизу
16. Вычислить ∆t = min (∆t, DT - t) // летим не дольше, чем нужно
17. Перейти к пункту 3

18. Здесь будет блок посадки

Всё супер, кроме расчёта нового dt, но про него я уже писал. Блок посадки встраивается в эту же схему почти без поблем, но об этом позже.

Вот и славно! Новую формулу для dt вставил, заодно добавил проверку на er≠0. Исправил ошибку с проверкой посадки. Одна лишь математика получилась раз в пять длинней стареньких «Лунолётов». Но МК-161 выдержит! :-) Главное, чтобы других ошибок в псевдокоде не было.

Переменные:
exy — желаемая ошибка координат на секунду полёта
ev — желаемая ошибка скорости на секунду полёта
dt_min — минимальный шаг по времени
R — радиус планеты
DT — время манёвра
∆t — шаг по времени
t — время от начала манёвра
p0 — вектор координаты и вектор скорости корабля
p1 — вектор координаты xy1 и вектор скорости v1 корабля на следующем шаге
pp1 — вектор координаты и вектор скорости корабля на половинном шаге
pp2 — вектор координаты xy2 и вектор скорости v2 корабля через два половинных шага
erxy — оценка ошибки координат корабля
erv — оценка ошибки скорости корабля
var1 — вспомогательная переменная

Алгоритм:
1. Начальная инициализация p0, DT, exy, ev, dt_min, gc и т.д. (вручную?)
2. Начальная инициализация t = 0, ∆t = DT

3. Вычисление p1 = f(p0, t, ∆t) // шаг
4. Вычисление pp1 = f(p0, t, ∆t/2) // первые пол-шага
5. Вычисление pp2 = f(pp1, t+∆t/2, ∆t/2) // вторые пол-шага
6. Вычисление erxy = 2 * |xy2-xy1| // оцениваем ошибку координат
7. Вычисление erv = 2 * |v2-v1| // оцениваем ошибку скорости
8. Если ∆t ≤ dt_min, перейти к пункту 10 // такой маленький шаг объявим точным
9. Если ∆t < max(erxy/exy, erv/ev), перейти к пункту 14 // если ошибка er грубее, чем e*∆t, уточним

10. Если y1 ≤ R, перейти к пункту 25 // врезались!
11. Вычислить p0 = p1 // полетели!
12. Вычислить t = t + ∆t
13. Если t + dt_min/2 > DT, остановиться // прилетели!

14. Если erxy≠0, перейти к пункту 18
15. Если erv=0, перейти к пункту 22
16. var1 = ev/erv
17. Перейти к пункту 21

18. Вычислить var1 = exy/erxy
19. Если erv=0, перейти к пункту 21
20. Вычислить var1 = min (var1, ev/erv)
21. Вычислить ∆t = 0.9 * var1 * ∆t^2 // коррекция шага
22. Вычислить ∆t = max (∆t, dt_min) // ограничиваем шаг снизу
23. Вычислить ∆t = min (∆t, DT - t) // летим не дольше, чем нужно
24. Перейти к пункту 3

25. Здесь будет блок посадки

Вообще говоря, значение с малой погрешностью является также значением с большой погрешностью. Обратное неверно. Например, число P=3.1416 это "пи" с точностью 5 знаков. Но также это "пи" с точностью три знака, т.к. 3.14

Ваше рассуждение отлично подходит для интегралов функций. Зачем брать грубое решение если есть более близкое к истине. Тоже верно для интегрирования с dt=DT/m и dt=DT/(2m) (постоянный шаг).

Но если вы адаптируете шаг на лету, то приходится брать менее точный вариант. Грубо говоря, он вносит меньше "шума" в дальнейшее решение. Это трудно объяснить. В общем, можно доказать что метод использующий более грубое приближение имеет более широкий "рабочий" диапазон - позже разваливается, более надёжен, хотя и даёт несколько меньшую точность. Там есть теоремы, которые говорят можно ли верить оценке точности. Для "грубого" их доказать легче. Извините за любительство, но это уже на пороге моих знаний.

В любом случае наша система так хороша, что можно брать любое решение. Я выбираю грубое чтобы не сбивать с толку тех кому метод потребуется в более сложных условиях.

Поддержку обещаю. Дело чести - получить на калькуляторах точные решения.

А в каком варианте с уточнённой математикой Вас прогон в моей программе интересует? Когда маневр вычисляется за 1 цикл (3 режим), или когда время маневра делится на отрезки равные двоичному логарифму времени в секундах, округлённому в меньшую сторону до целого (4 режим)?

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Тот, что максимально близок к физической реальности. Наверное четвёртый.

высота       в.скор        г.скор       угол  кг  сек
0            0             0             0    10  1
2.37630348   4.750758789   0             87   500 12
46.54724252  2.930772207   333.1039629   87   500 14
71.2499754   1.305940244   699.6348483   88   500 15
54.21931672  -2.423833907  1107.386899   88   250 8
38.385351    -1.088302786  1329.688752   88   250 9
44.30659262  3.020615185   1566.365673   89   100 4
58.43979912  4.15973399    1665.531307

Ну и так далее...

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Спасибо. В субботу полетаю по вашему маршруту.
Существует ли симулятор МК-152 на компьютере под Linux?

Существует: http://mk.semico.ru/vk6.htm

Исходники от AtH и ELF-файл Linux (i586, UTF-8) можно скачать отсюда: http://mk.semico.ru/prog.htm

Существует также несколько кроссплатформенных оболочек, посмотрите здесь:

http://pmk.arbinada.com/node/667
http://pmk.arbinada.com/node/683

http://mk-152.livejournal.com/31178.html
http://mk-152.livejournal.com/32581.html
http://mk-152.livejournal.com/32803.html

Кстати, ВК-6 компилируется под MacOS X без каких-либо специальных настроек. Пока так его и буду использовать. :-)

Ну и как полёт?
На сколько велики расхождения?
Я, кстати, пробовал менять в расчете деления интервалов времени основание логарифма с 2 на 1,26 (примерно втрое чаще) и высота в конце разгона была примерно метров на 100 больше.

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Я, к сожалению, ещё не пробовал. Сегодня вечером организую себе окно и полетаю.

100 метров - это много. Смотрели ли вы как быстро "сходится" ваше решение при уменьшении шага (в какое число раз уменьшается ошибка при уменьшении шага в два раза)?

Сколько интервалов в секунду успевает отсчитать калькулятор?

Точно отсчитать число интервалов в секунду я не пробовал, но точно больше 10 на 9750 и соответственно больше 80 на 9860. Кстати, вы можете скачать с этого сайта эмулятор fx-9860G SD (в разделе "Эмуляторы") и погонять там мою программу даже в русифицированном варианте.

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Вот полёт.

Схема использует постоянный расход топлива. Проверено, что формула Цолковского выполняется (в отсутствие тяготения Лунолёт-3.2 и Лунолёт-3.1 набирают одинаковые скорости). Печатается ускорение в конце манёвра (максимальное).

Разошлись и скорости и высота. Лучшее совпадение дала горизонтальная скорость.

Попробуйте, уменьшая шаг по времени, выяснить к каким величинам сходится ваше решение. Я хочу выяснить в чём причина расхождений, а с вашей помощью это сделать проще всего.

В принципе для проверки математики много команд вводить не надо. Достаточно примера одного манёвра, достаточно длительного и чтобы обе скорости были отделены от нуля. Взлёт под 45 градусов с большим расходом хорощо подходит.
Вот пример.
Здесь два полёта. Один разбит на десять манёвров, второй выполнен в один манёвр. Результаты совпадают с точностью до трёх знаков. Компьютер потратил на рассчёт каждого манёвра не более 20 итераций. Это я к тому, что калькулятору такое вполне под силу.

> два полёта. Один разбит на десять манёвров, второй выполнен в один манёвр. Результаты совпадают с точностью до трёх знаков.

Впечатляет! Точность, кстати, ПМКшная — все восемь знаков точны.

Кстати, предлагаю ввести графу "время полёта". Будет удобней ориентироваться в "простыне" результатов.

С простынёй я сейчас борюсь. И время полёта полёта тоже запланировано.

Попробовал уменьшить основание логарифма до 1.1 и получил в результате описанного взлёта такие параметры:
высота 99.84190091 м
верт. ск. 5,929834386 м/с
гор.ск. 1665.4996 м/с
пройденный путь 1,710084941 градуса 51873.41289 м

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Удивительно, что ошибки по горизонтали намного меньше, чем по вертикали. Правда ли что ваша пограмма считает вертикальную скорость с использованием горизонтальной скорости в конце манёвра?

То есть: правда ли что:

new_vx = f(y, vx, vy, dt, ax)
new_vy = g(y, new_vx, vy, dt, ay)

?

обратите внимание, что new_vx зависит только от значений в начале отрезка интегрирования, а new_vy зависит от new_vx - скорости в конце манёвра. Правда ли это?

Извините, что три раза один вопрос задал.

Ничего удивительного в этом нет. Постоянная тяга двигателя реализована математически точно, а вот плавное изменение центростремительного (и центробежного) ускорения и ускорения Кориолиса просчитать в один проход трудно. Не известны скорость координаты корабля в конце маневра. Поэтому, как и у Пухова, Кориолис и центростремительное ускорение на каждом элементарном интервале вычисляется в начале интервала и остаётся до конца интервала фиксированным. Разумеется ошибки по горизонтали - это ошибки в Кориолисе, а ошибки по вертикали - это ошибки центростремительного ускорения, которое значительно больше Кориолиса по модулю. Можно было бы вычислить маневр на каждом элементарном интервале путём нескольких последовательных приближений (т. е. вычислить "рывок" по Кориолису и по центростремительному), пока ошибка в координатах у двух последовательных расчётов не станет меньше заданной. Однако я это сделать поленился в своё время. Модель и так оказалась гораздо точнее Пуховской. А делить манёвр на 10000 элементарных интервалов для калькулятора (как сделали Вы для компьютера), я думаю, избыточно долго. Я именно поэтому делю маневр на интервалы в логарифмической зависимости, чтобы полёт по инерции на орбите не ел слишком много времени. Если логарифм по основанию два, например, то маневр в 15 секунд разделится на 4 интервала, а в 1000 секунд на 10 интервалов. Кстати, разница в высоте орбиты через пол витка при основании логарифма 1,26 примерно на 200 метров ниже, чем при основании 1,1 (то есть подъём с примерно 100 м до 4 км за 3200 секунд орбитального полёта по инерции). Сейчас у меня особо нет времени для разработки пятого варианта Лунолёта 3, но может я это как-нибудь сделаю.
А формулы, по которым я считал, можно посмотреть здесь: http://pmk.arbinada.com/node/214#comment-883 и далее по тексту

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

> А делить манёвр на 10000 элементарных интервалов для калькулятора (как сделали Вы для компьютера), я думаю, избыточно долго.

10000 - это предел числа итераций на один манёвр. При его привышении компьютер сообщает, что расчёт не укладывается в указанное количество. Это нужно, чтобы JavaScript не зацикливался с потерей данных пользователя. На деле шаг интегрирования занимает около секунды при работающем двигателе, около 10 секунд на эллиптической орбите и несколько десятков, если орбита близка к круговой.

В нормальной игре калькулятор и компьютер будет делать около 10 итераций на манёвр. Вообще число итераций пропорционально длительности манёвра.

Кстати, я тоже встроил счётчик итераций. Когда понял, что при ограничении шага снизу зацикливания не происходит, счётчик убрал.

Счётчик итераций всё равно нужен, так как у меня блок посадки встроен в основной цикл. Верхний уровень алгоритма такой: после очередного шага определяется состояние: полёт, заглубление или мы попали в тонкий приповерхностный слой. В последнем случае выдаётся сообщение о посадке, и манёвр прерывается (автоматика выключает двигатель при касании поверхности). Если произошло заглубление то шаг отменяется, DT обрезается (DT = <t после шага> - t0), dt делится пополам и попытки продолжаются. На деле это соответствует поиску корня уравнения (нулевой высоты) методом двоичного поиска. Ясно что в таком варианте вероятность превысить число итерации возрастает.

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

Сделал доделку, чтобы вычислять "рывок" по горизонтали и вертикали для кориолиса и центорбежной-центростремительной силы. для основания логарифма 1,1 получилось

высота       в.скор         г.скор          пройденный путь (рад) и (м)
103.4892952  6.075283793    1665.510818     1.7100913096(рад) 51873.60605(м)
110.2462834  6.306145104    1665.531254     1.710167907(рад)  51875.92956(м)

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Какой набор чисел к какому варианту относится?

Первая строка к логарифму 1.1, а вторая к логарифму 1,26 (первый вариант точнее, поскольку делит на более мелкие интервалы времени)
Кстати, нашёл у себя ошибку. При расчёте рывка забыл поделить на время. Исправленные данные:

выс           в. скор       г. скор       пройденный путь (рад)
107.7995089   6,269098967   1665.538966   1,7101666633     (log1.26T)
107.5520863   6,264084888   1665.537792   1.71010141       (log1.1T)

Возможно ещё линейный рывок заменить его квадратичной интерполяцией, но из 28 регистров уже 26 использованы. Можно ещё задействовать списки 6х255 переменных, но на мой взгляд - это уже слишком. Думаю на этом я остановлюсь.

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Что у вас за калькулятор? Где про него почитать? Число регистров? Быстродействие?

Как ни печально, но выловил у себя ещё одну ошибку. При вычислении рывка по кориолису забыл в одном месте (-) поставить и в результате рывок вычислялся не из разности а из суммы кориолисов.
итак для деления интервала времени T > 1 сек. на количество отрезков = log1.26(T)+1:

высота        в.скор        г.скор       пройденный путь
107.7888969   6.268390666   1665.49325   1.710147678(рад) = 51875.31594(м)

О моём калькуляторе:
casio cfx-9850GB Plus. Инструкция здесь: http://www.leningrad.su/calc/docs/rus_9850.rar
точность мантиссы - 15 десятичных цифр. Тактовая частота процессора - 20 МГц

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

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

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

28 переменных должно хватить для замены Эйлера (первый порядок) на модифицированный Эйлер (второй порядок). Это даст значительное улучшение точности.

Учитывая, что сила притяжения убывает пропорционально квадрату расстояния, второй порядок даст практически точное решение. Сегодня погонял калькулятор в транспорте, по дороге на работу. Кажется в программе больше нет ошибок. Но придётся оптимизировать весь код, поскольку для дополнительного отсчёта надо запомнить три параметра, а у меня осталось два регистра из 28 доступных (кстати у 7400 только 26 регистров доступно), или залезать в списки.

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Я ориентировался на старый третий "лунолёт". Там было занято 12 регистров. Некоторые использовались несколько раз, увеличим это число до 16. Остаётся в запасе ещё 10. Я не знаю вашей программы, где мои расчёты неверны?

Второй порядок, к сожалению не даст почти точного решения (его и четвёртый не выдаёт), но будет намного лучше первого.

При делении манёвра на фрагменты часть переменных дублируется, типа общее время маневра и текущее время и т.п.
Отдельно рассчитываются реактивные манёвры, отдельно гравитационно-кориолисовы, а затем результаты объединяются и по скоростям и по координатам (основываясь на линейности Ньютонова пространства), что увеличивает число независимых переменных. Некоторые регистры у меня так-же используются многократно. Ну можно ещё некоторые константы (уск. своб. пад., радиус Луны и т.п.) прописать явно, вместо записи в регистры. Но программа от этого удобнее не станет.

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Что такое списки?

Списки, в этой модели, - это одномерные массивы, которые могут содержать от 1 до 255 переменных (задаётся пользователем). В моделях поновее список может быть длиннее.

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Для оценки точности сделал два орбитальных облёта Луны. Один двумя командами по пол-витка (3250 секунд) и один - четырьмя командами по четверть витка (1650 секунд) разница высот периселения и апоселения примерно 4 км. Высота орбиты изменялась от витка к витку примерно на 4 метра в одной и той же точке орбиты. А как у Вас?

Мои программируемые калькуляторы:
Б3-21, Б3-34, МК-61, МК-52, МК-85
CASIO: cfx-9850GB+, fx-9750G+, fx-9750GII, fx-9860G, Algebra fx-2.0, fx-5800P, fx-7400G+
HP: 50G, 48G, 35s
TI: Nspire-CAS, Voyage-200, 89Titanium
SHARP EL-9600G

Вот мой вариант. Точность 10^-9 в секунду, время полёта вариант. Точность снижена до 10^-6, достаточно 24 итераций. Ошибка вползла в миллиметры (должна быть не более 10 мм).

И наконец такой вариант. Точность снижена до 10^-3, достаточно 7 итераций. Ошибка не более 10 метров.

А что уважаемое сообщество думают про идею обсчета не Пути к Земле, а Пути от Земли? Я имею в виду нынешнюю миссию Нью Хорайзонс http://en.wikipedia.org/wiki/New_Horizons а также предыдущие Вояджеры.

Например, известен тот факт, что Нью Хорайзонс никогда не догонит Вояджеры, но не потому, что они были запущены раньше, а потому, что они летят быстрее, так как они применяли гравитационный маневр для ускорения.
Вот, к примеру вопросы типичного домохозяина, к коим я причисляю и себя:
а) почему Нью Хорайзонс не применял гравитационный маневр для ускорения;
б) как вообще работает механизм гравитационного маневра для наращивания скорости?

---------------------------
Истина где-то рядом
aldebaran.ru/author/samurov_vitaliyi/kniga_dozvonitsya_do_devyi/

Чем сильней разгонишься, тем быстрей мимо Плутона пролетишь ;) Небольшой гравитационный маневр был у Юпитера, я пару лет назад изучал траекторию NH в Орбитере.
У Вояджеров была серия этих маневров - сильно повезло с расположением планет-гигантов на тот момент.

Вот тут есть анимация, хорошо иллюстрирующая откуда берется приращение скорости при грав.маневре:
http://ru.wikipedia.org/wiki/%D0%93%D1%80%D0%B0%D0%B2%D0%B8%D1%82%D0%B0%...

Хорошая статья в Вики, спасибо!

---------------------------
Истина где-то рядом
aldebaran.ru/author/samurov_vitaliyi/kniga_dozvonitsya_do_devyi/

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

Не думаю, что здесь возможности калькулятора играют особую роль. Скорее это то, что человек, начинающий интересоваться космонавтикой, может понять и переварить.

Когда только-только садишься за баранку космолёта, одномерной модели вполне хватает.

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

Нет, двумерная модель не слишком упрощена. Искуственный спутник планеты или луны вращается вокруг неё в одной плоскости. Правильно выбрав систему координат можно большинство орбитальных манёвров свести к плоской "калькуляторной" модели. Физическая достоверность при этом почти не пострадает. Перелёт Луна-Земля правильно описывается плоской моделью, если будет проходить в плоскости эклиптики (в плоскости вращения пары Луна-Земля). И даже межпланетные аппараты далеко от плоскости эклиптики не уходят (так требуется меньше топлива и времени на полёт).

Пилот освоивший "плоские" полёты не смутится при решении "пространственных" задач. Основные энергетические зависимости, скорости, траектории и манёвры будут ему знакомы по тренировкам в плоскости.

Например, стартуете с Байконура, как планируете попасть в плоскость орбиты Луны? Или находитесь на полюсе Луны, какую "плоскую" систему координат выберете, чтоб вернуться на Землю? Большинство пилотов, знакомых только с "плоскими" полетами, растеряется даже в выборе направления при старте - тут не два направления, а все 360 градусов.

З.Ы. Мой предыдущий пост навеян challenge-м, который еще в начале года выложил на форуме Орбитера один опытный пилот. Там даны первоначальные условия: корабль уже на орбите Луны, причем наклон этой орбиты почти 90 градусов к плоскости орбиты Луны вокруг Земли. И надо попасть на Марс. Топлива всего на 1671 м/с (!!!) плюс в маневровых двигателях еще запаса на 30-40 м/с. Так никто до сих пор не решил. Я наткнулся на эту задачку недели три назад, лучший мой результат около 1722 м/с, топлива все-равно чуть-чуть не хватает.

Рад увидеть здесь живого пилота Орбитера! Ваш опыт здесь очень пригодится. Сложно передать словами, насколько мы вам здесь рады.

Не подскажете, как в Орбитере вводится манёвр? Пилот задаёт расход топлива или ускорение — что из них постоянно? Относительно чего задаётся угол реактивной тяги?

Все выполняется в реальном времени. Ориентируешь корабль с помощью ориентационных двигателей (attitude thruster) в нужном направлении и даешь тягу. Сила тяги регулируется, задается в ньютонах. Дальше, в зависимости от текущего веса корабля и скорости истечения продуктов сгорания программа высчитывает ускорение и расход топлива. Максимальная тяга двигателей ограничена. Пилоту, в основном, приходится учитывать только сколько времени понадобится на маневр, и сколько он сожрет топлива.
У основных двигателей тяга направлена строго вперед, выхлопная струя - назад. На некоторых кораблях есть дополнительные двигатели вроде hover thrusters - тяга направлена вверх (удобны при взлете/посадке, особенно на спутниках без атмосферы), или retro thrusters - тяга назад (я их почти не использую, и особого смысла в них не вижу). Также небольшие коррекции можно выполнять слабенькими ориентационными двигателями - они умеют не только вертеть корабль, но и позволяют давать небольшую тягу по всем 6 направлениям (удобны при стыковках).

Сила тяги постоянна во времени (то есть во время манёвра ускорение растёт)?

Конечно, ускорение растет с уменьшением веса корабля в связи со сжиганием топлива. Вообще, физика в Орбитере очень точная, там даже кувырки гайки Джанибекова можно смоделировать, если закрутить корабль. Я как-то ради интереса написал к нему плагин, который используя NAIF-овскую библиотеку Spice добавляет в Орбитер объекты точно следующий расчетной НАСА-овской траектории для разных реальных миссий. Проверял на New Horizons и Cassini. И запускал параллельно кораблики, траектории которых считает Орбитер - корабли летели практически рядом. Четко было видно, где НАСА выполняло коррекции курса. А, например, при гравманевре New Horizons около Юпитера разница в положнение корабликов, после выхода из поля притяжения Юпитера была вообще смешная по космическим меркам - несколько сотен километров. Но вот при сравнении НАСА-вской траектории Кассини вблизи Сатурна и посчитаной Орбитером разница была уже больше. Подозреваю, что модель несфирического поля тяготения Сатурна (из-за его приплюснутости) в Орбитере не очень точная.

Манёвр, на языке «Лунолётов», это одна команда на двигатель. Например ∆m = 65 кг за ∆t = 3 секунды.

1. Стандартные «Лунолёты» обеспечивали постоянный расход горючего q = ∆m/∆t . Например в данном случае q = 25,667 кг/с в течении всех 3 секунд. При этом топливо сжигалось, масса корабля уменьшалась и реактивное ускорение росло. Всё, как вы говорите.

2. Некоторые авторы современных «Лунолётов» встроили блок с постоянным реактивным ускорением. В этом случае расход горючего q = M×a/c изменялся так, чтобы обеспечить постоянное реактивное ускорение корабля a, компенсируя потерю полной массы M.

Любопытно, какой вариант реализован в «Орбитере»? Ну и современных космических аппаратах, раз вы летали с ними наперегонки. :-) Правильно ли я понял, что заданная вами реактивная сила (в Ньютонах) F = c×q сохраняется в течении всех 3 секунд манёвра — следовательно в течении этих трёх секунд будет, как в «Технике-Молодёжи», постоянен расход ∆m/∆t и изменяться реактивное ускорение a = F/M.

Я налетал на Лунолетах в детстве на своем МК-52 не одну тысячу километров ;)
В Орбитере (да и в реальности) немного по-другому. Двигатель выдает определенную тягу, и т.к. ускорение зависит от массы корабля которая в течении маневра уменьшается, то и ускорение будет расти.
Орбитер берет ∆t исходя из вычислительных возможностей компа и заданной скорости моделирования. На большинстве современных компьютеров при расчете маневра в реальном времени ∆t выходит очень маленькое. При увеличении скорости моделирования до x100000 шаг увеличивается, точность моделирования уменьшается, и для компенсации ошибок приходится применяеть более точное интегрирование, вплоть до Рунге-Кутты 8-ого порядка.

> В Орбитере (да и в реальности) немного по-другому.

Так где отличия? Только в расчётной части, которую мы сейчас поправим.

Главное, что принцип управления (постоянный расхода топлива) тот же самый. Недаром консультантом был Глазков! :-)

Согласен, по-сути это одно и то же, только задается немного по другому. Зависимость простая: 65 ПП 3 С/П - это 65/3*3600 = 78кН тяги в течении 3 секунд ;)
Только, возможно, для реализма максимальную мощность двигателей стоит ограничить.

Нам надо принять решение о модели расхода топлива во время манёвра для "точных" "лунолётов". Варианты:

1. Постоянный расход, ускорение растущее по простой формуле.
2. Постоянное ускорение, меняющийся по формуле Циолковского расход.

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

Ссылки на оригинальные "лунолёты" не подходят - там ошибки ускорения или расхода так велики, что непонятно, какая модель расхода имелась в виду с самого начала.

3. Сделать оба варианта и переключаться между ними.

Против третьего варианта сложность для начинающих.

Сейчас делаем, как в журнале. Благо там и тестовые прогоны, и задания. Потом посмотрим, что делать с получившейся точностью и куда развивать. К тому же сейчас «Энергия» внедряет цифровые бортовые компьютеры, возможно и деньги считать научится.

Исходные «Лунолёты» имели постоянный расход. Это просто доказывается по ИПД. :-)

Журнальные задания вам не помогут. Разница между журнальным и "точным" "лунолётами" изрядно велика (у меня на сайте есть журнальные "лунолёты", полёты на которых совпадают с полётами на калькуляторах, так что знаю о чём говорю). Давайте пока делать постоянную тягу / меняющееся ускорение. А сравниться мы сможем только друг с другом и ещё с Орбитером и ещё проверить формулу Циолковского и ещё законы сохранения энергии и момента при безмоторном полёте.

Постоянная тяга ещё хороша тем, что формат входных данных совпадает с журнальным. А ИПД, извините, не аргумент :).

> Давайте пока делать постоянную тягу / меняющееся ускорение.

Постоянная тяга это постоянный расход? Так сейчас и реализовано. Ваш точный Лунолёт, к сожалению, основан на постоянном ускорении и для прямого сравнения с ЭКВМ придётся его немного обработать напильником.

> Постоянная тяга ещё хороша тем, что формат входных данных совпадает с журнальным.
> А ИПД, извините, не аргумент :).

Ну, если для вас точный расход заказанного топлива (ИПД) за точно заданное время (ИПС) не аргумент, то и формат входных данных можно поставить под сомнение. Вдруг хитрый Михаил Пухов подделывает реактивное ускорение так, чтобы тютелька-в-тютельку выполнить ваш манёвр ∆m/∆t — изменяя в процессе ускорение, скажем, по синусоиде? :-)

Вариант, обработанный напильником уже есть: полёт.

Про ИПД: Вот логическая цепочка: Я опирался на "Лунолёт-2" в котором формулы движения абсолютно точные, если принять постоянное ускорение. "Лунолёт-3" при небольших горизонтальных скоростях и выстоах должен быть близок ко второму. А это возможно только при постоянном ускорении. Тогда я решил, что ошибку в "точном" третьем "лунолёте" выгоднее перенести на расход топлива. Теперь мой третий "точный" "лунолёт" летает как второй из журнала при малых скоростях и высотах. Но это всё рассуждения почти бесплотные. Летать будем на постоянной тяге.

Здорово, про напильник в кустах! Я же сдуру на шагах 74-84 закодировал 90,349 ПП 5 С/П. Конечно, можно вручную нулями хвост забить, даже на ВК-6.

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

Судя по интерфейсу, программе, опубликованному алгоритму и опыту полётов замышлялся именно постоянный расход. Что же до точности, задания ЦУПа КЭИ можно попробовать перерешать на высокой точности и посмотреть, что и как выполнимо. Есть вероятность, что они давались там не от балды, просто с учётом достижимой тогда точности на ПМК.

Вы плясали от расхода всех "лунолётов", я от небесной механики второго "лунолёта". Топлива в "журнальном" третьем "лунолёте" уходило меньше - ускорение считалось для массы корабля в конце манёвра. А в первом и втором "лунолётах" масса бралась в начале манёвра. Первый и второй летали "тяжело", все остальные "легко". Точное решение где-то между ними.

На высокой точности все манёвры Коршунова выполнимы. Коршунов, честно говоря, не очень оптимально летал. Виток вокруг Луны делается с изрядным остатком топлива даже на "точном" "лунолёте", который тянет "хуже" журнального.

Ещё про рекордные полёты: на втором "лунолёте" выгоднее разгон выполнять небольшими порциями (следующая порция расходуется при меньшей массе корабля), а на третьем - большими (порция расходуется на более лёгком корабле). Такие дела.

Новый "точный" "лунолёт" тоже будет требовать разбиения манёвра на части для постановки рекордов - надо держать ускорение близким к максимальному, а оно растёт всё время.

Со всеми утверждениями согласен. Отмечу лишь, что зря вы упомянули журнальный вариант для унижения легендарного пилота. Коршунов вряд ли летал на журальном Лунолёте, т.к. жил в мире, где уже создана ЭКВМ. :-P Летал на журнальном Лунолёте Михаил Пухов, проверяющий отчёт Перепёлкина.

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

Коршунов летал на таком лунолёте:

"""
— Вот я и думаю, — сказал Михаил Коршунов. — Простите, Эдуард, если не ошибаюсь? Вы говорили, что много летали на ракетах. Не знаете, где можно раздобыть корабль, хотя бы плохонький?
Эдик поднял лицо. Краска с него уже схлынула, а в глазах появилось выражение, которое мне не очень понравилось. У него такое бывает. Что-то нехорошее, мстительное.
— Плохонький? — повторил он.
— Меня устроит любой.
— Тогда я вам помогу. У меня есть именно то, что вам нужно.
Так сказал Эдик. Я-то знал, что у него нет ничего, кроме старого лунолета, вроде того, на каком мы с сыном совершили свое невероятное путешествие.
"""

Рассказчк (Александр Перепёлкин) летал на лунолётё который по математике совпадает с "калькуляторным" первым и вторым. А Коршунов, стало быть летал на чём-то похожем :).

В реале, согласен. Запас топлива нужен больше, чем в математической модели. Но схожесть «Лунолётов» замечена Александром Перепёлкиным. На тот момент его стаж в штурманах составлял от силы пару минут и если в электронике он разбирался, то особенности космонавтики ему пришлось изучить позднее. :-)

Возможно, на «Кон-Тики» был похожий киберштурман, но его-то Лунный Коршун и выкинул. В «Истинной правде» же слишком много цифр, людям такое запомнить сложно. Скорее всего, компьютерщик Перепёлкин подготовился к встрече с Пуховым заранее, используя для подготовки своей презентации ту технику, которую сумел раздобыть в XX веке. Только из-за этого там цифры и сходятся.

Попробовал отлаживать. Идёт туго. Нужен останов после каждого шага dt. Нужны тесты в отсутствие тяготения только по вертикали и только по горизонтали. Нужен ввод команд с клавиатуры - тестировать так удобнее. Продолжаю попытки.

ВК-6 позволяет делать точки останова, даже запоминает их. Синтаксис, как у ДОСовского DEBUG'а. Справка есть в readme'шке.

Для шага по ∆t рекомендую точки останова в начале пунктов 8 (на PRM43) или 10 (на PRM42) — адреса 0231 и 0256 соответственно. Вообще, распечатки исходного кода и листинга очень помогают.

Ввод команд с клавиатуры можно легко встроить вместо прошитой команды, адреса 0074-0084. Места там более, чем достаточно. Ладно, мне надо идти на поминки двоюрдного брата. Люди ждут. :( Поздно вечером встретимся.

Баг: подпрограмма "физика" не знает текущего времени внутри манёвра. В результате ускорение постоянное.

Баг: с моей подачи формула неверна:
ay = ray + vx * vx / y - gc * y^2
надо так:
ay = ray + vx * vx / y - gc / y^2

Баг: условие посадки неверное.

Не исправили, к сожалению. Блок посадки не даёт взлететь. Первая итерация Эйлера не меняет высоту сразу после старта, только скорость. Блок посадки объявляет посадку.

Внёс обработку старта в блок посадки. Теперь он даёт добро на взлёт. Также заменил константы на ваши и ввёл предложенный вариант расчёта ошибок.

Выложено там же, проверяйте.

Обязательно. Только вечером.

Баг: в этих строках остались лишние умножения на 2:
+ FSQRT 2 * PM44 ; erxy = 2 * sqrt(..)
+ FSQRT 2 * PM45 ; erv = 2 * sqrt(..)

Давайте уточним работу с ошибками:

// относительные ошибки, безразмерные
rerxy = |xy2-xy1| / (exy * dt)
rerv = |v2-v1| / (evxy * dt)

// максимум и коэффициент (зависит от метода)
rer = 2 * max(rerxy, rerv)

// делаем шаг, если можно
if rer

Убрал множитель, изменил проверку на возможность шага и переставил вычисление относительных ошибок до неё.

Заодно сделал переключатель метода в R57. Пока должен работать лишь метод Эйлера (R57=0). Если оно по-прежнему работает, вставка Рунге-Кутта по R57=2 не составит труда.

Качать там же.

Я с опозданием, но...

Давайте уточним работу с ошибками:

// относительные ошибки, безразмерные
rerxy = |xy2-xy1| / (exy * dt)
rerv = |v2-v1| / (evxy * dt)

// максимум и коэффициент (зависит от метода)
rer = 2 * max(rerxy, rerv)

// делаем шаг, если можно
if rer
;6. Вычисление rerxy = |xy2-xy1| / (exy * ?t)

PRM36 PRM37 * PRM50 PRM51 * - Fx^2 ; (x1 * y1 - x2 * y2)^2
PRM37 PRM51 - Fx^2 ; (y1 - y2)^2
+ FSQRT ; erxy = sqrt(..)
PRM34 PRM54 * / ; rerxy = erxy / (exy * ?t)

;7. Вычисление rerv = |v2-v1| / (ev * ?t)

PRM38 PRM52 - Fx^2 ; (vx1 - vx2)^2
PRM39 PRM53 - Fx^2 ; (vy1 - vy2)^2
+ FSQRT ; erv = sqrt(..)
PRM35 PRM54 * / ; rerv = erv / (ev * ?t)

;7. Вычисление rer = 2 * max(rerxy, rerv)

Kmax 2 * PM55 ; rer = 2 * max(rerxy, rerv)

;8. Если ?t ? dt_min, перейти к пункту 10 // такой маленький шаг объявим точным

PRM43 PRM54 - ; dt_min - ?t
Px<0 Tochno

;9. Если ?t ? max(erxy/exy, erv/ev), перейти к пункту 14 // если ошибка er грубее, чем e*?t, уточним

PRM55 1 - ; rer > 1
Px<0 KorShag

;*** Точность нас устроила, шаг верен ***

;10. Если y1 ? R, перейти к пункту 21 // врезались!

Tochno:
PRM42 PRM37 - ; R - y1
Px<0 Posadka
Vzlet:

;11. Вычислить p0 = p1 // полетели!

PRM36 PM26 PRM37 PM27 PRM38 PM28 PRM39 PM29 ; p0 = p1

;12. Вычислить t = t + ?t

PRM19 PRM54 + PM19 ; t = t + ?t

;13. Если t + dt_min/2 > DT, остановиться // прилетели!

PRM18 PRM19 PRM43 2 / + - ; DT - (t + dt_min/2)
Px<0 KorShag
PRM16 PRM17 - PM16 ; m = m - mass
PGOTO Manevr ; "1": Манёвр завершён!

;17. Коррекция шага
KorShag:

; if (rer != 0)
; then new_dt = 0.9 * ?t / rer
; else new_dt = DT

PRM18 ; DT
PRM55 ; rer
Px!=0 LblZZ
0,9 PRM54 * PRM55 / ; 0.9 * ?t / rer
<->
LblZZ:
<-> PM54 ; ?t

Проверьте, что в заголовке нового «Лунолёта» стоит версия 0.03 — я только что исправил остроконечный баг с псевдокомандой .DD

Сейчас посмотрю ваше новое предложение.

У вас уже написан предложенный мной вариант. К сожалению, мне надо идти спать :(. До завтрашнего вечера.

Понял ваш вариант и учёл в новой выложенной версии. Вы хотите вообще избавиться от старых ошибок. Это можно. Но совсем избавиться от регистров R44 и R45 так не получится, т.к. до команды Kmax в вычислении rer первый результат будет выкинут из 4-х элементного стека, если за этим специально не следить.

Вообще, из-за роскошных возможностей МК-161 я отложил оптимизацию на потом. Её будет много. :-))) Сейчас при кодировании каждый пункт записывается отдельной фразой на языке МК-161 — хотя в большинстве случаев в стеке осталось что-то полезное от предыдущего пункта. Эта независимость от контекста позволяет легко менять пункты местами. Что, кстати, и позволило вам переставить относительные ошибки в другое место чуть ли не быстрее, чем это сделал автор программы.

Спасибо за математику, она у вас лучшая. Я же сейчас безработный, так что посижу ещё немножко с Рунге-Куттом. Вдруг заработает? :-)

Я пересмотрел свой вариант и не нашёл места где теряется значение. Я, правда, не отлаживал.

Часто формулы можно закодировать так, что предыдущее значение не потеряется — используя лишь 3 регистра стека. При оптимизации мы займёмся и этим тоже, сэкономим кучу регистров. Тогда это ваше умение будет востребовано, две пары глаз всегда лучше. Просто при написании больших программ для МК-161 (а я успешно писал программы, в два раза длиннее этого Лунолёта) такой подход вызывает проблемы для развития программы.

Если зависимость от предыдущего содержимого стека мешает свободно перемещать фразы ЯМК, то необходимость хранить в стеке предпоследние данные мешает свободно изменять фразы ЯМК. Таких вещей в момент активного развития программы лучше избегать — мы приближаемся к барьеру в сотню фраз и скоро запоминать, какие из них должны иметь глубину 2 или 3 начнёт становиться сложно.

Например, мне не удалось избежать стековой зависимости между первым двумя фразами блока "ФИЗИКА". Я решил сэкономить регистр под fm и специально отметил этот опасный момент в комментариях, чтобы кто-нибудь (например я сам, не выспавшись) не вставил что-нибудь между этими фразами. Скорее всего это был единственный такой момент — но даже я без текста не уверен в этом на 100%.

Практика программирования под ЭКВМ, с её роскошными 10000 шагами и 1000 регистрами, немножко осаждает лихость ПМКшников. Зато старый опыт помогает экономить байты и секунды там, где это действительно нужно. И когда это нужно.

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

Мы с вами написали почти одно и то же!

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

Залил программу в МК-161, запустил. Думает долго, ошибок не выдаёт. Через 10 минут стало интересно, что же она там делает? :-) Остановил. Счётчик команд 0424. РИП19 (t) выдал 1,09 — летает, но жутко медленно. Видимо, для Эйлера и неоптимизированной программы точность 1 ВП -09 немножко перебор. :-)

Смотрю на R54 (∆t), там 1 ВП -03. Другими словами, алгоритм почему-то идёт по самому маленькому шагу. Который, судя по черепашьей скорости, действительно мелкий. Остальные регистры:
R26 = 3,84460850472 ВП -05 (x)
R27 = 1738120,36345 (y)
R28 = 72,6620629757 (vx)
R29 = 120,903035194 (vy)

Не знаю уж, насколько точные значения. Но осмысленные. Запускаю ваш Лунолёт:
расстояние 66.831
высота 120.374
гор.скорость 72.662
верт.скорость 120.903

Некое сходство, вроде, есть. Запустил сначала, уже настроившись на 46-минутное ожидание окончания 5000 итераций. :-) Жаль, аккумулятор практически убит, более полугода пролежав в полностью разряженном состоянии.

Манёвр посчитан ровно за 50 минут, по 10 минут счёта на одну игровую секунду! Остановился на 0311, как и нужно. РИП19 (t) выдал ровно 5 — столько же, сколько и РИП18 (DT). РИП16 (m) выдал 100 (надо исправить). РИП54 (∆t) по-прежнему равен РИП43 (dt_min), то есть 1 ВП -03.

Остальные регистры:
R26 = 2,94630505258 ВП -04 (x)
R27 = 1738742,08358 (y)
R28 = 155,601135306 (vx)
R29 = 197,583262397 (vy)

Ваш вариант:
расстояние 512.122
высота 742.133
гор.скорость 155.602
верт.скорость 197.584

Не совсем, конечно, заказанная точность в 5 ВП -09… но уже есть надежда, что код имеет какое-то отношение к космическим путешествиям.  :-)) Ладно. Раз перешли на этот этап, добавил уменьшение массы топлива (без проверки) и рудиментарный интерфейс с подсказками для ввода команд «Лунолёта».

Получившуюся программу выложил на привычное место, но ещё не проверял на МК-161.

Он и будет ходить по минимальному шаги. А раз упёрлись в минимальный шаг, то ошибка выросла. Минимальный шаг - это средство от зацикливания при экстермальных входных данных.

Мы методом Эйлера (первого порядка) решаем диффур второго порядка (вторые производные). Я тестировал на точности 0.1. Видимо, 0.001 - наш предел по точности пока Эйлер работает. dt_min у меня был 1е-6.

Наша задача сейчас убедится, что метод устойчиво выбирает шаг dt и колебания dt невелики. А так же что ошибка не превосходит заданной оценки. Давайте тестировать короткие манёвры. Я пока только по вертикали полетал и без гравитации. 0/2/1 - вот такой старт делал.

Про обработку ошибок: этот кусок надо будет переделывать под методы высшего порядка. А он сейчас слишком сложен для переделки. Предлагаю такой вариант:

// относительные ошибки, безразмерные
rerxy = |xy2-xy1| / (exy * dt)
rerv = |v2-v1| / (evxy * dt)

// максимум и коэффициент (зависит от метода)
rer = 2 * max(rerxy, rerv)

if (rer != 0)
then new_dt = 0.9 * dt / rer
else new_dt = DT

А дальше ограничения сверху и снизу:

22. Вычислить ∆t = max (∆t, dt_min) // ограничиваем шаг снизу
23. Вычислить ∆t = min (∆t, DT - t) // летим не дольше, чем нужно

Математика та же. Только один условный переход, должно стать проще. И будет только одно место, где в будущем придётся править под Рунге-Кутту и Дорманда-Принса.

Она сходится! Я летал без гравитации, только по вертикали, 0,000001 dt_min, 0,001 exy ev, 0/2/1.
На остальные тесты нет времени. Вечером буду продолжать.
Готовтесь к Рунге-Кутте.

Можно продолжать.

А Рунге-Кутт нужен. Эйлер на МК-161 всё-таки тормозит. 10 минут на игровую секунду мне не нравится!

Держите ссылку на метод. Надеюсь, она вам знакома.
Использовать так: считаем для dt и для два раза по dt/2. Считать ошибку также. Коэффициент для ошибки и расчёт new_dt немного поменяются. Но это позже.
Пожалуйста, сделайте так чтобы методы можно было легко заменять. Например, можно сделать вызов процедуры, которая экстраполирует и считает ошибку (без коэффициента).

Предлагаю установить традицию — сначала написать и проверить псевдокод, а уже потом кодировать его в ЯМК и отлаживать полученную программу. Должно сэкономить время.

Вот Эйлер в наших обозначениях, из последнего текста программы:
1. Вычисление p1 = f(p0, t, ∆t) // шаг
2. Вычисление pp1 = f(p0, t, ∆t/2) // первые пол-шага
3. Вычисление pp2 = f(pp1, t+∆t/2, ∆t/2) // вторые пол-шага

После вычисления трёх точек ошибки оцениваются так:
4. Вычисление erxy = 2 * |pp2.xy2 - p1.xy1| // оцениваем ошибку координат
5. Вычисление erv = 2 * |pp2.v2 - p1.v1| // оцениваем ошибку скорости

Мы встраиваем Рунге-Кутта как альтернативу пунктам 1-3 или же это другой способ вычисления f(p, t, dt) ?

Страничку посмотрел. Чистая математика примерно ясна, но мне трудно установить соответствие этих x, y, y', f, h и k с нашими реалиями. pdf'ка по Рунге-Кутту, при внимательном рассмотрении, кроме опечаток основывается на зависимости a(v) и v(r). У нас же основным параметром, шаг по которому делается, является время.

Например, пусть y и x это наше p и t; h это наше dt, а их f() это наш блок "ФИЗИКА", который теперь обозначим буквой F. Также пусть мы повторно используем код пунктов 1-3 и встраиваем Рунге-Кутта туда, где раньше был вызов f(). Тогда пункты 1-3 выше вызывают вместо блока "ФИЗИКА" вот такую «прокладку»:
1. Вычисляем k1 = dt * F(p, t, dt)
2. Вычисляем k2 = dt * F(p + k1 / 2, t + dt / 2, dt / 2)
3. Вычисляем k3 = dt * F(p + k2 / 2, t + dt / 2, dt / 2)
4. Вычисляем k4 = dt * F(p + k3, t + dt, dt)
5. Вычисляем p_new = p + (k1 + 2 * k2 + 2 * k3 + k4) / 6

При таком подходе, кстати, нарушаются размерности. F() нам выдаёт четыре значения (координаты и скорости по обоим осям), мы их домножаем на секунды шага по времени и добавляем опять к координатам и скоростям.

В-общем, выбор метода в мой неопубликованный код встроил, вместе с заглушками для РК и ДП. А сейчас иду читать умные книжки и ждать вашего псевдокода.

Предложение принято.

В Эйлере у нас три раза вызывалась функция f. Эти три вызова и их аргументы сохранятся, но именится их содержимое. Итак, f состояла из двух этапов - вычисления производной и экстраполяции:
f(p, t, dt) {
p' = g(p, t)
return p + dt * p';
}

Базовым блоком для РК будет функция g:
f(p, t, dt) {
p1' = g(p , t ) // начало
p2' = g(p + dt/2 * p1', t + dt/2) // середина
p3' = g(p + dt/2 * p2', t + dt/2) // ещё раз
p4' = g(p + dt * p3', t + dt ) // конец шага!
// взвешенная сумма производных
return p + dt * (p1' + 2 * (p2' + p3') + p4') / 6;
}

Рунге-Кутт добавлен, Эйлер сохранён. Переключение по R57. Программу брать на прежнем месте.

Я провёл один тест, оба метода дали близкие результаты. На глаз особой разницы в точности или скорости не заметил. Но нужно лучше программу погонять, конечно.

Уже девятую страницу заполняем, приближаясь к 10% ресурса ЭКВМ. Кстати, я первый раз на МК-161 (да и вообще на ПМК) использую такое количество десятичных регистров. Причём здесь они индивидуально адресуемые, без косвенной адресации. Надеюсь, Дорманд-Принс всё-таки уложится в 100 регистров — тогда нам не придётся использовать трёхбайтовые команды доступа к ним.

Любопытно, сколько названий регистров способен запомнить мозг человека? Мы рискуем определить это на собственном опыте. :-) После оптимизации всё будет намного проще.

Думаю, что для ДП нужна будет косвенная адресация или кодогенерация. Если МК-161 поддерживает извлечение из регистра с декрементом его адреса, то мы напишем общий интегратор для табличных методов. Если нет, то код придётся генерировать программой на компьютере.

Поддерживает, конечно — всё, что умела ваша МК-52. Младшие регистры я специально оставил нетронутыми. Именно для таких случаев, ну и для циклов.

Табличка коэффициентов в Википедии выглядит не такой громоздкой. Да и вряд ли на машинке нам понадобится более сильный метод. Но напишем, конечно, ради вашей тайны золотого ключика ДП — раз уже написали Эйлера и Рунге-Кутта.

(перечитав Википедию) Ого! В Дорманд-Принсе мы можем оставить более точный прогноз пятого порядка. Более того, метод специально заточен именно на такой подход.

Правильно я понял, что регистрами с 0 по 3 можно идексировать регистры с большими номерами? И при этом будет автодекремент?

Более сильный нам не понадобится. Чем выше порядок метода, тем выше требования к разрядности компьютера. Большое число операций с разными по модулю операндами (чем больше порядок, тем больше разница в слагаемых) накапливают ошибки округления. И в результате точность тонет в оруглениях без следа.

Да прогноз ДП будет пятого порядка, но сам метод - четвёртого. Как РК. Надо будет сравнить, у кого шаги будут больше.

Поняли правильно, автодекремент соответствует МК-52. Более того, если к любому из них обращаться с помощью префикса P, будет эффект, как КИП↑ у Б3-34. Только этих бесценных регистров всего четыре, а нам могут понадобиться циклы FL0/FL3.

Кстати, к регистрам с номером больше 99 экономнее всего обращаться именно с помощью косвенной адресации. Часто, при оптимизации по размеру кода выгоднее обращаться к регистру функций исключительно с помощью косвенной адресации — даже с учётом занесения его номера в младший регистр.

Правильно ли я понял, что для ДП я могу принять coeff = 16/15, а для вычисления нового ∆t использовать k=4 ?

А что мешает использовать регистры R4-R6 с автоинкрементом? При этом R0-R3 освобождаются для циклов. Мне кажется, не имеет значения, в каком порядке будут расположены переменные в памяти.

Механизм предусматривался. Если чисел много и место важнее, чем скорость выполнения программы, то для упаковки и распаковки можно использовать байтовый буфер. Например, для ввода из таблицы достаточно побайтно считать упакованное число из памяти программ (R9044), записать его в буфер (R9034) и перенести число в RX из R9037 (по указателю R9030).

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

Начнём всё равно с регистрового файла, а там посмотрим.
За советы спасибо, учту на будущее.

Да, вы поняли правильно.

AtH, дабы не испытывать возможности человеческого мозга, можно попробовать заглянуть в документацию. ;)

В операторах записи и чтения, в том числе косвенных, допускается использовать не только номера регистров, но и идентификаторы, определяемые псевдооператором ".EQU".

К примеру, в любом месте программы (обычно в начале) записывается:

ay .EQU 15 ; R15 ay — проекция полного ускорения на местную вертикаль "ФИЗИКА"

далее в тексте программы:
P M ay (и т.п.) - записать в регистр,
P RM ay (и т.п.) - считать из регистра.

Пробелы и прописные буквы можно расставлять по вкусу, запись pmay или PrMay компилятор тоже поймёт.

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

Михаил Борисович — рад видеть, что вы читаете наше реалити-шоу. :-)

Про .EQU знаю, но стараюсь обходиться без него, пока это возможно. У этого есть три причины.

Во-первых, от ручного распределения регистров она не избавляет и документацию по регистрам всё равно приходится поддерживать в начале программы. Она с одной стороны большая, и удваивать её блоком .EQU не хочется. С другой стороны, порой регистр заковыристый и описать его одной строчкой не так уж и просто. Занимать ценное место таким избыточным кодом не хочется. Но это возможно решить, например, встраиванием такого синтаксиса в компилятор:

.R<номер> <id> ; произвольный комментарий

Вторая причина сложнее. Сейчас ошибку допустить проще, чем без .EQU . Но если она допущена, её просто отследить глазом, сравнив с табличкой регистров в начале текста программы. Более того, по номеру регистра его найти в списке проще, чем по алфавиту. Особенно если стараешься группировать регистры осмысленно, блоками — это будет в релизе «Лунолёта».

Третья же, главная причина в том, что идентификаторы не доходят до ПМК. И в момент отладки исходный код, где явно указаны регистры, несравнимо полезней — именно для разгрузки мозга и экономии времени. Под рукой готовые и комментированные блоки макросов, которые нужно лишь набить в режиме F АВТ, чтобы получить отладочную информацию.

Когда я использовал .EQU, код в режиме F ПРГ был совершенно незнакомый. Ассоциировать его с местом в исходном тексте было значительно сложнее. Даже автор с таким справляется плохо, а отладка превращается в кошмар. Когда я впервые с этим столкнулся, вручную перебил все .EQU обратно на номера регистров и перешёл на тот стиль, который вы наблюдаете. Создавать сложные программы стало если не просто, то несравненно удобней.

Но это хорошо, что есть куда переходить, когда будет сложно поддерживать номера регистров вручную. Пока что справляюсь неплохо. А там выкинем Эйлера и Рунге-Кутта, оптимизируем расчёты и регистровый бум закончится.

Страницы