Размещено на http://www.allbest.ru/

РОССИЙСКАЯ АКАДЕМИЯ НАУК

СИБИРСКОЕ ОТДЕЛЕНИЕ

ИНСТИТУТ КРИОСФЕРЫ ЗЕМЛИ

УДК 519.6+521.1

осесимметричные многослойные вращающиеся структуры

И.И. Смульский

Тюмень 2013

Аннотация

Рассматриваются точные решения задачи ньютоновского гравитационного взаимодействия N материальных точек, находящихся на N2 концентрических окружностях. На каждой окружности осесимметрично расположено N3 тел с одинаковыми массами. Вся структура как единое целое вращается вокруг оси симметрии. Введена классификация параметров задачи и унифицированы их обозначения. Определена сила воздействия всех тел структуры на любое ее тело. Записаны дифференциальные уравнения их движения, которые сведены к линейным алгебраическим уравнениям для масс тел. Разработан алгоритм и программа RtCrcSt2.for для их численного решения. Кроме того, программа позволяет подобрать необходимые параметры структуры и подготовить файл начальных условий для исследования ее динамики с помощью системы Galactica.

Рассчитаны структуры с числом слоев N2 = 1, 2, 3, 4, 5, 15, 30, 100, 103, 1000 и числом тел N3 = 2, 5, 8, 10, 29, 30, 999. Приведены рисунки конфигураций некоторых из них и описаны их свойства. При малом количестве тел такие структуры известны в литературе как конфигурации гомографической динамики. Они представляют собой вложенные друг в друга воображаемые многоугольники, в вершинах которых находятся тела. Осесимметричность некоторых из конфигураций, например, вложенные друг в друга ромбы, или треугольник, вложенный в шестиугольник, не является очевидной. В работе рассчитаны также эти структуры и представлены на рисунках. В п. «Дополнительные уточнения» рассмотрены другие примеры решения этой задачи. Отмечен первоначальный метод ее решения при сведении задачи к системе нелинейных уравнений. Рассмотрены проблемы устойчивости структур и представлены примеры их изучения с помощью программы Galactica. В Заключении обсуждаются вопросы использования рассматриваемых структур в задачах космической и небесной механики.

1. Введение

Существуют точные решения гравитационного ньютоновского взаимодействия N-материальных точек при определенной их конфигурации. Например, при осесимметричном расположении N-тел на окружности эта задача решена в полном объеме [1] - [2]. При этом в центре окружности может находиться центральное тело. В зависимости от начальной скорости периферийные тела могут двигаться по эллипсу, параболе или гиперболе. Кроме такой однослойной конфигурации в ряде работ [3] - [7] рассматриваются точные решения конфигураций из нескольких слоев. Обычно их рассматривают как систему вложенных друг в друга воображаемых многоугольников, вращающихся с угловой скоростью как единое целое. В вершинах многоугольников расположены материальные точки, которые взаимодействуют между собой по закону тяготения Ньютона. В обобщающей работе [3] приведены примеры решения задач для вложенных друг в друга треугольников, ромбов, квадратов, пятиугольников и шестиугольников. При этом вершины соседних многоугольников могут лежать либо на одном радиусе, либо на радиусах, проходящих через середины сторон соседних многоугольников.

Эти задачи решены [3] для нескольких слоев таких многоугольников: для треугольников - до 4-х слоев, для квадратов - до 3-x, для пятиугольников и шестиугольников - до 2-х слоев. Наибольшее число взаимодействующих тел, находящихся в вершинах этих многоугольников, было 12.

В рамках гамильтоновой динамики эта задача сводится к системам алгебраических уравнений, которые решаются методами компьютерной алгебры [3], [8]. В каждом случае задача требует особого рассмотрения, и необходимо затратить немало усилий для ее решения.

Эти задачи взаимодействующих материальных точек, находящихся в вершинах вложенных друг в друга многоугольников, рассматривают как задачи гомографической динамики [3]. В свою очередь гомографическая динамика является новым разделом космической динамики. На наш взгляд, рассматриваемые задачи являются необходимой ступенью, ведущей к решению проблем существования колец планет, шаровых звездных скоплений и галактик.

В настоящей работе рассматривается несколько иной подход к решению этих задач. Он является продолжением метода, использованного для решения осесимметричной задачи [1]- [2] и при рассмотрении многослойных кольцевых структур [9]. Движение тел рассматривается на основе сил взаимодействия между ними. Вместо многоугольников, воображаемые стороны которых соединяют тела, мы рассматриваем окружности, на которых находятся тела. Для решения систем алгебраических уравнений используются численные методы.

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

Рассмотрим (см. рис. 1) многослойную осесимметричную структуру материальных точек, взаимодействующих по закону тяготения Ньютона. Она состоит из N2 окружностей, на каждой из которых расположено N3 тел. Совокупность тел, центры которых расположены на одной окружности, будем называть кольцом тел или слоем. Номера колец обозначим как j = 1, 2 … N2, а номера тел на каждом кольце - l = 1, 2 … N3. В плоскости xoyo, в которой тела располагаются, обозначим rj,l и j,l - полярные радиус и угол тела с массой mj,l. С целью упрощения в дальнейшем символом mj,l будем также обозначать и само тело. Все тела на одном кольце имеют одинаковый радиус rj,l = rj, где rj - радиус кольца, и массы их одинаковы, т.е. mj,l = mj. Угол первого тела на каждом кольце j,1 определяет вид структуры. В дальнейшем он будет задаваться. А углы остальных периферийных тел определяются по формуле

j,l = j,1 + (l - 1)0, (1)

где 0 = 2/N3 - угол между телами на кольце.

Итак, геометрия осесимметричной многослойной структуры определяется количеством колец N2, количеством тел на каждом кольце N3, радиусами колец rj и углами положения первых тел j,1. Масса каждого тела на кольце j равна mj и при наличии центрального тела с массой m0 масса всей системы будет

. (2)



Рис. 1. Геометрические характеристики осесимметричной многослойной структуры с параметрами: N2 = 5; N3 = 8; углы j,1 первых тел на соседних кольцах чередуются; исходная масса mi центрального тела и всех тел первого кольца немного превышает массу Солнца; линиями возле тел показаны вектора скорости; радиусы изображений тел пропорциональны их массам.

гравитационный многослойный ньютоновский взаимодействие

Вся система вращается с угловой скоростью . При задании параметров многослойной структуры: N2, N3, j,1, m0, неизвестными являются радиусы колец rj и массы mj тел.

3. Силы взаимодействия тел

Рассмотрим силы воздействия всех тел на первое тело на кольце j, масса которого - mj,1. С телом mj,1 связываем траекторную систему координат (n,), где n - нормаль к траектории, а - касательная к ней (см. рис. 1). Сила Fj,1,i,l гравитационного воздействия тела mi,l, находящегося на кольце i, на тело mj,1 равна: Fj,1,i,l = Gmj,1mi,l/rj,1,i,l2, где G - гравитационная постоянная; rj,1,i,l - расстояние тела mi,l от тела mj,1. Тогда проекции силы Fj,1,i,l на оси и запишутся так:

; (3)

, (4)

где nj,1,i,l и j,1,i,l - проекции расстояния rj,1,i,l на оси координат n и .

Отметим, что здесь проекции сил (3) и (4) направлены по векторам осей и .

В треугольнике Omi,lmj,1 (рис. 1) угол между радиусами тел ri и rj будет

j,1,i,l = i,l - j,1, (5)

а расстояние между ними согласно теореме косинусов запишется так:

rj,1,i,l2 = rj 2 + ri 2 - 2rirjcos j,1,i,l. (6)

Тогда проекции этого расстояния на оси n и будут

nj,1,i,l = -(ri cos j,1,i,l - rj); j,1,i,l = ri sin j,1,i,l. (7)

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

,

где согласно (7) при ri =0 для центрального тела nj = rj.

После подстановки (6) и (7) в выражения (3) и (4) и после суммирования по всем телам системы получаем выражения для проекций сил воздействия на тело mj,1 всех остальных тел:

; (8)

, (9)

где mj - масса тел mj,1, а mi - масса тел mi,l.

Чтобы исключить из рассмотрения силу воздействия тела mj,1 на себя, в выражениях (8) и (9) воздействие остальных тел j-того кольца извлечено из общего выражения и записано последним слагаемым. Оно легко получается при замене i на j в предыдущем слагаемом. В пределах суммирования исключение j-того кольца обозначено как i j.

Будем рассматривать такие конфигурации вращающихся структур, для которых выражения для силы (8) и (9) будут давать одну и ту же величину для каждого тела j-го кольца. Это возможно только в том случае, если при прохождении оси n через любое тело кольца j геометрические положения воздействующих тел относительно него не изменится. Последнее условие будет выполняться при условии, что начальный угол тел на кольцах будет принимать значение j,1 = 0 либо j,1 = 0.50. На рис.1 представлен вид структуры, где значения начального угла j,1 последовательно чередуются на соседних кольцах. Вышеуказанному условию удовлетворяют также структуры с произвольным порядком чередования начального угла j,1.

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

Для рассмотренных конфигураций нормаль n является осью симметрии (см. рис. 1). Поэтому углы отклонения j,1,i,l воздействующих тел от оси n имеют, согласно (5), попарно одинаковые по величине и обратные по знаку значения. Следовательно, в выражениях (9) синусы в числителях также попарно одинаковы по величине и обратные по знаку. Так как косинусы этих углов в знаменателях одинаковы, то касательные силы равны нулю. При четном количестве тел N3 еще одно тело будет находиться на оси n симметрично относительно центра O. Так как угол j,1,i,l этого тела равен р, то сила его воздействия в (9) также равна нулю. Итак, проекции всех сил на касательную ось равны нулю, т.е. . Поэтому сила воздействия всех тел осесимметричной многослойной структуры на любое тело на кольце с номером j направлена по нормали n к траектории, т.е. к центру O, и определяется выражением (8).

Для кольца j разность углов, согласно (1) будет

j,1,j,l = j,l - j,1 = 2(l - 1)/N3. (10)

Тогда выражение в знаменателе последнего слагаемого формулы (8) запишется

2[1 - cos(2(l - 1)/N3)] = 4sin2((l - 1)/N3). (11)

После подстановки (11) в (8) направленная к центру O сила воздействия всех тел на любое тело на кольце j будет

, (12)

где ; (13)

ri,,j = ri/rj - отношение радиусов колец i и j.

4. Уравнения движения вращающейся структуры

При воздействии с силой (12) на тело mj,1 с массой mj (рис. 1) оно будет совершать ускоренное движение. В траекторной системе координат (n, ) сила воздействия (12) имеется только вдоль одной оси n, по которой направлено нормальное ускорение wn = v2/, где v - тангенциальная скорость движения тела mj,1, а - радиус кривизны его траектории. Поэтому дифференциальное уравнение его движения запишется так

. (14)

Мы рассматриваем вращающуюся структуру с угловой скоростью и с неизменными радиусами траектории. Поэтому для тела mj,1 радиус кривизны траектории = rj, а скорость v = rj. После подстановки этих величин и силы (12) в уравнение (14) дифференциальное уравнение движения тела mj,1 получаем в виде

, (15)

где j = 1, 2, N2.

Итак, движение тел вращающейся структуры описывается N2 уравнениями (15). Это алгебраическая система уравнений. Как отмечалось ранее, неизвестными являются радиусы колец rj и массы тел mj. При задании радиусов колец rj уравнениями (15) определяются массы тел mj. При необходимости можно задать массы тел mj, а из (15) определить радиусы колец rj.

С целью создания одного алгоритма решения задачи для многослойной структуры с центральным телом и без него, введем исходную массу mi центрального тела и всех тел первого кольца. Обозначим долю в ней массы центрального тела коэффициентом pm0. Тогда масса центрального тела определится как m0 = mi•pm0. Структура без центрального тела задается коэффициентом pm0 = 0.

Анализ возможных применений результатов этой задачи показывает, что целесообразно задавать геометрию осесимметричной многослойной структуры, в том числе и радиусы колец rj, а затем в результате решения уравнений (15) определять массы mj. Поэтому перепишем уравнения (15) в другом виде. Выделим из второго слагаемого безразмерное ускорение тела mj,1, обусловленное воздействием единицы массы тела mi,1 на i-том кольце

, (16)

где ri,,j = ri/rj - является безразмерным отношением радиусов.

Тогда безразмерное ускорение тела mj,1, вызванное всеми телами на i-том кольце, будет

. (17)

Так как тела кольца каждому своему телу сообщают безразмерное ускорение fn3, то обозначим

Aj,j = fn3. (18)

Массы всех тел отнесем к исходной массе mi, и безразмерные массы тел обозначим:

moj = mj/mi; mo0 = m0/mi. (19)

Тогда с учетом (16) - (19) уравнения (15) запишутся в виде следующей системы линейных алгебраических уравнений:

; (20)

где

Bj = Ciroj3 - mo0; roj = rj / r1 ; Ci = r132/(miG). (21)

В системе линейных алгебраических уравнений (20) задаются параметры: N2, N3, j,1, mi, pm0, , rj, а неизвестными являются безразмерные массы moi периферийных тел. Они отнесены к первоначальной массе mi. Масса центрального тела остается неизменной, т.е. m0 = const. А безразмерные массы всех периферийных тел moj определяются в результате решения уравнений (20). Абсолютные массы затем рассчитываются как

mj = mojmi. (22)

5. Решение уравнений

5.1 Описание программы

Система уравнений (20) решается методом Гаусса. Для ее решения разработана программа RtCrcSt2.for на языке FORTRAN. Параметры структуры задаются в файле данных RtCrcSt2.dat. Эти файлы приведены в Приложении.

Программа RtCrcSt2.for предназначена для расчета периферийных масс mj осесимметричной многослойной вращающейся структуры. Кроме того, она позволяет подобрать необходимые параметры, определяющие конфигурацию структуры. И дополнительно, программа RtCrcSt2.for подготавливает файл исходных и начальных данных для программы Galactica, которая имеется в свободном доступе [10]- [11]. Система Galactica позволяет рассчитать динамику многослойной вращающейся структуры и исследовать ее эволюцию.

Как видно из Приложения, в файле данных RtCrcSt2.dat задаются параметры многослойной структуры: N2, N3, smi = mi, pm0, Prd. Для удобства вместо угловой скорости задается период Prd в сидерических годах Yrs = 365.25636042 дня. Как видно из программы RtCrcSt2.for угловая скорость рассчитывается так

= 2/(PrdYrs243600). (23)

В файле данных RtCrcSt2.dat предусмотрено два варианта задания начального угла j,1 первого тела на j-том кольце: в 1-м варианте j,1 = 0 и во 2-м варианте задается последовательное чередование значений 0 и 0.50. Кроме того, возможно задать произвольное чередование этих значений в файле RtStrfi0.dat. При задании значения ключа Kl3 = 1 значения j,1 будут считываться с этого файла.

В файле данных RtCrcSt2.dat радиусы колец задаются с помощью параметра okr. Его значение не может быть меньше 0.5. Радиус j-того кольца рассчитывается так

rj = jokrr1, (24)

где r1 - радиус первого кольца.

Радиус первого кольца определяется из условия его существования при заданной массе центрального тела m0 и первоначальной массе m1 тела первого кольца [1] - [2]. Это условие существования следует из уравнений движения (15) при отбрасывании остальных колец

. (25)

В файле RtCrcSt2.dat при задании ключа Kl1 = 1 радиусы колец будут считываться из файла FIMA2, имя которого можно задать. В этом файле могут быть заданы любые радиусы колец.

Итак, в файле данных RtCrcSt2.dat могут быть заданы параметры некоторых алгоритмов, по которым будет образовываться геометрия многослойной вращающейся структуры, а именно, начальные углы j,1 и радиусы rj колец. Кроме того, при задании ключей Kl1 = 1 и Kl3 = 1 любая геометрия структуры может быть задана в дополнительных файлах.

При большом количестве колец и тел точность решения уравнений (20) методом Гаусса может давать отличные от нуля невязки Vj, которые определяются из уравнений (20) так:

, (26)

где moi - полученные решения в результате работы программы RtCrcSt2.for. Невязки используются для дальнейшего уточнения масс. Число шагов уточнения задается параметром Ltk, а уточнение происходит до тех пор, пока среднее значение модулей невязок не достигнет величины EPS.

Следует отметить, что во всех рассмотренных далее примерах однократное решение уравнений (20) методом Гаусса было достаточным и уточнения не потребовалось.

Последние пять параметров в файле данных RtCrcSt2.dat: eps0, Rob, dt, FIMA5, Idta необходимы для создания файла исходных данных и начальных условий для решения дифференциальных уравнений движения тел многослойной структуры с помощью программы Galactica. Углом eps0 в радианах задается наклон плоскости структуры к некоторой плоскости, которая названа плоскостью экватора. С помощью этого угла двумерная структура пересчитывается в трехмерную. Поэтому в программе Galactica динамика многослойной структуры решается по трем координатам. С помощью плотности Rob определяются радиусы тел, которые затем используются для расчета взаимодействий при столкновении тел. Параметром dt в сидерических столетиях задается номинальный шаг интегрирования в программе Galactica. Параметрами FIMA5 и Idta задается имя файла начальных условий и время его создания. Структура этого файла приведена в работе [10], а также - в Описании системы Galactica [11].

5.2 Работа программы

После запуска исполняемого модуля программы RtCrcSt2.exe и завершения его работы создается 5 файлов: RtStrfi1.dat, ErRtCrSt.dat, 1i, v2103b29.dat (как пример) и RtNcRMEr.dat. В файле RtStfi1.dat записаны номера колец j и углы первых тел j,1. Этот файл может быть взят за основу файла RtStfi0.dat, в котором может быть задана любая конфигурация углов первых тел. Файл ErRtCrSt.dat является файлом ошибок и сообщений. Эти сообщения выдаются одновременно на дисплей и в этот файл. Сообщения об ошибках работы процессора выдаются только на дисплей. Самой распространенным сообщением является сообщение об отрицательной массе тела с номером k. Оно формируется частью программы подготовки файла начальных условий для программы Galactica. Файл начальных условий не создается при отрицательной массе одного из тел. Номер тела k в программе Galactica рассчитывается по алгоритму:

, (27)

где J - номер кольца, на котором находится тело с номером L.

При отсутствии ошибок в файл ErRtCrSt.dat записывается масса центрального тела и начальная масса тела на первом кольце. Далее приводятся номера колец j и их радиусы rj.

В файле типа 1i приводится шаг уточненного решения системы уравнений (20) методом Гаусса и средний модуль невязки при этом решении. Как уже отмечалось, во всех рассмотренных случаях был только один шаг решения. При этом модуль невязки равнялся нулю.

В файле V2103b29.dat приводятся исходные данные и начальные условия для исследования динамики образованной структуры с помощью системы Galactica. Детальное описание этого файла и системы Galactica даны в отмеченных ранее источниках [10]-[11]. В названии этого файла отмечены свойства его структуры: второй вариант начального угла первого тела кольца, 103 кольца и 29 тел на каждом кольце.

В файле RtNcRMEr.dat приведены рассчитанные параметры осесимметричной вращающейся структуры в виде номера кольца j, его радиуса rj, массы одного тела mj на кольце и расчетной невязки массы dmj. В случае однократного решения системы уравнений (20) методом Гаусса невязки масс dmj не рассчитываются и приравниваются нулю.

Файл RtNcRMEr.dat может быть принят за основу при необходимости задания произвольных радиусов rj колец. Во всех файлах программы RtCrcSt2.for радиусы колец приведены в метрах, а массы тел - в кг. В файле, например, V2103b29.dat для программы Galactica все параметры приведены в безразмерных единицах, которые отнесены к масштабу длины Am и массе всей системы mSS. Безразмерное время получено умножением времени на коэффициенты времени kt, а безразмерная скорость - на коэффициент скорости kv. При этом безразмерное время выражено в сидерических столетиях: 100Yrs. В работах [10] - [11] показано, что полученные решения для одного масштаба могут быть с помощью параметров Am, mSS и kt пересчитаны на другой масштаб. Это позволяет однажды полученные решения использовать для задач различных масштабов.

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

Приведенный в Приложении образец программы RtCrcSt2.for допускает N2 = N3 = 1000, т.е. 106 тел. При необходимости количество тел может быть увеличено заданием большего значения l2 в операторе PARAMETER.

6. Примеры решенных задач

В процессе решения этой проблемы были рассчитаны многослойные вращающиеся структуры с числом колец 1, 2, 3, 4, 5, 15, 30, 100, 103, 1000. На кольцах задавалось число тел: 2, 5, 8, 10, 29, 30, 999. Общее число тел достигало одного миллиона. Рассматривались конфигурации с разными вариантами угла j,1 первого тела на каждом кольце, а также с разным чередованием его на соседних кольцах. Были рассчитаны структуры, которые формируются при разных исходных массах. В качестве исходных масс задавались массы Солнца и Земли. Рассматривались структуры с центральным телом и без него.

На рис.1 представлена структура с N2 = 5 кольцами по N3 = 8 тел на каждом кольце и периодом вращения Prd =1 год. Номера колец отсчитываются от центра О, номера тел от оси xО. Углы первых тел на кольцах j,1 последовательно принимают значения 0; 0/2; 0; 0/2; 0. Первоначальная масса mi = 1.98912·1030 кг задана немного больше массы Солнца. Радиусы колец и массы одного из тел на них, отнесенные, соответственно, к радиусу первого кольца и массы тела на нем, равны = 0; 1; 2.005; 2.985; 4.082; 4.980 и = 1.256; 1; 2.957; 2.250; 7.712; 2.973. Здесь кольцом с радиусом = 0 обозначено центральное тело. При этом радиус первого кольца и масса тела на нем соответственно, равны r1 = 1.481·1011 м и m1 = 1.568·1030 кг. Масса всей этой структуры составляет mSS = 2.138·1032 кг. То есть она в 108 раз превышает массу Солнца.

Далее рассмотрим три структуры с N2 = 15 колец и N3 = 30 тел на каждом кольце, которые отличаются чередованием угла j,1 первых тел на кольцах. Созданный программой RtCrcSt2.exe файл исходных данных и начальных условий использовался для интегрирования дифференциальных уравнений движения тел структуры с помощью программы Galactica. В ней имеется опция графического изображения взаимодействующих тел. На рис. 2 показана выдача на экране монитора результатов задачи при углах первых тел на всех кольцах j,1 = 0 после первого шага счета с T = 10-10. Линиями у периферийных тел представлены их вектора скорости. Следует отметить, что показанная на рис. 2 выдача несколько модифицирована. Модификация заключается в изменении цвета изображения и направления вращения.

В верхней части рис. 2 приведена числовая выдача программы Galactica: время счета в сидерических столетиях T = 110-10; наибольшая безразмерная масса momax = 0.00323…; модуль наибольшей безразмерной скорости vmax = 10.927…; исполненный шаг счета Tp = 110-10; количества движения всей системы: Px, Py, Pz; моменты количества движения: Mx, My, Mz. О последующих 7 числах буде сказано ниже. Последними двумя числами представлены величина следующего шага счета T = 10-7 и относительное изменение z-составляющей суммарного момента количества движения Mz = 0, где Mz = (Mz - Mz1)/Mz1, а Mz1 - момент количества движения системы на первом шаге счета. Более детальное описание этих и остальных параметров задачи имеется в [10] и [11]. Эта числовая информация представляет значительный интерес при предварительном анализе динамики структуры, а также точности интегрирования дифференциальных уравнений движения ее тел.



Рис. 2. Вид на экране компьютера осесимметричной многослойной вращающейся структуры при интегрировании дифференциальных уравнений движения программой Galactica за время 0.001·ДT, где ДT = 1·10-7 - шаг интегрирования: N2 = 15; N3 = 30; углы j,1 = 0; Prd =1 год; масса центрального тела равняется массе Солнца. Дополнительно на рис. 2 вертикальными столбцами приведены радиусы колец и массы одного тела на них, отнесенные, соответственно, к радиусу первого кольца r1 = 1.4938371011 м и массе m1 = 8.6849661029 кг тела на нем. Общая масса всей структуры mSS = 3.2184891033 кг.

Как видно из рис. 2, эта структура имеет радиально-лучевое строение. В гомографической динамике [3] ее назвали бы конфигурацией из 15-и вложенных друг в друга концентрических равносторонних 30-угольников. Из второго столбца на рис. 2 видно, что наименьшие массы имеют тела первого кольца. Если масса m0 центрального тела равна массе Солнца, то масса m1 в 2.27 раза меньше m0. С первого кольца по десятое массы тел почти монотонно увеличиваются, затем они уменьшаются. В этой структуре наибольшая масса тела m10 = 12.3·m1.

В структуре на рис. 2 и в двух последующих радиусы колец, согласно (24), увеличиваются прибавлением радиуса первого кольца r1. В структуре, представленной на рис. 1, радиусы колец изменяются по другому закону. Это обусловлено тем, что в этом случае при пропорциональном увеличении радиусов колец решение уравнений (20) приводит к отрицательным массам.

Рис. 3. Осесимметричная многослойная структура при N2 = 15; N3 = 30; углы j,1 = 0 и j,1 = 0.50 чередуются; Prd =1 год; r1 = 1.4938371011 м; m1 = 9.1720581029 кг и mSS = 3.9773021033 кг. Остальные обозначения см. рис. 2.

На рис. 3 представлена многослойная вращающаяся структура с последовательным чередованием углов j,1 первых тел на соседних кольцах. В гомографической динамике [3] ее назвали бы конфигурацией из 15 вложенных друг в друга равносторонних 30-и угольников, повернутых друг относительно друга на угол р/15. По сравнению со структурой, показанной на рис. 2, общая масса этой структуры возросла в 1.23 раза. Масса тел на кольцах в этом случае также растет, начиная от центра. Однако после 10 кольца монотонность изменения масс тел нарушается, и тела с наибольшей массой находятся на последнем кольце.

В структуре, представленной на рис. 4, угол j,1 = 0/2 повторяется через три слоя. В этом случае общая масса системы также увеличилась по сравнению с радиально-лучевым строением (см. рис. 2), однако она меньше чем в предыдущей структуре на рис. 3. Масса тел на кольцах также увеличивается при удалении от центра. При этом на кольцах с углом j,1 = 0/2 массы тел преимущественно больше, чем в соседних кольцах. И наибольшая масса тел находится на таком кольце с номером 12. В случаях трех рассмотренных структур это наибольшая масса тел. Она в 20 раз превышает массу тел на первом кольце.

Большинство гомографических задач представлено вложенными друг в друга правильными многоугольниками [3]. Однако имеются и неправильные многоугольники, например, вложенные друг в друга ромбы [3], [12]. Такая многослойная ромбическая структура (см рис. 20 [3]) может быть создана в виде ряда колец.



Рис. 4. Осесимметричная многослойная структура при N2 = 15; N3 = 30; углы j,1 = 0 и j,1 = 0.50 чередуются через три слоя; Prd =1 год; r1 = 1.4938371011 м, m1 = 8.8890511029 кг и mSS =3.5701641033 кг. Остальные обозначения см. рис. 2.

На рис. 5 она показана в виде N2 = 4 колец с N3 = 2 телами на каждом кольце. Радиусы колец и массы одного тела на них, отнесенные к радиусу первого кольца и массе тела на нем, равны, соответственно, = 0; 1; 1.697; 1.703; 2.40 и = 1161; 1; 193.47; 490.54; 23.76. В отличие от предыдущих структур на первом кольце (см. рис. 5) угол первого тела 1,1 = 0.50. Поэтому оно находится на оси ординат. Эта структура из двух вложенных ромбов идентична структуре в работе [3] на рис. 20.



Рис. 5. Осесимметричная четырехслойная структура при N2 = 4; N3 = 2; углы j,1 = 0.50 и j,1 = 0 чередуются; Prd =2 года; r1 = 1.4916511011 м, m1 = 1.6958541027 кг и mSS = 4.3731511030 кг. Остальные обозначения см. рис. 2.

В четырехслойной кольцевой структуре, представленной на рис. 5 относительные радиусы второго и третьего колец, как показано выше, = 1.697 и = 1.703, т. е. почти равны. Однако углы первых тел на них сдвинуты на угол 0.5•0. Отсюда следует, что можно создавать осесимметричные вращающиеся кольцевые структуры, на отдельных кольцах которых может находиться 2N3 тел. Для этого необходимо, чтобы радиусы двух соседних колец были одинаковы, а углы первых тел на них - разные. В качестве примера на рис. 6 представлена трехслойная структура, у которой радиусы второго и третьего слоя равны. Относительные радиусы и массы в этом случае, соответственно = 0; 1; 1.375; 1.375 и = 21844; 1; 175.8; 193.8. Отсюда видно, что при равных радиусах второго и третьего колец массы тел на них отличаются. Эта структура идентична гомографической конфигурации состоящей из правильного шестиугольника, который содержит внутри себя концентрический равносторонний треугольник (см. рис. 26 [3]).

В примере, представленном на рис. 6, наблюдается некоторая особенность. Математически эта структура задается тремя кольцами, т. е. N2 = 3, а физически существуют два кольца или два слоя. Однако во втором слое тела имеют два вида масс. Итак, представленные примеры показывают, что широкий набор конфигураций гомографической динамики может быть выражен осесимметричными многослойными вращающимися структурами.

Рис. 6. Осесимметричная трехслойная структура при N2 = 3; N3 = 3; углы j,1 = 0; 0.50 и 0; Prd = 1.5 года; r1 = 1.4919901011 м, m1 = 9.0185541026 кг и mSS = 2.9718831030 кг. Остальные обозначения см. рис. 2.

7. Дополнительные уточнения

Как уже отмечалось, были получены структуры с числом колец N2 = 103 и числом тел N3 = 29; с N2 = 1000 и N3 = 999. В первой структуре общее число тел было N = 2988, а во второй - N = 999001. Кроме того, рассчитывались структуры с исходной массой mi равной массе Земли. Число колец в этих структурах было N2 = 6; 15 и 103. Созданные файлы исходных данных и начальных условий для программы Galactica совпали с файлами при исходной массе mi немного большей массе Солнца. Это означает, что результаты этой программы, представленные на рис. 2 - 6 с массой Солнца, являются такими же в случае структур с массой Земли.

Как уже отмечалось, в ряде случаев при решении систем линейных уравнений (20) отдельные массы тел могут быть отрицательными. Поэтому необходимо изменением параметров структуры добиваться того, чтобы все массы mj были положительными. Например, в структуре N2 = 4 и N3 = 2 при втором варианте углов j,1 отрицательная масса была у тел последнего слоя. При увеличения периода вращения Prd в два раза масса у всех тел стала положительной. Во втором примере с N2 = 5 и N3 = 8 отрицательные массы были у тел 3-его и 5-ого колец. При увеличении радиуса 3-его кольца массы на нем стали положительными. Затем при увеличении периода Prd в 1.4 раза все массы стали положительными.

В структурах с большим числом колец для устранения отрицательных масс требуется приобрести определенный опыт при работе с программой RtCrcSt2.for. Следует отметить, что эта программа является вторым вариантом решения рассмотренной задачи. При первоначальном ее решении мы свели уравнения (15) к системе нелинейных алгебраических уравнений. И в первом варианте программа RtCrcStr.for задача решалась методом итераций. Первоначально массы задавались по одному из 4-х алгоритмов. Затем после каждой итерации они уточнялись. Время счета одной итерации существенно меньше времени решения программой RtCrcSt2.for. Однако количество итераций для задач с большим числом колец требуются десятки тысяч. Поэтому время решения всей задачи в первом варианте значительно больше. В связи с этим задачи решались на суперкомпьютерах Сибирского суперкомпьютерного центра СО РАН в г. Новосибирске.

При первом методе решения задачи программой RtCrcStr.for после последней итерации производилось уточнение радиусов колец. Это позволяло использовать следующий прием для устранения отрицательных масс. Задача запускалась несколько раз на счет с малым количеством итераций Ltk. С каждым запуском программы происходило уточнение радиусов колец. Это позволяло в дальнейшем устранить отрицательные массы. Затем счет с большим количеством итераций приводил к окончательному решению. При этом массы определялись с относительной средней погрешностью EPS, которая задавалась в файле данных RtCrcStr.dat.

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

В вышеупомянутых работах по гомографической динамике, особенно в работе [3], вопросы устойчивости вращающихся структур рассматриваются аналитическими методами в рамках гамильтоновой динамики. Как отмечено в [3], существует более ста определений понятия устойчивости, а в космической динамике можно встретить до тридцати таких определений. К сожалению, эти методы не могут дать прямого ответа на то, как поведет себя конкретная структура. Поэтому необходимо выполнять численное интегрирование дифференциальных уравнений движения тел структуры, чтобы установить динамику структуры и ее эволюцию. При этом точность интегрирования должна быть такой, чтобы погрешность решения задачи не оказывала влияние на поведение структуры. Этому требованию удовлетворяет программа Galactica.

Рис. 7. Начало разрушения осесимметричной структуры при N2 = 5; N3 = 8; углы j,1 = 0 и j,1 = 0.50 чередуются; Prd =1 год; структура совершила 1.6 оборота, а в начальном состоянии она показана на рис. 2.

С позиций силового взаимодействия все многослойные вращающиеся структуры являются неустойчивыми. Однако время существование структуры может изменяться в широких пределах. На рис. 7 представлена уже рассмотренная на рис. 1 структура после 1.6 оборота. Изменения в структуре начались из-за разрушения первого, т. е. внутреннего, кольца. Его восемь тел попарно объединились. Дальнейшее их движение приводит к разрушению всей структуры.

С помощью программы Galactica были выполнены исследования динамики осесимметричных структур в работах [9], [13] - [16], а также в находящейся в печати работе: Смульский И. И. «Осесимметричное кулоновское взаимодействие и неустойчивость орбит». Эти исследования показали, что существует широкий спектр поведения осесимметричных структур. Наибольший интерес представляет время существование структуры. Это время тем больше, чем меньше безразмерная скорость вращения системы. При определенной скорости на время существования влияет точность расчета параметров структуры. Например, в первом варианте программы RtCrcStr.for массы определялись с заданной невязкой EPS. С уменьшением величины EPS время существования структуры может увеличиваться. Если осесимметричная структура входит в состав других тел, то время ее существования также может существенно увеличиваться. В этом случае тела структуры совершают колебания вдоль орбиты. Это нарушает их тенденцию к сближению, что может приводить к неограниченному времени существование структуры [13].

8. Заключение

Выполненные исследования показывают, что разработанные алгоритм и программа позволяет создавать осесимметричные многослойные вращающиеся структуры с любым количеством тел. Представляет интерес изучения их динамических свойств с помощью программы Galactica. Эти исследования могут привести к использованию рассматриваемых структур в различных задачах небесной и космической динамики. Например, однослойные осесимметричные структуры были использованы для составных моделей вращения Земли [13] и Солнца [15] - [16]. При этом были получены качественные представления об эволюции Земной оси. Кроме того, колебания периферийных тел в этой модели навели на мысль о возможном колебании континентов в широтном направлении [13]. Составная модель вращения Солнца позволила получить избыток вращения перигелия Меркурия, которого недоставало при описании динамики Солнечной системы с помощью закона тяготения Ньютона.

Представляет интерес использования рассматриваемых структур для моделирования дисков планет и дисковых галактик. После создания модели этих объектов их динамика может быть исследовано с помощью программы Galactica.

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

Литература

1. Смульский И.И. Теория взаимодействия. - Новосибирск: Из-во Новосиб. ун-та, НИЦ ОИГГМ СО РАН, 1999 г. - 294 с. http://www.ikz.ru/~smulski/TVfulA5_2.pdf.

2. Смульский И.И. Осесимметричная задача гравитационного взаимодействия N-тел// Математическое моделирование, 2003, т. 15, № 5, с. 27-36. http://www.smul1.newmail.ru/Russian1/IntSunSyst/Osvnb4.doc.

3. Гребеников Е.А. Математические проблемы гомографической динамики. М.: МАКС Пресс, 2010. - 256 с.

4. Диарова Д.М., Земцова Н.И., Ихсанов Е.В. Существование центральных конфигураций в одной модели для ньютоновой проблемы восьми тел. / Труды ИСА РАН «Динамика линейных и нелинейных систем», т. 25(1). - 2006, с.6471.

5. Гуцу В.Д., Диарова Д.М., Земцова Н.И. Исследование устойчивости стационарных решений ромбоподобной ограниченной задачи десяти тел. / В сб. Теоретические и прикладные задачи нелинейного анализа. М.: ВЦ РАН. - 2007. С. 99-109.

6. Гребеников Е.А., Диарова Д.М., Земцова Н.И. Существование устойчивых ромбоподобных центральных конфигураций в смысле Уинтнера для ньютоновой модели девяти тел. / В сб. Теоретические и прикладные задачи нелинейного анализа. М.: ВЦ РАН. - 2006. С.65-76.

7. Силушик А. Konieczne i wystarczajace warunki istnienia homograficznych rozwianzan w specjalnym zagadnieniu 7-u i 10-u cial. / Proceeding of the international workshop. Brest: Изд-во БрГУ. - 2003. С. 206-210.

8. Гребеников Е.А., КозакСковородкина Д., Якубяк М. Методы компьютерной алгебры в проблеме многих тел. - М.: Изво РУДН, 2002. - 212с.

9. Смульский И. И. Многослойные кольцевые структуры// Письма в ЭЧАЯ, 2011, т. 8, No. 5(168), с. 737-743. http://www.ikz.ru/~smulski/Papers/MnsKoStr4c.pdf.

10. Smulsky J.J. The System of Free Access Galactica to Compute Interactions of N-Bodies. I.J.Modern Education and Computer Science. - 2012, 11, pp. 1-20. http://www.mecs-press.org/ doi:10.5815/ijmecs.2012.11.01.

11. Смульский И.И. Система Galactica. http://www.ikz.ru/~smulski/GalactcW/.

12. Diarova D., Zemtsova N.I. The Instability of the Rhombus-Like Central Сonfigurations in Newton 9-body Problem / Computer Algebra in Scientific Computing (CASC 2006), LNCS 4194, Springer-Verlag, Berlin Heidelberg, Chisinau, September 11-15. - 2006, рp. 141-148.

13. Мельников В. П., Смульский И.И., Смульский Я.И. Составная модель вращения Земли и возможный механизм взаимодействия континентов // Геология и Геофизика, 2008, №11, с. 1129-1138. http://www.ikz.ru/~smulski/Papers/RGGRu190.pdf.

14. Smulsky J. J. Galactica Software for Solving Gravitational Interaction Problems // Applied Physics Research, 2012, Vol. 4, No. 2, pp. 110-123. http://dx.doi.org/10.5539/apr.v4n2p110.

15. Смульский И.И. Численное моделирование эволюции спутника вращающегося тела / В сб. Теоретические и прикладные задачи нелинейного анализа. Российская Академия Наук: ВЦ им. А.А. Дородницына. М.: ВЦ РАН А.А. Дородницына. - 2008. С. 100-117. http://www.ikz.ru/~smulski/Papers/ModSun07c.pdf.

16. Smulsky J.J. New Components of the Mercury's Perihelion Precession // Natural Science, 2011, Vol. 3, No.4, pp. 268-274. doi:10.4236/ns.2011.34034. http://www.scirp.org/journal/ns.

Приложение

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

1. Содержание файла исходных данных RtCrcSt2.dat

No	Описание параметра	Пара-метр	Величина, начиная с 60 позиции
1	Количество колец 10.06.2013	N2	103
2	Количество тел в кольце	N3	29
3	Коэффициент начальных радиусов колец (>0.5)	okr	1.00000D+0
4	Ключ ввода радиусов из файла FIMA2	Kl1	0
5	Имя файла с исходными радиусами	FIMA2	RtNcJR10.dat
6	Первоначальная масса центр.тела и 1-го кольца	smi	1.98912D+30
7	Часть массы центрального тела	pm0	0.99000D+0
8	Вариант угла первых тел в кольцах	Ivr	2
9	Углы первых тел на кольцах из RtStrfi0.dat.	Kl3	0
10	Период обращения в сидерических годах	Prd	1.00000D+0
11	Относительная средняя погрешность	EPS	1.00000D-8
12	Количество шагов итераций	Ltk	10
13	Угол наклона колец к плоскости экватора	eps0	0.40900D+0
14	Абсолютная плотность тел, kg/m^3	Rob	5000.
15	Шаг интегрирования по времени	dt	1.00000D-7
16	Имя файла начальных условий	FIMA5	v2103b29.dat
17	Дата создания файла начальных условий	Idta	20130723

2. Текст программы RtCrcSt2.for

C************************************************************

C* ***** RtCrcSt2.for ******

C Вращающиеся кольцевые структуры

C Joseph J. Smulsky

C Последнее изменение 18 ч. 04.06.2014 г. !!-1str

C************************************************************

C---------------------------------------------------

C Нужно задать значения параметров!

PARAMETER (l2=1000,l4=6000) !

C ---------------------------------------------------

IMPLICIT DOUBLE PRECISION (A-H,O-Z)

DIMENSION A1(l2,l2+1),A2(l2,l2+1),fi(l2,l2),cm(l2),

*R(l2),Ro(l2),cmo(l2),dmo(l2),dmm(l4),

*omo(l2*l2+1),Xo(l2*l2+1,3),Uo(l2*l2+1,3),X1(l2*l2+1,3),

*U1(l2*l2+1,3),Pns(3),Xc(3),AM1(3),Pn(3),X(l2*l2+1,3),

*U(l2*l2+1,3),Sp(l2*l2+1,3),Et(l2*l2+1),Ra(l2*l2+1),dUp(l2*l2+1,3)

INTEGER*4 N,Nstr,N2,N3,Lt,Ltk,Lt4,Lt5,Idta,q

CHARACTER*12 FIMA1,FIMA2,FIMA3,FIMA4,FIMA5,FIMA6,FIMA7

CHARACTER*3 Lt5Ima

CHARACTER*4 ITER

FIMA1 ='RtCrcSt2.dat'

FIMA3 ='ErRtCrSt.dat'

FIMA4 ='RtNcRMEr.dat' ! Номер кольца, радиус, масса тела, ошибка

FIMA6 ='RtStrfi0.dat'

FIMA7 ='RtStrfi1.dat'

C--- Считывание файла RtCrcSt2.dat.

OPEN(UNIT=1,FILE=FIMA1,STATUS='OLD')

READ(1,39,ERR=40)N2

READ(1,39,ERR=40)N3

READ(1,54,ERR=40)okr

READ(1,38,ERR=40)Kl1

READ(1,'(59X,A12)',ERR=40)FIMA2

READ(1,54,ERR=40)smi

READ(1,54,ERR=40)pm0

READ(1,38,ERR=40)Ivr

READ(1,38,ERR=40)Kl3

READ(1,54,ERR=40)Prd

READ(1,54,ERR=40)EPS

READ(1,39,ERR=40)Ltk

READ(1,54,ERR=40)eps0

READ(1,54,ERR=40)Rob

READ(1,54,ERR=40)dt

READ(1,'(59X,A12)',ERR=40)FIMA5

READ(1,39,ERR=40)Idta

CLOSE(1)

38 FORMAT(59X,I5)

39 FORMAT(59X,I12)

53 FORMAT(59X,F11.4)

54 FORMAT(59X,D11.5)

GO TO 1

40 WRITE(*,*)'*01Rt Error reading file RtCrcSt2.dat *** ',

*N2,N3,okr,Kl1,FIMA2,smi,pm0,Ivr,Kl3,Prd,EPS,Ltk,eps0,

*Rob,dt,FIMA5,Idta

OPEN(UNIT=3,FILE=FIMA3)

WRITE(3,*)'*01Rt Error reading file RtCrcSt2.dat *** ',

*N2,N3,okr,Kl1,FIMA2,smi,pm0,Ivr,Kl3,Prd,EPS,Ltk,eps0,

*Rob,dt,FIMA5,Idta

CLOSE(3)

GO TO 150

1 pi = 3.1415926535897932385D+00

p2 = 2.0D+00*pi

G = 6.67259D-11

Yrs = 365.25636042D+00

dsc = 24D+00*3600D+00

DZI = 0.0D+00

PND = 1.5D+00

PND2 = 1.0D+00/3.0D+00

dfi0 = p2/N3

Prdsc = Prd*Yrs*dsc

w = p2/Prdsc

w2 = w*w

fn3 = DZI

Vs = DZI

N = N2*N3 + 1

Nstr = l4/2

Lt4 = 1

Lt5 = 1

k1 = 1

C okm = 0.05D+00

DO l = 2,N3

fn3 = fn3 + 0.25D+0/DSIN((l-1)*pi/N3)

END DO ! l

C--- Углы первых тел на кольцах

IF(Ivr.EQ.1) THEN

DO I = 1,N2

fi(I,1) = DZI

END DO ! I

ELSE

I2 = 0

DO I = 1,N2

I2 = I2 +1

fi(I2,1) = DZI

I2 = I2 +1

fi(I2,1) = 0.5D+0*dfi0

IF(I2.EQ.N2) GO TO 3

END DO ! I

3 END IF ! Ivr

C--- Углы первых тел на кольцах из файла RtStrfi0.dat.

IF(Kl3.EQ.1) THEN

OPEN(UNIT=6,FILE=FIMA6)

DO I = 1,N2

READ(6,*)I1,fi(I,1)

END DO ! I

CLOSE(6)

ELSE

OPEN(UNIT=7,FILE=FIMA7)

DO I = 1,N2

WRITE(7,*)I,fi(I,1)

END DO ! I

CLOSE(7)

END IF !

C--- Масса центрального тела и предварительная масса тела первого кольца

cm0 = smi*pm0

cmo0 = pm0

sm = smi - cm0

cm(1) = sm/N3

C--- Радиусы колец

IF(Kl1.EQ.1) GO TO 27

R(1) = (G*(cm0 + cm(1)*fn3)/w2)**PND2

OPEN(UNIT=3,FILE=FIMA3)

WRITE(3,*) 'Initial masses of the central and the body of the',

*' 1-st circle ',cm0,cm(1)

WRITE(3,*) 'Radii of the circles'

WRITE(3,*) ' 1 ',R(1)

DO I = 2,N2

R(I) = okr*I*R(1)

WRITE(3,*) I,R(I)

END DO ! I

CLOSE(3)

GO TO 11

27 OPEN(UNIT=2,FILE=FIMA2)

DO J = 1,N2

READ(2,*)J1,R(J)

END DO ! J

CLOSE(2)

OPEN(UNIT=3,FILE=FIMA3)

WRITE(3,*) 'Initial masses of the central and the body of the',

*' 1-st circle ',cm0,cm(1)

WRITE(3,*) 'Radii of the circles'

WRITE(3,*) ' 1 ',R(1)

DO I = 2,N2

WRITE(3,*) I,R(I)

END DO ! I

CLOSE(3)

11 okt = 1/(100*Prdsc)

wo = w/okt

Am = R(1)

C--- Относительные радиусы тел

DO I = 1,N2

Ro(I) = R(I)/Am

END DO ! I

TMP = Am*Am*w2

TMP = TMP*Am

Ci = TMP/(G*smi)

C--- Углы всех тел на кольцах

DO J = 1,N2

DO l = 2,N3

fi(J,l) = fi(J,l-1) + dfi0

END DO ! l

END DO ! J

C--- Cила Qs воздействия на тело кольца j всех тел всех колец

DO J = 1,N2

DO I = 1,N2

Qs = DZI

RIJ = R(I)/R(J)

DO l = 1,N3

dfi = fi(I,l) - fi(J,1)

C--- Предотвращение деления на нуль в случаях

C--- сил воздействия тела на себя.

IF(l.EQ.1.AND.I.EQ.J) dfi = pi

C--- Расчет сил воздействия на тело данного кольца тела другого кольца

TMP = RIJ*DCOS(dfi)

Q1 = (1-TMP)/(1+RIJ*RIJ-2*TMP)**PND

C--- Обнуление сил воздействия тела кольца на себя.

IF(l.EQ.1.AND.I.EQ.J) Q1 = DZI

Qs = Qs + Q1

END DO ! l

A1(J,I) = Qs

A2(J,I) = Qs

END DO ! I

TMP = Ro(J)*Ro(J)

A1(J,N2+1) = Ci*TMP*Ro(J) - cmo0

A2(J,N2+1) = A1(J,N2+1)

END DO ! J

CALL GAUSS(N2,A1,cmo,S,l2)

C--- Уточнение масс

DO Lt = 1,Ltk

DO J = 1,N2

Qs = DZI

DO I = 1,N2

Qs = Qs + cmo(I)*A2(J,I)

A1(J,I) = A2(J,I)

END DO ! I

A1(J,N2+1) = A2(J,N2+1) - Qs

Vs = Vs + DABS(A1(J,N2+1))

C WRITE(*,*) '3',' Lt = ',Lt,J,cmo(J),A1(J,N2+1)

END DO ! J

Vs = Vs/N2

IF(Vs.LT.EPS) GO TO 5

CALL GAUSS(N2,A1,dmo,S,l2)

dmm0 = DZI

DO J = 1,N2

dmoo = dmo(J)/cmo(J)

dmm0 = dmm0 + DABS(dmoo)

END DO ! J

dmm(Lt4) = dmm0/N2

IF(dmm(Lt4).LT.EPS) GO TO 5

DO J = 1,N2

C Ослабление приращения dmo(J)

cmo(J) = cmo(J) + 0.1D+00*dmo(J)

END DO ! J

C--- Запись разностей итерационного процесса

Lt4 = Lt4 + 1

IF(Lt.EQ.Ltk) GO TO 7

IF(Lt4.LE.Nstr) GO TO 9

7 WRITE(Lt5Ima,'(I3)') Lt5

ITER = Lt5Ima//'i'

OPEN (UNIT=2,FILE=ITER)

DO in = 1,Lt4-1

inLt = (Lt5-1)*Nstr + in

WRITE(2,*)inLt,dmm(in)

END DO ! in

CLOSE(2)

Lt5 = Lt5 + 1

Lt4 = 1

9 END DO ! Lt

GO TO 15

5 WRITE(Lt5Ima,'(I3)') Lt5

ITER = Lt5Ima//'i'

OPEN (UNIT=2,FILE=ITER)

DO in = 1,Lt4

inLt = (Lt5-1)*Nstr + in

WRITE(2,*)inLt,dmm(in)

END DO ! in

CLOSE(2)

15 ssm = cm0

DO J = 1,N2

cm(J) = cmo(J)*smi

ssm = ssm + cm(J)*N3 !!-1str

END DO ! J

C--- Запись номера кольца, радиуса, массы тела

C--- и ошибка в файл RtNcRMEr.dat.

17 OPEN(UNIT=4,FILE=FIMA4)

DO I = 1,N2

WRITE(4,*)I,R(I),cm(I),dmo(I)

END DO ! I

C--- Массы тел в программе Galactica и параметры в плоскости орбиты

23 omo(1) = cm0/ssm

TMP = G*ssm/(okt*okt)

Am = TMP**PND2

Robo = Rob*Am*Am*Am/ssm

k = 1

DO J = 1,N2

Ro(J) = Ro(J)*R(1)/Am

DO l = 1,N3

k = k+1

omo(k) = cm(J)/ssm

Xo(k,1) = Ro(J)*DCOS(fi(J,l))

Xo(k,2) = Ro(J)*DSIN(fi(J,l))

Xo(k,3) = DZI

Uo(k,1) = -Ro(J)*wo*DSIN(fi(J,l))

Uo(k,2) = Ro(J)*wo*DCOS(fi(J,l))

Uo(k,3) = DZI

END DO ! l

END DO ! J

C--- Координаты и скорости периферийных тел в экваториальной плокости

DO k = 2,N

X1(k,1) = Xo(k,1)

X1(k,2) = Xo(k,2)*DCOS(eps0) - Xo(k,3)*DSIN(eps0)

X1(k,3) = Xo(k,2)*DSIN(eps0) + Xo(k,3)*DCOS(eps0)

U1(k,1) = Uo(k,1)

U1(k,2) = Uo(k,2)*DCOS(eps0) - Uo(k,3)*DSIN(eps0)

U1(k,3) = Uo(k,2)*DSIN(eps0) + Uo(k,3)*DCOS(eps0)

END DO ! k

C--- Координаты и скорости центрального тела в экваториальной плокости

DO q = 1,3

X1(1,q) = DZI

U1(1,q) = DZI

END DO ! q

C--- Параметры центра масс.

30 DO q = 1,3

Pns(q) = DZI

Xc(q) = DZI

AM1(q) = DZI

END DO ! q

DO k = k1,N

IF(omo(k).EQ.DZI) GO TO 32

DO q = 1,3

Xc(q) = Xc(q) + omo(k) * X1(k,q)

Pn(q) = omo(k) * U1(k,q)

Pns(q) = Pns(q) + Pn(q)

END DO ! q

32 END DO ! k

C--- Приведение к центру масс.

DO k = k1,N

IF(omo(k).EQ.DZI) GO TO 16

DO q = 1,3

X(k,q) = X1(k,q) - Xc(q)

U(k,q) = U1(k,q) - Pns(q)

END DO ! q

AM1(1) = AM1(1)-omo(k)*U(k,2)*X(k,3)+omo(k)*U(k,3)*X(k,2)

AM1(2) = AM1(2)+omo(k)*U(k,1)*X(k,3)-omo(k)*U(k,3)*X(k,1)

AM1(3) = AM1(3)+omo(k)*U(k,2)*X(k,1)-omo(k)*U(k,1)*X(k,2)

16 END DO ! k

C WRITE(*,*)' 1. AM1(1) = ',Am1(1),Amn,vkn

DO q = 1,3

Pns(q) = DZI

END DO ! q

omm = DZI

DO k = k1,N

DO q = 1,3

Pn(q) = omo(k) * U(k,q)

Pns(q) = Pns(q) + Pn(q)

dUp(k,q) = DZI

Sp(k,q) = DZI

END DO ! q

Et(k) = DZI

TMP = 3*omo(k)/(4*pi*Robo)

IF(TMP.GE.DZI) GO TO 19

OPEN(UNIT=3,FILE=FIMA3)

WRITE(3,*)'*05Rt Computing Error of Programm RtCrcSt2 *** ',

*' after terminal step the mass of body k = ',k,' omo(k) = ',

*omo(k),' is negative.'

WRITE(*,*)'*05Rt Computing Error of Programm RtCrcSt2 *** ',

*' after terminal step the mass of body k = ',k,' omo(k) = ',

*omo(k),' is negative.'

CLOSE(3)

GO TO 150

19 Ra(k) = TMP**PND2

IF(omo(k).LT.omm) GO TO 13

omm = omo(k)

13 END DO ! k

C WRITE(*,*)' 2. Ra(1) = ',Ra(1),Amn,vkn

A = 1.0D+0

B = 0.5D+0

C = 0.25D+0

TMP = (N+1)**PND2

Mu = NINT(TMP)

C WRITE(*,*)TMP,' MU = ',MU

T = DZI

Um = DZI

Spsx = DZI

Spsy = DZI

Spsz = DZI

E = DZI

Em = DZI

Ett = DZI

dtp = DZI

C--- Создание файла начальных условий FIMA5 для программы Galactica.

OPEN(UNIT=5,FILE=FIMA5)

WRITE(5,*)T,omm,Um,dtp,Pns(1),Pns(2),Pns(3),AM1(1),AM1(2),

*AM1(3),Spsx,Spsy,Spsz,E,Em,Ett,dt,1,1,N,A,B,C,Mu

DO k = k1,N

WRITE(5,*) omo(k),(X(k,q),q=1,3),(U(k,q),q=1,3),

*(dUp(k,q),q=1,3),(Sp(k,q),q=1,3),Ra(k),Et(k)

END DO ! k

C--- Information line at end of file

WRITE(5,*)Idta,N,ssm,Am,okt,0,0,0,0,0,0,0,0,0,0

CLOSE(5)

150 CONTINUE

END

C --- Подпрограмма решения систем линейных уравнений методом Гаусса GAUSS

SUBROUTINE GAUSS(N,A,X,S,l2)

DIMENSION A(l2,l2+1),X(l2)

DOUBLE PRECISION A,X,S,S1,R

N1 = N+1

DO 25 K = 1,N

K1 = K+1

S = A(K,K)

J = K

DO 21 I = K1,N

R = A(I,K)

IF(DABS(R).LE.DABS(S)) GOTO 21

S = R

J = I

21 CONTINUE

IF(S.EQ.0) RETURN

IF(J.EQ.K) GOTO 23

DO 22 I = K,N1

R = A(K,I)

A(K,I) = A(J,I)

22 A(J,I) = R

23 DO 24 J = K1,N1

24 A(K,J) = A(K,J)/S

DO 25 I = K1,N

R=A(I,K)

DO 25 J = K1,N1

25 A(I,J) = A(I,J) - A(K,J)*R

X(N) = A(N,N1)

DO 27 I = N-1,1,-1

S1 = A(I,N1)

DO 26 J= I+1,N

26 S1 = S1 - A(I,J)*X(J)

27 X(I) = S1

RETURN

END

Размещено на Allbest.ru

...