Создание собственного псевдотрехмерного пространства для моделирования задач небесной механики и релятивистской физики с применением пакета OpenGL

Орлов Иван, Зисман Михаил,

г. Челябинск, лицей №37, 11 класс

Введение

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

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

Трудности, связанные с разработкой подобных программ заключались в отсутствии сведений об использовании системной библиотеки OpenGL в среде программирования Delphi. Зачастую OpenGL использовался в игровой индустрии и в профессиональных трехмерных редакторах, которые нуждались в больших системных возможностях. В настоящей работе нам удалось разработать программу, использующую собственное псевдотрехмерное пространство с применением OpenGL, не предъявляющую высоких требований к ПК.

Использованные ресурсы

В качестве языка программирования мы выбрали Pascal в среде Delphi 6. Для вывода графического изображения мы выбрали библиотеку OpenGL, которая более пригодна для создания трёхмерного изображения, чем низкоуровневые функции вывода Windows – функции GDI (Graphic Device Interface).

НЕ нашли? Не то? Что вы ищете?

OpenGL – это стандартная библиотека для всех 32-разрядных операционных систем, в том числе и для системы Windows [1].

OpenGL – не отдельная программа, а часть операционной системы. Это означает, что откомпилированное приложение, использующее OpenGL, применяет только стандартные модули, содержащиеся на любом компьютере с установленной операционной системой Windows 95/NT версии OSR2 и выше.

В настоящей работе мы использовали стандартные библиотеки opengl32.dll и glu32.dll находящиеся в подкаталоге \SYSTEM\ операционной системы Windows. Библиотека opengl32.dll отвечает за осуществление взаимодействия программы с акселератором или программной эмуляцией ускорителя за счет центрального процессора.

Представленная программа разрабатывалась на компьютере AMD Duron 900 Mhz. Минимальные требования к системе, на которой будет осуществляться запуск программы: Pentium 100, Windows 95 OSR 2, 16 MB RAM, клавиатура, мышь, рекомендуется установленный 3ds Max v4-5.

Основные принципы построения программы

Рассмотрим основной алгоритм работы нашей программы. Он состоит из следующих этапов:

Определение формата пикселя и установка контекста воспроизведения. Процедура записи параметров системы тел для каждого кадра в файл - скрипт Процедура чтения данных из файла-скрипта и прорисовка сцены. Процедура прорисовки траекторий движения тел.

Теперь рассмотрим эти пункты подробнее.

Прежде чем получить контекст воспроизведения, сервер OpenGL должен получить характеристики формата пикселя. Эти характеристики хранятся в специальной структуре типа TpixelFormatDescriptor [1]. Этот формат определяет конфигурацию буфера цвета и вспомогательных буферов. Для определения сервером OpenGL ссылки на контекст воспроизведения, используется переменная типа HGLRC (Handle OpenGL Rendering Context). Значение этой переменной автоматически определяется при запуске программы.

Заполнение полей структуры TPixelFormatDescriptor осуществляется в процедуре SetDCPixelFormat [1]. Функция ChoosePixelFormat осуществляет запрос системе, поддерживается ли на данном устройстве вывода выбранный формат пикселя. И вызовом функции SetPixelFormat устанавливается формат пикселя в контексте устройства. Наши установки корректируются сервером OpenGL применительно к реальным характеристикам системы.

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

Файл-скрипт – это файл, содержащий информацию о положении тел в каждый момент времени. Запись этой информации осуществляется расчетной программой Trajectories (а для решения задачи специальной теории относительности – RelatSTO).

Чтение из файла-скрипта производится последовательно по кадрам. Другими словами, положение тела в следующий момент времени заменяет его предыдущее положение. Каждое тело рисуется строго в начале системы координат OpenGL. Затем начало системы координат смещается в точку положения тела, для этого используется команда glTranslatef. В качестве аргумента вектора смещения выступают координаты положения тела. Для сохранения текущего состояния системы координат, перед прорисовкой следующего по порядку тела, используется команда glPushMatrix. После прорисовки тела, систему координат следует восстановить, для этого вызывается команда glPopMatrix [1]. Для прорисовки последующих тел вышеперечисленные действия повторяются.

В течение просмотра анимации, пользователь имеет возможность изменять положение “камеры” (“камера” – точка, из которой ведется просмотр сцены). Для этого мы используем команду glRotate для поворота сцены, и команду glScale – для ее масштабирования. glClearColor – отвечает за изменение фона сцены. В программе есть опция использования графического редактора 3ds Max, из которого также можно просматривать анимацию. Преимущество 3ds Max заключается в более качественной прорисовке тел. Недостаток в том, что создается файл *.avi, в котором невозможно изменить положение камеры.

Прорисовка траектории ведется для каждого тела отдельно (последовательно). Для установки непрерывности линии устанавливаем параметр GL_LINE_STRIP для текущей командной скобки, в которой идет чтение положений тел.

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

Рассмотрим некоторые примеры применения программы для задачи n тел.

1) Для начала рассмотрим систему четырех тел, расположенных в вершинах квадрата. Если массы всех тел будут равны, то, в случае нулевых начальных скоростей, результирующая сила для каждого тела будет направлена по диагонали квадрата, и их траектории движения сойдутся в центре, где расположен центр масс этой системы. Если массы будут неравными, то траектория тел сместится от центра квадрата. Теперь зададим начальные скорости всех тел, причем векторы скоростей направим под одним углом к сторонам квадрата (модули скоростей равны): тело №1 – (0;0;0) Vy=–10 м/с, m=1,5*1017кг; №2 – (105;0;0) Vx=–10 м/с, m=1,5*1017кг; №3 – (105;-105;0) Vy=10 м/с, m=1,5*1017кг; №4 – (0; –1e5;0) Vx=10 м/с, m=1,5*1017кг (параметры, стоящие в скобках, означают координаты тела в начальный момент времени; в дальнейшем запись вида 1,5*1017 будет заменяться эквивалентной записью вида 1,5e17). При этих параметрах получатся следующие траектории (см. приложение, рис.1).

Казалось бы: мы задали одинаковые расстояния и скорости тел, а, в конце концов, система распадается. Данный эксперимент показал, что, несмотря на симметричность системы, она является нестабильной.

2) Теперь проверим на стабильность модель куба. Тела верхней грани куба пусть движутся «по часовой стрелке», а нижней – наоборот. Размеры квадратов и скорости тел пусть будут теми же, что и в предыдущей модели. На рисунке видно, что такая система тел распадается гораздо раньше, чем предыдущая, траектории движения тел более сложные (см. приложение, рис.2).

3) В ходе исследований мы убедились в существовании относительно стабильных систем тел (двойных звезд). Рассмотрим одну из таких систем. В центре масс мы расположили еще одно тело. Его предназначение – показать, что объект, расположенный в центре масс такой системы, тоже обладает относительной стабильностью, при симметричности модели. Параметры системы: №1 (центральное) – (0;0;0) Vxyz=0 m=1e10кг; №2 – (100;100;0) Vy=–5e-2 м/с, m=1e7кг; №3 – (-100;-100;0) Vy=5e-2 м/с, m=1e7кг (см. приложение, рис.3).

Для окончательной проверки всей программы решено было провести серию «контрольных» тестов на известных примерах, с целью выявления ошибок и погрешностей. Для тестов была выбрана наша Солнечная система. Рассмотрим некоторые из «контрольных» тестов:

4) В первом тесте нашей задачей было получить траектории движения нижних планет Солнечной системы. Все параметры планет задавались в соответствии с реальными, за исключением их размеров. При реальных радиусах планет, их невозможно будет увидеть, рассматривая всю систему целиком. Результат расчета траекторий получился близким к реальному. При шаге, равном 40000 секунд, относительная погрешность составила ≈2,24% (см. приложение, рис.4). На рисунке видно отклонение орбиты Меркурия от окружности (в действительности, эксцентриситет орбиты Меркурия равен 0,21).

5) Чтобы окончательно убедиться в правильности построения математической модели программы, мы проверили траекторию движения Луны вокруг Земли, относительно Солнца. Результат отличался от реального на ≈4,3% (см. приложение, рис.5). На рисунке видно, что траектория Луны представляет собой траекторию знакопеременной кривизны (в реальности это действительно так).

В приложении (рис.6) приведены наиболее часто встречающиеся траектории в задаче двух тел.

Теперь рассмотрим некоторые тесты программы для задач СТО.

1) В первом тесте мы решили проверить, будет ли отклоняться частица, летящая с субсветовой скоростью, в гравитационном и электростатическом поле. Параметры системы тел: частица – (1e7;1e5) m=1,6e-27 Vx=-2e8 Vy=0 q=1,67e-19 центральное тело – m0=2e30, Q=5 при шаге в 1e-5 секунд и времени просмотра в 0,1 сек. Получаем (см. приложение рис.7). В данном примере мы рассмотрели модель отклонения заряженной частицы вблизи звезды, подобной Солнцу, но саму частицу расположили в 1000 раз ближе, чем радиус Солнца, чтобы показать отклонение (если частица будет пролетать на расстоянии радиуса Солнца, то отклонение будет в 1,75 секунды). Отклонение же в нашем случае равно ≈2,860.

2) Во втором тесте скорость частицы много меньше скорости света. Параметры системы тел: частица – (-1,5e8;0) m=1,6e-27 Vx=0 Vy=3e6 q=3e-16 центральное тело – m0=2e30, Q=1 при шаге в 3e-1 секунд и времени просмотра в 8000 сек. получаем спираль (см. приложение рис.8).

3) В третьем тесте увеличим массу центрального тела: частица – (-1,5e8;0) m=1,6e-27 Vx=0 Vy=3e6 q=3e-16 центральное тело – m0=1,1e32, Q=1 при шаге в 3e-2 секунд (см. приложение рис.9).

4) В четвертом тесте увеличим заряд центрального тела и уменьшим заряд у частицы: частица – (-1,5e8;0) m=1,6e-27 Vx=0 Vy=3e6 q=2e-16 центральное тело – m0=2e30, Q=4 при шаге в 3e-1 секунд (см. приложение рис.10).

Рисунки для задачи n тел в небесной механике (см. приложение, рис1 – рис6) выполнены с помощью программы 3ds Max, поскольку в ней тела показываются сплошными, и возможна одновременная прорисовка тел и их траекторий движения.

Погрешность, вносимая компьютером, составляет около 10-11%. Это происходит из-за невозможности учета всех знаков числа после запятой. По сравнению с ошибкой, вносимой шагом, ошибкой компьютера можно пренебречь. Ошибка, внесенная шагом (шаг равен 78840 секунд) в системе Солнце – Земля – Луна составляет ≈4,4%. Для исследования изменения ошибки, мы уменьшили шаг в 2 раза. При этом погрешность стала равна ≈2,2%. Она никогда не станет равной нулю, потому что для этого нам придется задавать шаг, стремящийся к нулю.

Основные выводы:

представленные в работе программы используют пакет OpenGL для создания псевдотрехмерного пространства, не нуждаясь в больших системных возможностях; с помощью представленных программ можно моделировать задачи релятивистской физики и n тел в небесной механике.

Литература:

OpenGL. Графика в проектах Delphi. – СПб.: БХВ-Петербург, 2002. – 352с. , Лифшиц поля. М.: Наука, 1988. – 510с.

Приложение

Математическая модель программы

Описанные нами программы мы решили использовать для решения задачи n тел в небесной механике и задач СТО. Для решения задачи n тел в небесной механике необходимо просчитывать положения всех тел для каждого кадра прорисовки. Для этого из второго закона Ньютона (ma=F) и формулы гравитационного взаимодействия F=GmM/r2 распишем проекции ускорения на оси x, y и z: mAx=–GmM*x/r3 , mAy =–GmM*y/r3 , mAz =–GmM*z/r3, где A – ускорение, x, y, z – проекции радиус-вектора на соответствующие оси, r – расстояние между телами. Так как ускорение есть производная скорости по времени, то заменим A на V′, где V′≈(Vi+1–Vi)/h, h – шаг по шкале времени, Vi – значение скорости в i-тый момент времени. Таким образом, получим окончательные уравнения для расчета скорости и координаты каждого тела в системе, состоящей из n тел:

Vi+1= Vi –hGM1(x-x1)/((x-x1)2+(y-y1)2+(z-z1)2)3/2 –hGM2(x-x2)/((x-x2)2+(y-y2)2+(z-z2)2)3/2–…–hGMn-1(x-xn-1)/ ((x-xn-1)2+(y-yn-1)2+(z-zn-1)2)3/2

Xi+1=Xi+hVi

Vi+1= Vi –hGM1(y-y1)/((x-x1)2+(y-y1)2+(z-z1)2)3/2 –hGM2(y-y2)/((x-x2)2+(y-y2)2+(z-z2)2)3/2–…–hGMn-1(y-yn-1)/ ((x-xn-1)2+(y-yn-1)2+(z-zn-1)2)3/2

Yi+1=Yi+hVi

Vi+1= Vi –hGM1(z-z1)/((x-x1)2+(y-y1)2+(z-z1)2)3/2 –hGM2(z-z2)/((x-x2)2+(y-y2)2+(z-z2)2)3/2–…–hGMn-1(z-zn-1)/ ((x-xn-1)2+(y-yn-1)2+(z-zn-1)2)3/2

Zi+1=Zi+hVi

Теперь получим уравнения расчета скорости и координаты для релятивистской частицы, движущейся в гравитационном и электростатическом полях [2]. Так как производная импульса по времени равна силе, действующей на тело, то dp/dt=F. Расписывая импульс и силу, получаем d(mV/(1-V2/c2)1/2)/dt=–(kQq+Gm0m)r)/r3(жирный шрифт означает векторную величину), где m и m0 – масса частицы и центрального тела, V – скорость частицы, Q и q – заряды центрального тела и частицы, r – расстояние от центрального тела до частицы, с – скорость света в вакууме. Если разложить расстояние и скорость на их проекции, то мы получим следующее: r2=x2+y2, V2=Vx2+ Vy2, или V2=(dx/xt)2+(dy/dt)2. Заменим dx/dt и dy/dt на x' и y' (первые производные от координат). То есть первоначальная система для получения основных формул имеет вид

d(mx'/(1-(x'2+y'2)/c2)1/2)/dt=–(kQq+Gm0m)x/(x2+y2)3/2

d(my'/(1-(x'2+y'2)/c2)1/2)/dt=–(kQq+Gm0m)y/(x2+y2)3/2

Проведя необходимые преобразования, в левой части уравнений получим выражения cm(x''(c2-y'2)+x'y'y'')/(c2-x'2-y'2)3/2, и cm(y''(c2-x'2)+y'x'x'')/(c2-x'2-y'2)3/2 ,где x'' и y'' – вторые производные координат по времени. Выразим из левой части первого уравнения x'', подставим его во второе и проведем преобразования. Получим окончательную систему уравнений для расчета скорости и координаты в следующий момент времени:

Vxi+1=Vxi–h*(kQq+Gm0m)(c2-Vxi2-Vyi2)3/2(x(c2-Vyi2)-xyVyi)/cm(x2+y2)3/2((c2-Vxi2)(c2-Vyi2)–xVyi2Vxi)

Vyi+1=Vyi–h*(kQq+Gm0m)(c2-Vxi2-Vyi2)3/2(y(c2-Vyi2)-xyVxi)/cm(x2+y2)3/2((c2-Vxi2)(c2-Vyi2)–yVxi2Vyi)

Xi+1=Xi+h*Vxi

Yi+1=Yi+h*Vyi

Задачу для релятивистской частицы мы решили рассматривать в плоскости, а не в объеме, руководствуясь законом сохранения момента импульса и сферической симметрией центрального тела [2].

Битовые флаги, используемые в нашей программе:

PSO (PFD_SUPPORT_OPENGL) – сообщает о том, что система поддерживает OpenGL. PDW (PFD_DRAW_TO_WINDOW) – осуществляет вывод в окно. PDB (PFD_DOUBLEBUFFER) – включает режим двойной буферизации, когда вывод осуществляется не на экран, а в память, затем содержимое буфера выводится на экран. Это очень полезный режим: если в любом примере на анимацию убрать режим двойной буферизации, то при выводе кадра будет заметно мерцание. PSG (PFD_SUPPORT_GDI) – включает опцию использования функции GDI. PGA (PFD_GENERIC_ACCELERATED) – устанавливает опцию использования графического акселератора (если компьютер оснащен таковым) [1].

СЛОВАРЬ

Контекст воспроизведения – структура, определяющая комплект графических объектов, связанных с ними атрибутов и графические режимы, влияющие на вывод.

Сервер OpenGL – режим работы библиотеки, при котором файл OpenGL32.dll обрабатывает команды, осуществляя вывод на устройство. При этом все приложения обслуживает одна копия сервера в памяти.

Формат пикселя – детальные характеристики используемого оборудования.