В связи с потеплением климата в последние 40 - 50 лет проблема его изменения стала одной из важных тем современной науки. Чтобы установить причины потепления нужно понять механизмы изменения климата за прошедшие эпохи. Они изучаются в Астрономической теории изменения климата. В 2011 г. интерес к этой теории был актуализирован В.А. Большаковым в его совместной с А.П. Капицей статье [7] (см. Литературу в конце статьи) в Вестнике РАН. Однако в этой статье и других работах В.А. Большакова звучали обвинения в адрес Милутина Миланковича, который почти 100 лет назад завершил создание Астрономической теории ледниковых периодов (так М. Миланкович ее назвал). При знакомстве с аргументами В.А. Большакова выяснилось, что существует недопонимание Астрономической теории изменения климата. В 2013 г. я выступил [17] с разъяснением отдельных ее положений, а также привел некоторые свои результаты, которые открывали новые перспективы этой теории.

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

       Эту статью, за исключением рис. 13 и рис. 14, я подготовил для журнала "Вестник РАН". Я считал, что проблема изменения климата, намеченная в статьях [7] и [17], в этой статье будет достаточно полно изложена, чтобы ученые могли в дальнейшем успешно использовать ее и развивать. Однако у редколлегии Вестника РАН я не встретил понимания, и статья была депонирована в ВИНИТИ. В конце статьи, в Приложении, представлена моя переписка с журналом.

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

       Смульский И.И. Основные положения и новые результаты астрономической теории изменения климата / Институт криосферы Земли СО РАН. – Тюмень, 2014. – 30 с.: ил: 16.– Библиогр.: 44 назв. – Рус. Деп . в ВИНИТИ РАН 30.09.2014, № 258-В2014.

------------------------------------------------------------------------------------------------------------------

------------------------------------------------------------------------------------------------------------------

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

       В истории Земли известны многократные чередования потеплений и похолоданий [1]-[2]. Десять тыс. лет назад сошел ледник, который покрывал северную и среднюю полосу Европы. С другой стороны, в современных, практически безжизненных полярных районах, некогда была обильная растительность, росли деревья, и обитал обширный животный мир: мамонты, шерстистые носороги, бизоны, лошади и другие животные. Какая причина таких колебаний климата на Земле?

       Еще в 19 веке Луи Агассис [3], Ж. Адемар [4], Джеймс Кроль [5] и др. вынашивали идеи о том, что изменение параметров орбиты Земли и оси ее вращения может приводить к изменению количества тепла Солнца, которое поступает на поверхность Земли на разных ее широтах. К концу 19 века достижения в небесной механики стали таковы, что можно было рассчитать изменения орбитальных и вращательных характеристик Земли, и в начале 20^ого века Милутин Миланкович [6] завершил создание Астрономической теории ледниковых периодов, которая по существу является астрономической теорией изменения климата. В этой теории на основании трех параметров: эксцентриситета e орбиты Земли, углового положения перигелия φ_pγ и угла наклона ε плоскости орбиты к плоскости экватора Земли рассчитывается инсоляция Земли на разных ее широтах.

       Процесс облучения Земли Солнцем и величину получаемого ее поверхностью тепла принято называть инсоляцией (in-sol, in - приставка к глаголам, означающим приносить, приводить; solis - солнце). В то же время необходимо помнить, что существуют разные стороны этого процесса, которые обозначены другими терминами: иррадиация, освещенность, облучение и т.д. Инсоляция может быть мгновенная, за сутки, за какой-то сезон или полугодие, за год.

       Так как параметры e, φ_pγ и ε изменяются и колеблются с периодами десятки тысяч лет, то от эпохи к эпохе может быть рассчитаны инсоляции, например, за летнее полугодие на широте 65°, и по изменению такой инсоляции можно судить об изменении климата. Однако амплитуда колебаний эксцентриситета e и угла наклона ε в теории М. Миланковича была небольшой. Например, угол ε колебался в диапазоне ± 1°. Эти колебания могли дать колебания температуры также порядка 1÷2°C. Поэтому, как во времена М. Миланковича, так и в последнее время астрономическая теория изменения климата вызывала сомнение ученых [7], занимавшихся изучением палеоклимата.

       После М. Миланковича, с периодичностью нескольких десятков лет, его исследования повторяли разные группы ученых [8]-[13]. Они уточняли эволюцию параметров e, φ_pγ, ε и первоначальный расчет за 600 тыс. прошедших лет [6] продолжали за более длительное время, например, до 30 млн. лет [10]. Однако главный результат: колебания эксцентриситета e и угла наклона ε оставались по-прежнему небольшими.

       Во второй половине 20^ого века выяснилось, что в морских осадках обнаруживаются колебания изотопа кислорода O¹⁸. Были проведены обширные исследования по всему мировому океану. Результаты исследований обобщены в виде стандартных зависимостей относительного содержания изотопа кислорода δO¹⁸ по толщине осадков, которая соотносится со временем T [14]. Эти зависимости названы по начальным буквам названия программы CLIMAP [15] или имен авторов LR-4 [16]. Предполагается, что более легкий изотоп O¹⁶ в составе воды испаряется и накапливается в ледниках. Поэтому было принято, что избыток изотопа кислорода O¹⁸, т.е. превышение δO¹⁸ над средним уровнем, пропорционален объему льда, накопленному в ледовых покровах Земли.

       Один из больших периодов кривой изменения δO¹⁸ равен 100 млн. лет. Он близок к одному из периодов изменения эксцентриситета орбиты Земли e. Поэтому было принято, что эксцентриситет e оказывает существенное влияние на климат Земли. Например, в работе [7] даже утверждается, что М. Миланкович не учитывает прямого влияния эксцентриситета на климат Земли. В действительности теория инсоляции М. Миланковича дает математически строгую зависимость инсоляции от параметров e, φ_pγ и ε [17]. Это с одной стороны. А с другой стороны, предстоит еще большая работа по выяснению обстоятельств и причин образования определенных свойств морских осадков, ледовых кернов современных ледников и других палеоклиматических наборов данных. Такая же работа должна быть выполнена по выяснению и надежности выводов астрономической теории климата. Этой работой мы занимаемся в Институте криосферы Земли на протяжении двух десятков лет.

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

       Вначале было отмечено, что решения М. Миланковича последовательно повторяли несколько поколений исследователей. Однако все они шли одним и тем же путем, который за столетия сложился в небесной механике. Мы пошли другим путем. Мы не копируем уравнения предшественников, а сами их выводим из первоначальных основ. Во-вторых, при выводе стремимся ввести минимальные упрощения. И, в-третьих, задачи решаем численными методами, стремясь использовать самые высокоточные их варианты или создать новые. В настоящей работе мы излагаем Астрономическую теорию изменения климата в результатах, полученных нами. Вначале речь пойдет об инсоляции Земли, затем об эволюции орбитального движения и об эволюции оси Земли. По первым двум задачам наши независимые исследования подтверждают исследования предшественников. А результаты вращательного движения - другие. Амплитуда колебаний угла наклона ε в семь раз больше по сравнению с прежними теориями. Эти колебания дают такие колебания инсоляции, которые могут объяснить имевшие место колебания климата в прошлом. Эта существенная разница по результатам задачи о вращении Земли требует всесторонней проверки. В заключительной части будет уделено внимание проверкам результатов задачи о вращении Земли.

       В центре небесной сферы 1 (см. рис. 1) находится наблюдатель M. Его горизонт пересекает небесную сферу по кругу HH'. Перпендикуляр к плоскости горизонта пересекает небесную сферу в точке зенита Z. Ось вращения Земли, обозначенная вектором угловой скорости Земли , пересекает небесную сферу в точке северного полюса N. Угол φ между и плоскостью горизонта является широтой наблюдателя. Напомним, что угол дуги большой окружности сферы равняется центральному углу между радиусами ее концов, например, дуга φ равна ∠HMN.

       Годовое движение Солнца S прочерчивает на небесной сфере 1 круг эклиптики EE' против часовой стрелки. Он пересекает круг экватора AA' в точке γ. Долгота Солнца λ отсчитывается от точки γ, которая является точкой весеннего равноденствия. Расстояние Солнца от экватора AA' определяется склонением δ.

       Вокруг оси MN Земля вращается против часовой стрелки. А небесная сфера совместно с Солнцем вокруг этой оси совершает суточное вращение относительно наблюдателя за часовой стрелкой. Поэтому суточное движение Солнца проходит по кругу SrMdSd, который параллелен кругу экватора. В точке Sr Солнце поднимается над горизонтом, в точку Md оно приходит в полдень, а в т. Sd заходит за горизонт. Часовой угол Солнца ω отсчитывается от точки полудня Md.

       Следует отметить, что на рис. 1 плоскости экватора AA' и эклиптики EE' названы подвижными, т.к. они со временем изменяются в пространстве.

Рис. 1. Основные геометрические характеристики Солнца S при облучении точки M на земной поверхности: 1 - небесная сфера; НН' - плоскость горизонта; N - северный полюс; AA' - плоскость подвижного экватора; ЕЕ' - плоскость подвижной эклиптики, а ε - угол между плоскостями AA' и ЕЕ'; Z - зенит точки M, а z = ∠ZMS - зенитный угол Солнца; Дуга HN = φ - географическая широта точки M; ω = ∠MdNS - часовой угол Солнца, отсчитываемый от полудня Md; δ = SB - склонение Солнца; λ = γS - долгота Солнца.

       Длительность дня пропорциональна длине дуги SrMdSd, а длительность ночи соответствует длине дуги SdMnSr. В показанном положении Солнца S день длиннее ночи. Если Солнце S будет на экваторе в т. γ или γ', то в течение суток оно будет перемещаться по кругу экватора AA'. В этом случае день равен ночи. Если Солнце S' будет находиться в южной части небесной сферы, то его путь под горизонтом будет длиннее пути над горизонтом HH', т.е. ночь будет длиннее дня.

       На рис. 2а показано положение наблюдателя, находящегося на большей широте φ. В этом случае суточный круг MdS движения Солнца S не пересекает горизонт и находится над ним. Мы имеем полярный день. При нахождении Солнца S' в южной части сферы суточный круг его движения, проходящий параллельно экватору, будет находиться под горизонтом. Для наблюдателя M наступает полярная ночь. При нахождении Солнца на круге ЕЕ' в положениях более близких к т. γ или γ' наблюдатель M будет иметь и ночь и день.

       В летнее время при нахождении наблюдателя на северном полюсе, φ = 90° (см. рис. 2б), суточное движение Солнца происходит по кругу, параллельному плоскости горизонта НН'. В этом случае круг горизонта на небесной сфере совпадает с кругом экватора AA'. На полюсе нет восходов Солнца (т. Sr) и его заходов (т. Sd). В зимнее время суточное движение Солнца S' на рис. 2б происходит под горизонтом, т.е. на полюсе стоит полярная ночь.

       При нахождении наблюдателя на экваторе φ = 0° (см. рис. 2в) северный полюс N лежит в плоскости горизонта HH', а точка зенита Z находится на плоскости экватора AA'. Суточное движение Солнца происходит по кругу SrMdSd, который параллелен кругу экватора AA'. В точках восхода Sr и заката Sd Солнце движется перпендикулярно горизонту, поэтому день начинается и заканчивается практически мгновенно.

Рис. 2. Геометрические характеристики Солнца S при нахождении наблюдателя M на разных широтах φ земной поверхности: а - в заполярных широтах; б - на Северном полюсе; в - на экваторе: 1 - небесная сфера; соответствующими стрелками показано годичное и суточное направление движения Солнца на небесной сфере; остальные обозначения см. рис. 1.

       Угол падения солнечных лучей z на плоскость горизонта HH (см. рис. 1) отсчитывается от линии зенита ZM. В полдень Md он наименьший z_min, а в точке восхода Sr и заката Sd зенитный угол z = π/2. Изменение этого угла в диапазоне от π/2 до z_min и от z_min до π/2 определяет изменение солнечной радиации в течение светового дня. Тогда количество солнечного тепла на 1 м² земной поверхности в единицу времени, т.е. мощность инсоляции запишется в виде [19] - [20]:

(1)

где J₀ = 1366.22 Вт/м² - поток солнечного тепла на расстоянии Солнца от Земли r, равным среднему радиусу земной орбиты a;

ρ = r/a – относительное расстояние до Солнца.

       Величина J0 называется солнечной постоянной, а зависимость (1) определяет инсоляцию Земли без учета влияния атмосферы. Если в формуле (1) зенитный угол z выразить через другие углы и солнечный поток dW/dt проинтегрировать за время длительности дня, то суточная инсоляция Земли запишется в виде:

(2)

где τ = 24×3600 – длительность суток в секундах, а часовой угол границы дня ω₀, определяемый дугами SrMd = SdMd (см. рис.1), зависит от углов φ и δ и рассчитывается так [19] - [20]:

(3)

       В выражения (2) и (3) для суточной инсоляции входят долгота λ Солнца и его расстояние r от Земли. Они зависят от параметров орбиты Земли (см. рис. 3). Земля (E) движется по орбите против часовой стрелки. В точке перигелия P_E она находится на ближайшем расстоянии R_p от Солнца, а в точке афелия A_phE – на наибольшем R_a. Эксцентриситет орбиты определяется этими расстояниями

(4)

Рис. 3. Схема движения Земли (E) по орбите вокруг Солнца (S): γ – точка весеннего равноденствия; P_E и A_phE – перигелий и афелий орбиты Земли, соответственно; φ_o – полярный угол движения Земли по орбите; φ_pγ – угол перигелия Земли.

       На рис. 1 плоскость орбиты Земли расположена в плоскости эклиптики EE', а линия Sγ орбиты совпадает с линией Mγ. По кругу эклиптики EE' перемещается Солнце S относительно Земли, находящейся в т. М. Поэтому на рис. 1 изображения Земли будут зеркальными Солнцу. В точке γ весной находится Солнце, а Земля в этот момент находится в т. γ', если в т. М расположено Солнце S. Плоскости экватора и орбиты Земли изменяются в пространстве, поэтому точка γ по орбите Земли (рис. 3) перемещается. Положение перигелия P_E независимо от точки Sγ также перемещается по орбите. Угол φ_pγ измеряется между этими двумя точками γ и P_E.

       На основании новых представлений об орбитальном движении [21] разработан алгоритм расчета параметров r и λ, от которых зависит суточная инсоляция (2). Была создана программа для ее определения в любой день года и для расчета инсоляции за сезоны, полугодия, за год, а также и др. ее составляющие [19] - [20]. На рис. 4 дано изменение по широте φ° летней Q_S, зимней Q_W и уменьшенной в 2 раза годовой Q_T инсоляции в современную эпоху. Инсоляция за год Q_T от экватора монотонно убывает к полюсам. При этом на полюсах она в 2.4 раза меньше, чем на экваторе. Зимняя инсоляция Q_W вблизи экватора имеет максимальное значение, а на полюсах стремится к нулю. Летняя инсоляция Q_S имеет максимальные значения вблизи тропиков (φ = ± ε° = ± 23.4°) и наименьшие значения на полюсах. При этом в современную эпоху летняя инсоляция Q_S на Северном полюсе, как видно из графика на рис. 4, меньше, чем на Южном. Также летняя инсоляция на северном тропике в 1.04 раза меньше, чем на южном.

Рис. 4. Распределение по широте Земли удельного тепла в ГДж/м² в современную эпоху (1950 г.): Q_S – за летнее калорическое полугодие; Q_W – за зимнее калорическое полугодие; Q_T – за весь год: на графике величина Q_T уменьшена в два раза; φ > 0 - северное полушарие; φ < 0 - южное полушарие; Mil - расчеты М. Миланковича.

       Программа для расчета всех составляющих инсоляции имеется в свободном доступе [22]. Она основана на нашем алгоритме расчета инсоляции, который более удобен для применения, чем алгоритм М. Миланковича. Количественно результаты этих двух теорий, как показано на рис. 4, совпадают.

Рис. 5. Параметры орбиты и оси Земли в неподвижных экваториальной xyz и эклиптической x₁y₁z₁ системах координат.

1 – небесная сфера; неподвижные плоскости в эпоху Т₀: 2 – экватора Земли, 3 – орбиты Земли (плоскость эклиптики); подвижные плоскости в эпоху Т: 4 – экватора Земли, 5 – орбиты Земли; единичные вектора: – оси Земли, – оси орбиты Земли; – вектор момента количества движения Солнечной системы; γ₀ – точка весеннего равноденствия в эпоху Т₀; B – положение перигелия на небесной сфере; φ_ω = γ₀γ₂ – угловое расстояние восходящего узла орбиты; φ_р = γ₂B – угловое расстояние перигелия; i – угол наклона плоскости орбиты к плоскости неподвижного экватора.

Угловые параметры и расстояние ρ, входящие в формулу суточной инсоляции (2) зависят от параметров орбиты Земли и оси ее вращения, и они, как уже отмечалось, изменяются со временем. На рис. 5 кругами на небесной сфере 1 показаны плоскости экватора Земли 2 и ее орбиты 3, застабилизированные на определенную эпоху, например 2000 года. Угол между этими плоскостями равен ε₀. В другую эпоху плоскость экватора переместится в положение 4, а плоскость орбиты – в положение 5, и угол между ними будет ε. Для определения законов изменения параметров, входящих в инсоляцию (2), необходимо решить две задачи об изменении: 1) орбитального и 2) вращательного движений.

       Согласно закону всемирного тяготения тело с номером k притягивает тело с номером i, и это притяжение выражается силой:

(5)

где G – гравитационная постоянная; – радиус-вектор от тела с массой m_k до тела с массой m_i.

       На i-тое тело все остальные k тела оказывают воздействие силой (5). Тогда просуммировав силы всех воздействий, и поделив на массу i-того тела, получаем его ускорение в виде:

(6)

где – радиус-вектор тела m_i относительно некоторого центра в неускоренной системе координат (в данном случае относительно центра масс Солнечной системы).

       Соотношение (6) представляет систему 3n нелинейных дифференциальных уравнений, где n = 11 (девять планет, Солнце и Луна). Для их решения разработан метод высокой точности, на основании которого создана система Galactica [23], которая имеется в свободном доступе [24]. Эффективность метода показана в результате решения ряда задач современной небесной и космической динамики [25]-[29].

       Эволюция орбиты рассматривается в неподвижной системе координат xyz, связанной с неподвижной плоскостью экватора 2 (см. рис. 5). Исследуется эволюция эксцентриситета орбиты e, угла ее наклона i к плоскости экватора 2, угла φ_Ω = γ₀γ₂, восходящего узла орбиты γ₂ и углового положения перигелия φ_р = γ₂B, где B – проекция перигелия на небесную сферу 1.

       На рис. 6 точками представлена динамика изменения элементов орбиты Земли: эксцентриситета e и углов φ_Ω, i и φ_р. На промежутке 7 тысячелетий (от -3.4 тыс. лет до 3.6 тыс. лет) эксцентриситет е и наклон орбиты i уменьшаются, а угол перигелия φ_р увеличивается. В то время как угол восходящего узла φ_Ω вначале этого интервала уменьшается, а затем возрастает. При этом минимум φ_Ω наступает за 1000 лет до 30.12.49 г. Угол перигелия φ_р в отличие от других параметров изменяется неравномерно. Линиями 2 и 3 на рис. 6 приведены средние изменения этих параметров, которые получены С. Ньюкомбом [30] и Дж. Л. Симоном и др. [31] по наблюдениям. Из графиков видно, что результаты расчетов подтверждаются наблюдениями.

Рис. 6. Вековые изменения орбиты Земли 1 и сопоставления с аппроксимациями 2 и 3 данных наблюдения С. Ньюкомба [30] и Дж.Л. Симона и др. [31], соответственно: e – эксцентриситет; i – наклон плоскости орбиты к плоскости экватора 2000.0 г.; φ_Ω – угловое положение восходящего узла орбиты от оси x на эпоху 2000.0 г.; φ_р – угловое положение перигелия в плоскости орбиты от восходящего узла; углы даны в радианах, а время Т – в столетиях от 30.12.1949 г.; интервал между точками – 200 лет; 1 cyr – столетие.

       Эволюция параметров орбиты Земли за 3 прошедших млн. лет показана на рис. 7. Эксцентриситет испытывает короткопериодические изменения со средним периодом Т_е1 = 94.5 тыс. лет вокруг его среднего значения е_m = 0.028. Кроме того, наблюдаются более длительные колебания с периодом Т_е2 = 413 kyr, которые приводят к достижению крайних значений эксцентриситета е = 0.0022 и е = 0.062. Долгота восходящего узла φ_Ω меняется со средним периодом Т_Ω = 68.7 тыс. лет вокруг среднего значения φ_Ωm = 0.068 радиан. На основной период накладываются колебания с большей длительностью.

       Угол наклона плоскости орбиты і испытывает колебания с таким же периодом Т_і = 68.7 тыс. лет вокруг среднего значения і_m = 0.402 радиан. Колебания угла i происходят в пределах от 0.36 до 0.45 радиан. Диапазон колебаний составляет 5°.

       Угол положения перигелия φ_р увеличивается со временем. Перигелий перемещается в направлении обращения Земли вокруг Солнца, совершая в среднем за 3 млн. лет один оборот за 147 тыс. лет. Вращение перигелия происходит неравномерно. Наряду со средним вращением против часовой стрелки существуют его возвратные движения по часовой стрелке.

Рис. 7. Эволюция орбиты Земли за 3 млн. лет в прошлое. Обозначения те же, что и на рис. 6. Величины с индексом m – средние значения параметров за 50 млн. лет; Т_e, Т_i, Т_Ω – основные периоды колебаний соответствующих параметров в тысячах и миллионах лет, а Т_p – средний за 3 млн. лет период обращения перигелия; интервал между точками – 10 тыс. лет; kyr – тысячелетие; Myr – 1 миллион лет.

       В результате анализа полученных решений было установлено, что ось орбиты Земли (см. рис. 5) вращается за стрелкой часов, т.е. против орбитального движения Земли с периодом T_S = -68.7 тыс. лет. Напомним, что ось орбиты – это перпендикуляр к ее плоскости 5. Благодаря этому вращению происходят колебания углов i и φ_Ω, показанных на рис. 7. Вращение оси орбиты происходит вокруг вектора (рис. 5), который является моментом количества движения всей Солнечной системы.

       Дифференциальные уравнения орбитального движения (6) были проинтегрированы с помощью программы Galactica за 100 млн. лет. На этом интервале представить результаты в прежних параметрах невозможно, т.к. их колебания сливаются в сплошную полосу. Поэтому на рис. 8 показано эволюция параметров орбиты Земли за 50 млн. лет в несколько ином виде. Скользящие средние эксцентриситета e_s осреднялись на интервале, равном 2Te₂ ≈ 850 kyr. Как видно из графика здесь наблюдается третий период изменения эксцентриситета Te₃ = 2.31 Myr. Вместо углов орбиты i и φ_Ω на рис. 8 показана эволюция углов Θ_Ss и ψ_S оси орбиты относительно вектора (см. рис. 5). Угол наклона Θ_Ss, также представленный в виде скользящего среднего, является углом между векторами и .

Рис. 8. Эволюция параметров орбиты Земли за 50 млн. лет впрошлое: e_s – скользящие средние эксцентриситета; φ_р – угол перигелия; ψ_s – угол прецессии и Θ_Ss – скользящие средние угла наклона оси орбиты . Углы φ_р и ψ_s даны в радианах, Θ_Ss – в градусах.

       Максимальное отклонение оси от момента составляет Θ_Smax = 2.94°. Периоды ее колебаний следующие: Т_Θ1 = 97.4 kyr, Т_Θ2 = 1.16 Myr и Т_Θ3 = 2.32 Myr. Колебания с этими периодами видны на рис. 8. Угол прецессии ψ_s идентичен углу φ_Ω, но только измеряется не в плоскости экватора 2, а в плоскости, перпендикулярной моменту количества движения . Как и угол перигелия φ_р, угол прецессии ψ_S изменяется неравномерно, и существуют его возвратные движения. Однако в таком масштабе, как показано на рис. 8, этих колебаний φ_р и ψ_S не видно. Вследствие неравномерности вращения средние периоды перигелия и оси орбиты на интервале -50 млн. лет немного другие и равны: T_p = 150 тыс. лет, и T_S = -72.9 тыс. лет, соответственно.

       На интервале от -50 млн. лет до -100 млн. лет изменение параметров не отличается от рис. 8 [32]. Аналогичные исследования были выполнены для остальных планет и Луны [33]-[34]. Параметры их орбит так же, как и Земли совершают установившиеся колебательные движения. Поэтому можно сделать вывод, что, кроме установленных периодов и амплитуд колебания параметров орбит, других не существует. Из полученных результатов также следует вывод об устойчивости Солнечной системы. Необходимо отметить, что в работах других исследователей [35]-[37], которые решали задачу другими методами, были получены изменяющиеся колебания. Поэтому они пришли к выводу о неустойчивости орбит и Солнечной системы в целом [37].

Эволюция вращательного движения Земли

       Вращающаяся Земля под действием центробежных сил растягивается в плоскости экватора (см. рис. 9). Рассмотрим воздействие тела B на две половины Земли: на ближнюю – силой и дальнюю – силой . Если бы Земля была центрально-симметричным шаром, то равнодействующая сил и прошла бы через центр Земли O. Для сплюснутой Земли центр масс ближней части Земли приблизится к телу B, а дальней – удалится. Поэтому сила увеличится, а сила уменьшится, в результате чего возникнет момент сил m_O, направленный по часовой стрелке. Этот процесс воздействия на Землю определяется теоремой изменения момента количества движения:

(7)

где   момент количества движения Земли относительно центра О в невращающейся системе координат x₁y₁z₁, – сумма моментов сил, действующих на Землю тел.

Рис. 9. Прецессия оси Земли при воздействии на нее тела B: 1 и 2   плоскости экватора Земли и ее орбиты, соответственно; 3 – плоскость орбиты воздействующего на Землю тела B.

       Задача о вращении Земли рассматривается (см. рис. 5) в невращающейся системе координат x₁y₁z₁, связанной с неподвижной плоскостью орбиты Земли 3. Подвижная плоскость экватора Земли 4 определяется углом наклона θ к плоскости 3 и углом прецессии ψ = γ₀γ₁. Кроме того, рассматривается скорость собственного вращения Земли относительно ее подвижной оси . На основании теоремы (7) мы вывели [38]-[39] дифференциальные уравнения вращательного движения Земли в следующем виде:

(8)

(9)

(10)

где J_x, J_y и J_z – моменты инерции вращающейся Земли на оси системы координат, связанные с вращающейся Землей;

E_d = (J_z-J_x)/J_z – динамическая эллиптичность Земли;

ω_E = const – проекция абсолютной скорости вращения Земли на ее ось (см. рис. 5);

n – количество действующих на Землю тел, а M_i – их масса, а x_1i, y_1i, z_1i – их координаты.

       Уравнениями (8) и (9) определяются углы положения ψ и θ подвижной плоскости экватора 4 относительно плоскости неподвижной орбиты 3 (см. рис. 5). В результате интегрирования уравнений (8)-(9) определяются углы наклона θ и прецессии ψ подвижного экватора 4 относительно неподвижной плоскости орбиты 3. А в результате интегрирования уравнений орбитального движения (6) рассчитываются углы наклона i и восходящего узла φ_Ω подвижной плоскости орбиты 5 относительно неподвижной плоскости экватора 2, а также угол перигелия φ_p. Эти параметры вращательного и орбитального движения позволяют определить углы наклона ε и перигелия φ_pγ подвижной плоскости орбиты 5 относительно подвижного экватора 4. Именно эволюция этих углов ε, φ_pγ и эксцентриситета e определяют эволюцию инсоляции.

       Рассмотрим основные результаты задачи (8)-(9) о вращении Земли. На рис. 10 показано изменения угла наклона ε на пяти разных интервалах времени I_n. На графиках указаны основные периоды T_ni и амплитуды (θ_ai и ε_a4) колебаний угла наклона: полумесячные T_n2, полугодовые T_n3 и T_n4 = 18.6 лет. Эти колебания называются нутационными. Угол прецессии ψ имеет аналогичные периоды колебаний, а амплитуды – в 2 - 3 раза большие.

       На интервале I_n = 0.1 года видны полумесячные колебания и просматриваются суточные; на интервале I_n = 1 год появляются полугодовые колебания; на интервале I_n = 10 лет появился тренд колебания с периодом 18,6 года, а на интервале I_n = 100 лет колебания с этим периодом преобладают.

       На интервале I_n = 100 лет видно, что рассчитанный 1 угол наклона ε колеблется вокруг среднего 2 угла наклона по С. Ньюкомбу [30] и Дж. Симону и др. [31]. Амплитуда колебаний ε_a4 = 9.2" периода T_n4 = 18.6 лет также совпадает с наблюдениями. Эта амплитуда в астрономии называется постоянной нутации. Рассчитанный угол прецессии ψ также колеблются относительно осредненного угла прецессии по наблюдениям, а средняя динамика ψ также совпадает с наблюдениями.

Рис. 10. Динамика угла наклона θ и ε (в радианах) плоскости экватора Земли к плоскости ее орбиты на пяти интервалах In времени: yr – год; Δθ ≈ ε - ε0; ε0 – угол наклона в начальную эпоху 30.12.1949 г. Tn2, Tn3, Tn4 и θa2, θa3, εa4 – периоды и амплитуды колебаний угла наклона ε; 1 – по нашим результатам численного интегрирования; 2 – аппроксимация данных наблюдений по С. Ньюкомбу [30] и Дж. Симону и др. [31]; 3 – по результатам интегрирования Дж. Ляскара и др. [12]; 4 – по результатам интегрирования Ш. Г. Шараф и Н. А. Будниковой [10].

Как видно из рис. 10, на интервале 10 тыс. лет такое совпадения с аппроксимацией данных наблюдения и результатами решений других авторов [10] и [12] имеется до 2000 лет. Далее рассчитанный нами угол наклона ε начинает отличаться от результатов решений [10] и [12].

       На рис. 11 показано, что со временем отличия нарастают, и дальнейшая эволюция рассчитанного нами угла наклона ε значительно отличается от эволюции, полученной другими авторами, которые решали упрощенную задачу о вращении Земли. Как видно из графиков, колебания угла ε по нашим решениям происходят в диапазоне от 16.7° до 31°, в то время как по прежним решениям от 22.26° до 24.32°, т. е. диапазон колебаний больше в 7 раз. Мы выполнили расчеты инсоляции Земли с новой эволюцией угла ε [17]. Ее колебания также в 7 раз превышают колебания инсоляции в прежних астрономических теориях палеоклимата.

       Второй угол в уравнениях (8) - (10), угол прецессии ψ, изменяется несколько по другому, по сравнению с углом θ. Он непрерывно уменьшается, но это уменьшение совершается с такими же колебаниями, как и угла θ, только амплитуда их в 2-3 раза больше. Такое изменение угла ψ свидетельствует, что ось Земли (см. рис. 5) вращается за стрелкой часов, и средний период ее вращения P_N = -25740 лет.

       Выполненные исследования по этим двум задачам позволяют представить следующую картину движений в Солнечной системе. В плоскости орбиты происходит неравномерное вращение перигелия против стрелки часов со средним периодом T_p = 147 kyr и деформация орбиты с колебаниями эксцентриситета с периодами 94.5 kyr, 413 kyr, 2.31 Myr. Аналогичные изменения орбиты со своими периодами испытывают Луна и другие планеты, за исключением Плутона. Перигелий Плутона вращается по часовой стрелке.

Рис. 11. Эволюция угла наклона ε (в радианах) плоскости экватора Земли к плоскости ее орбиты на интервале 200 тыс. лет: 1 – по нашим результатам численного интегрирования; 2 – по результатам интегрирования Дж. Ляскара и др. [12]; 3 – по результатам интегрирования Ш. Г. Шараф и Н. А. Будниковой [10]. В градусах приведены максимальные и минимальные значения угла ε.

       Ось орбиты Земли (см. рис. 5) вращается по часовой стрелке вокруг вектора момента со средним за 3 млн. лет периодом T_S = 68.7 kyr. Ось также совершает колебания с разными периодами 97.4 kyr, 1.16 Myr и 2.32 Myr. Ось орбиты Луны (на рис. 5 не показана) вращается по часовой стрелке вокруг подвижной оси орбиты Земли. Период этого вращения равен 18.6 лет. Кроме того, ось орбиты Луны совершает колебательное движение с периодом 0.4745 года.

       Ось вращения Земли (см. рис. 5), как уже отмечалось, вращается по часовой стрелке, как и ось орбиты Луны, вокруг подвижной оси орбиты Земли со средним периодом 25740 лет. При этом ось совершает колебания с периодами: полумесячным, полугодовым и 18.6 лет. Кроме того, существуют колебания оси большой амплитуды (см. рис. 11) и с периодами десятки и сотни тысяч лет. В процессе этих колебаний угол наклона ε оси Земли изменяется от 16.7° до 31°.

6. Эволюция инсоляции

       На рис. 10 и рис. 11 было представлено изменение угла наклона оси Земли в будущем. На рис. 12 линией 1 показана эволюция угла наклона ε за прошедшие 200 тыс. лет. По сравнению с расчетами по прежним теориям (линия 2) наши решения так же, как при расчетах в будущее, дают в 7 раз большие амплитуды колебаний. В начальный период несколько тысяч лет в прошлое, начиная с T = 0, угол наклона 1, как и угол наклона 2, увеличивается. Затем он начинает уменьшаться и достигает минимума, в то время как угол наклона 2 по прежним теориям находится в максимуме. На остальном интервале времени также максимумы и минимумы колебаний углов наклона 1 и 2 не совпадают. Однако более существенным являются величины этих экстремумов. По прежним теориям на этом интервале времени угол наклона изменяется от 22.21° до 24.43°. А по нашим решениям наклон плоскости экватора Земли к плоскости ее орбиты изменяется от 14.8° до 32.1°. Такой же диапазон изменений угла ε получен при решении задачи за будущие 200 тыс. лет (см. рис. 11).

       В астрономических теориях палеоклимата рассматривается инсоляция не за астрономические, а за равные калорические полугодия. Начало и конец летнего калорического полугодия определяется так, чтобы инсоляция за любой его день была больше инсоляции за любой день зимнего полугодия. Далее будем рассматривать инсоляцию на широте 65° северного полушария, которое обозначается буквой N. Мы рассчитали изменение инсоляции Q_S^65N за 200 тыс. лет в прошлое как по определенным нами параметрам: e, ε и φ_pγ (линия 1 на рис. 12), так и по тем же параметрам, рассчитанных Дж. Ляскаром и др. [12] (линия 2). Как видно из графиков, инсоляция Q_S^65N за летнее калорическое полугодие на северной широте 65° по нашим решениям имеет амплитуду колебаний также в 7 раз большую, нежели по прежним теориям. Моменты потепления и похолодания по нашим расчетам 1 и по прежним теориям 2 также не совпадают. Вначале, начиная с T = 0, как видно из Q_S^65N на рис. 12, летняя инсоляция растет в течение 4-х – 5-и тысяч лет, затем начинает падать и достигает минимума 16 тысяч лет назад. После этого минимума идет потепление, которое завершается большим максимумом инсоляции 31 тыс. лет назад.

       Средний период колебаний угла наклона по прежним теориям (см. 2 на рис. 12) равен 41.1 тыс. лет. Этот же период имеется в колебании инсоляции Q_S^65N (линия 2). Период колебаний 41.1 тыс. лет является наиболее упоминаемый в астрономической теории палеоклимата. Как видно из новой зависимости для угла наклона ε (линия 1) характерный период колебаний имеет в 1.5 - 2 раза меньшую длительность.

Рис. 12. Эволюция угла наклона ε, а также летних инсоляций Q_S^65N и I за 200 тыс. лет в прошлое. Сравнение новых результатов 1 с результатами прежних теорий 2 на примере работы Дж. Ляскара и др. [12]. ε – угол наклона в радианах экватора Земли к плоскости ее орбиты; Q_S^65N – инсоляция в ГДж/м² за летнее калорическое полугодие на северной широте 65°; I° – инсоляция в эквивалентных широтах за летнее калорическое полугодие на северной широте 65°. В градусах приведены максимальные и минимальные значения угла ε.

       По нашим решениям получены в 7 раз большие колебания инсоляции. Насколько они существенны? На этот вопрос может дать ответ представление инсоляции в эквивалентных широтах I, которая рассчитывается следующим образом. Если в эпоху T летняя инсоляция на широте φ была такая, как в современную эпоху на широте φ₀, то инсоляция в эквивалентных широтах будет I = φ₀. На рис. 12 представлена инсоляция I в эквивалентных широтах, рассчитанная нами для широты 65° северного полушария как по нашим данным (линия 1) так и по данным Дж. Ляскара и др. [12] (линия 2). Начиная с T = 0, инсоляция I по нашим данным 1 с широты 65° уменьшается на несколько градусов, т.е. на широте 65° становится теплее. Затем после максимума величина I начинает уменьшаться до широт 80° и 90°. При T = -15 тыс. лет летняя инсоляция на широте 65° меньше, чем современная летняя инсоляция на полюсе, поэтому на графике I проведена горизонтальная площадка. Таким образом, горизонтальная площадка примерно 15 тыс. лет назад означает, что инсоляция на широте 65° была меньше, чем сейчас на полюсе. Такое малое количество тепла могло привести к оледенению территорий на широте 65°.

       При дальнейшем изменении времени до T = -30 kyr инсоляция I в эквивалентных широтах достигает широт 50°, 40°, 30°, т.е. на широте 65° становится значительно больше солнечного тепла. Горизонтальная площадка в эпоху T = -30 kyr означает, что на широте 65° поступает больше тепла, чем сейчас на экваторе.

       Линией 2 приведена инсоляция I в эквивалентных широтах по прежним теориям. Как видно, на широте 65° летняя инсоляция I на рассмотренном отрезке времени в 30 - 50 тыс. лет изменяется в пределах от 60° до 70°. Маловероятно, что изменения количества тепла на широте 65° до таких величин, которые сейчас имеются на широтах 60° и 70°, могут привести к существенным потеплениям климата или его похолоданиям. Такие небольшие изменения инсоляции всегда вызывали сомнения [7].

Рис. 13. Распределение по широте φ° Земли удельного количества тепла в ГДж/м² в самую жаркую в высоких широтах эпоху T = 31.28 тыс. лет в прошлое: Q_s – за летнее калорическое полугодие; Q_w – за зимнее калорическое полугодие; Q_T – за весь год: на графике величина Q_T уменьшена в два раза; φ > 0 – северное полушарие; φ < 0 – южное полушарие.

       Выше мы исследовали изменение инсоляции во времени на широте 65° северного полушария. Теперь рассмотрим изменение инсоляции по широте в отдельные моменты времени. На рис. 4 была представлено распределение инсоляции по всему земному шару в современную эпоху. На рис. 13 показано изменение по широте φ° летней Q_s, зимней Q_w и уменьшенной в два раза годовой инсоляции Q_T в эпоху T = -31.28 тыс. лет. В это время на широте 65° северного полушария была наибольшая за 200 тыс. лет летняя инсоляция (см. рис. 12). Как и в современную эпоху (см. рис. 4), инсоляция за год Q_T – наибольшая на экваторе и убывает к полюсам. Зимняя инсоляция Q_w имеет максимальное значение вблизи экватора, а на полюсах стремится к нулю. Летняя инсоляция Q_s имеет максимумы на тропиках, которые находятся на больших широтах по сравнению с современной эпохой (см. рис. 4). Кроме того, летняя инсоляция Q_s на высоких широтах близка к максимальной и существенно превышает летнюю инсоляцию на экваторе.

       Инсоляция за год Q_T в эпоху T = -31.28 тыс. лет на высоких широтах больше (рис. 13), чем в современную (рис. 4). В то же время годовая инсоляция на экваторе меньше, чем в современную эпоху. Но это уменьшение не такое значительное, как увеличение в высоких широтах. Зимняя инсоляция Q_w в эпоху T = -31.28 на всех широтах меньше по сравнению с современной эпохой (рис. 4).

Рис. 14. Распределение по широте φ° Земли удельного количества тепла в ГДж/м² в самую холодную в высоких широтах эпоху T = 46.44 тыс. лет в прошлое. Обозн. см рис. 13.

       На рис. 14 представлено изменение этих же составляющих инсоляции в наиболее холодную за 200 тыс. лет эпоху T = -46.44 тыс. лет. В высоких широтах, например, на широте φ° = 65°N летняя инсоляция Q_s = 4.72 ГДж/м² существенно уменьшилась по сравнению с Q_s = 7.43 ГДж/м² в эпоху T = -31.28 тыс. лет и Q_s = 5.92 ГДж/м² в современную эпоху T = 0. Инсоляция за год, например Q_T = 3.58 ГДж/м² на Северном полюсе, уменьшилась с Q_T = 7.43 ГДж/м² в эпоху T = -31.28 тыс. лет и с Q_T = 5.58 ГДж/м² в современную эпоху T = 0. Зимняя инсоляция Q_w в эту холодную эпоху немного возросла на всех широтах.

       Итак, от теплой эпохи к холодной уменьшение летней инсоляции Q_s на широте 65° происходит в 1.57 раза. А на полюсе годовая инсоляция Q_T уменьшается еще сильнее, а именно – в 2 раза. В то же время в экваториальной зоне происходят обратные по знаку изменения, но величина их в 4 раза меньшая. При этом на широте φ = 45° годовая инсоляция Q_T практически не изменяется, т.е. изменения летней инсоляции Q_s компенсируются изменениями зимней инсоляции Q_w.

       Полученные в наших расчетах колебания инсоляции могут вызвать наблюдавшиеся изменения климата. Уменьшение летней инсоляции I 19 - 12 тыс. лет назад до значений меньших, чем на полюсе (линия 1 на рис. 12), могло привести к оледенению. Это, по-видимому, подтверждается существованием в это время покровного ледника в северных частях Евразии и Америки. А увеличение летней инсоляции I 35 - 28 тыс. лет назад до значений больших, чем в экваториальной области, могло привести к благоприятным условиям существования мамонтовой фауны. По-видимому, это подтверждается преимущественной датировкой ископаемых представителей этой фауны.

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

7. Проверка достоверности решения задачи о вращении Земли

       Как уже отмечалось во Введении, все три задачи: инсоляция Земли, ее орбитальное движение и вращательное движение Земли мы решали по-другому, начиная с вывода уравнений и заканчивая обработкой и анализом результатов. Например, дифференциальные уравнения (6) для орбитальной задачи ранее непосредственно не решали. Начиная с И. Ньютона, эту задачу рассматривали в рамках так называемой теории возмущений. Исследовалось движение одного тела под воздействием материнского тела: Солнца – для планет, и Земли – для Луны, а воздействия остальных тел рассматривались как возмущающиеся факторы. В итоге, в результате работы поколений таких выдающихся ученых как Эйлер, Лагранж, Лаплас и мн. др. были получены уравнения возмущенного движения для шести элементов орбиты относительно неподвижной эклиптики (см. 3 на рис. 5). Эти уравнения приближенные. А затем они также решались приближенными аналитическими методами разложением в ряды. В настоящее время эти методы переведены на численный счет, благодаря чему количество членов в рядах достигает несколько сотен. Тем не менее, приближенный характер уравнений при решении задач за сотни тысяч и миллионы лет всегда оставляет сомнение в справедливости результатов. Поэтому полученные нами решения неупрощенной орбитальной задачи (6) позволили подтвердить многие результаты, полученные прежними методами, а некоторые результаты, как неустойчивость орбит и Солнечной системы в целом, были опровергнуты.

       Задача о вращении Земли (8)-(10), при ее решения за большие промежутки времени в работах прежних авторов упрощалась более существенно, по сравнению с орбитальной. В уравнениях (8)-(9) отбрасывались вторые производные и , а также отбрасывались произведения первых производных и . Поэтому такое решение задачи не позволяло получить короткопериодических колебаний, которые получены нами и показаны на рис. 10. А полученные колебания угла ε с периодом 41 тыс. лет (см. линия 2 на рис. 11 и на рис. 12) не имели других подтверждений.

       Задачу (8)-(10) о вращении Земли за сотни тысяч лет мы решили впервые. В рамках принятой технологии ее решения все необходимые проверки были выполнены. Например, задача решалась последовательно, при воздействии одного из тел [38]. Полученные периоды колебаний оси Земли подтвердились общетеоретическими выводами, а также результатами других авторов [39]. При воздействии всех тел задача решалась за разные интервалы времени, и полученные результаты подтвердились наблюдениями. Интегрирование уравнений (8)-(10) за 200 тыс. лет проводилось с разными начальными условиями и с разными шагами интегрирования. Вид результатов от этого не изменился. Однако полученные амплитуды колебаний оси Земли в 7 раз превышали амплитуды по прежним теориям. Поэтому требовалось подтверждение полученных результатов. По-видимому, нужно решать эту задачу каким-то другим способом. На протяжении ряда лет мы занимались поиском таких способов и их реализацией.

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

       При исследовании орбит выяснилось, что эволюция орбиты Луны подобна эволюции оси вращения Земли. Этот результат привел нас к составной модели Земли (см. рис. 15), в которой часть массы Земли равномерно распределена между периферийными телами, обращающимися вокруг центрального тела по круговой орбите. Под воздействием Луны, Солнца и планет орбиты периферийных тел начинают изменяться. Эволюция оси орбиты одного из этих тел моделирует эволюцию оси вращения Земли. В первоначальной серии исследований [26] было исследовано три модели, и была подтверждена возможность моделирования эволюции Земной оси. В этих моделях периоды прецессии осей орбит моделей были 170 лет и 2604 года, в то время как средний период прецессии оси Земли T_prE = 25740 лет. В последующем было разработано еще 11 моделей, пока не был достигнут необходимый период прецессии.

Рис. 15. Составная модель вращения Земли. Масса Земли распределена между центральным телом и периферийными телами: a - радиус орбиты периферийных тел.

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

       На рис. 16 точками представлены результаты 13-ой составной модели Земли на интервале 300 лет. Видно, что угол наклона ε и угол прецессии ψ колеблются с периодом 18.6 лет. Решения на меньших интервалах времени дали колебания с полумесячными и полугодовыми периодами, т.е. эти результаты совпали с результатами интегрирования уравнений (8)-(10). Амплитуды этих колебаний также совпали. Прямыми линиями на рис. 16 нанесены средние изменения углов ε и ψ. Они также совпали с результатами прямой задачи о вращении Земли, а также с данными наблюдений. Идентичность результатов видна также при сопоставлении рис. 16 с графиком для угла ε на интервале 100 лет на рис. 10.

       Такое совпадение результатов модельной задачи с результатами прямой происходит до 3 тыс. лет. В процессе интегрирования уравнений (6) программой Galactica происходит накопление погрешностей для тел модели, вследствие чего ее размеры начинают изменяться. Изменяется также период прецессии орбиты. Например, к концу интервала интегрирования 13.763 тыс. лет период прецессии с 25740 лет уменьшился до 14840 лет. Поэтому модель перестает представлять эволюцию Земной оси. Этот процесс обусловлен весьма напряженными динамическими параметрами модели. Например, радиус орбиты периферийных тел (см. рис. 15) a равен радиусу Земли, их период обращения 0.142 часа, а взаимодействие между телами модели усилено в 9.6 раза по сравнению с гравитационным. Таким образом, тела модели в 170 раз вращаются быстрее Земли. Поэтому шаг интегрирования задачи (6) приходится уменьшить в 1000 раз по сравнению с шагом при решении орбитальной задачи. Это приводит к большому времени счета. Например, решение задачи за 13.763 тыс. лет заняло 2.13 месяца. Чтобы модель не изменялась на интервале 200 тыс. лет, нужно уменьшить шаг интегрирования до таких величин, при которых требуется невыполнимое время счета.

Рис. 16. Эволюция углов наклона ε и прецессии ψ оси вращения составной модели Земли № 13 за 300 лет. Точками приведены результаты интегрирования уравнений (6) с помощью программы Galactica. Расстояние между точками 3 года. Линиями показаны средние изменения углов угла ε и ψ со скоростями = - 2.28·10^-4 1/cyr = - 0.470 "/год. = - 2.42·10^-2 1/cyr = - 49.9 "/год.

       Итак, на данном этапе составная модель Земли на интервале 3000 лет подтвердила полученные результаты интегрирования дифференциальных уравнений (8)-(10) вращения Земли. Это свидетельствует, что принятые допущения и упрощения при выводе уравнений (8)-(10), их вывод, метод решения и преобразование результатов интегрирования в конечный вид также подтвердился.

       Вторая независимая проверка заключалась в использовании более точного метода при интегрировании уравнений (8)-(10), а именно метода Рунге-Кутты 8-ого порядка в реализации Дорманда и Принса Хайрер Э., Нёрсетт С.П., Ваннер Г. Решение обыкновенных дифференциальных уравнений. – М.: Мир, 1990. – 513 с. До этого использовался метод интегрирования Рунге-Кутты 4-ого порядка в реализации [41]. Мы использовали его в течение нескольких десятков лет для решения разнообразных задач, и всегда получали удовлетворительные результаты. При интегрировании уравнений вращательного движения (8)-(10) за интервалы времени порядка 200 тыс. лет, мы столкнулись с неожиданным недостатком метода. В решениях этих уравнений существуют суточные колебания производных и . Амплитуда этих колебаний к концу вышеупомянутого интервала интегрирования увеличивалась на несколько порядков. Несмотря на то, что были проведены проверки по влиянию суточных амплитуд и на конечные результаты и приняты определенные меры по их устранению, оставалась опасность непредусмотренного их влияния.

       Интегрирование уравнений (8)-(10) было выполнено новым методом в реализации Дорманда и Принса на разных интервалах времени, включая 200 тыс. лет. Все полученные ранее результаты подтвердились. При этом амплитуды суточных колебаний производных и не увеличиваются и остаются на одном уровне. Итак, метод интегрирования уравнений не влияет на полученные результаты, и применение более точного метода подтверждает их.

       Третья независимая проверка заключалась в смене технологии решения задачи. В дифференциальные уравнения вращательного движения (8)-(9) входят привязанные к плоскости орбиты координаты x_1i, y_1i, z_1i воздействующих на Землю тел. При решении орбитальной задачи (6) с помощью программы Galactica мы получаем координаты тел x_i, y_i, z_i, привязанные к плоскости неподвижного экватора. Они пересчитывались в координаты x_1i, y_1i, z_1i. Однако при интегрировании задачи (8)-(9) за большие промежутки времени массив этих данных будет занимать нереализуемый объем памяти. Поэтому мы создали математическую модель Солнечной системы [21], которая выдает в необходимый момент времени координаты тел: x_1i, y_1i, z_1i на основании результатов задачи двух тел: материнского тела и его спутника. При этом параметры орбиты тела: e, i, φ_Ω, φ_р, R_р и др. в каждый момент определяются по данным, предварительно рассчитанных с помощью программы Galactica. В процессе решения этой задачи математическая модель Солнечной системы была всесторонне проверена. Тем не менее, оставалась вероятность того, что на больших интервалах времени несущественные отличия результатов математической модели Солнечной системы от значений координат, полученных с помощью программы Galactica, могут повлиять на эволюцию параметров вращательного движения ε и ψ.

       В течение работы над задачей о вращении Земли мы неоднократно возвращались к вопросу о смене вышеописанной технологии. Однако все попытки были безуспешными, пока не появилась определенная идея эффективного объединения этих двух задач в одну. В итоге была разработана новая программа совместного решения орбитальной задачи и задачи о вращении Земли. В ней за один шаг по времени происходит решение методом программы Galactica орбитальной задачи (6), затем методом Дорманда и Принса за этот шаг по времени решается задача о вращении Земли (8)-(10). С помощью новой программы были выполнены решения этих двух задач за разные интервалы времени, включая интервал 200 тыс. лет. Все ранее полученные результаты подтвердились. Этой проверкой была также подтверждена математическая модель Солнечной системы на большом интервале времени.

Табл. 1. Сравнение результатов трех методов интегрирования уравнений вращательного движения (8)-(9) за 200 прошедших тыс. лет: РГ-4 – метод Рунге-Кутты 4-ого порядка; ДП-8 – метод Рунге-Кутты 8-ого порядка в реализации Дорманда-Принса; Gal – координаты тел, входящих в уравнения (8)-(9) определяются по программе Galactica, а уравнения решаются методом ДП-8.

       Полученные первоначальным методом на рис. 12 графики полностью повторились при решении задачи двумя последними методами. В табл. 1 приведены количественные сравнения периода прецессии P_pr, и минимального ε_min и максимального ε_max углов наклона с точностью по пятую значащую цифру. С этой точностью первые два метода полностью совпали. Результаты третьего метода, как видно из табл. 1, по периоду прецессии P_pr отличаются в 4-ом знаке, а по углу наклона ε – в 5-ом. Так этот метод является боле точным, то последние величины являются уточнением результатов, полученных первыми двумя методами.

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

       Обратимся еще раз к новым решениям по инсоляции Земли. Как видно из рис. 12, как летняя инсоляция Q_s^65N, так и инсоляция в эквивалентных широтах I по нашим расчетам (линия 1) уже на начальном участке 0 ÷ -50 тыс. лет существенно отличается от инсоляции по прежним теориям (линия 2). Согласно новым решениям ледниковый период мог быть в интервале 19 - 12 тыс. прошедших лет. Ему предшествовал очень теплый период 35 - 28 тыс. лет назад. Эти изменения климата подтверждаются исследованиями Архипова А.С. (1928-1998 гг.) по позднему плейстоцену Западной Сибири [42]. Автор пишет, что "Сартанский гляциокомплекс состоит из морен максимальной гыданской и двух рецессионных: ньяпанской и норильской стадий. Кульминация оледенений приходится на 20 - 18 тыс. лет назад, ньяпанская стадия на 15 - 13 и норильская – на 11.5 - 10.4 тыс. лет назад". А вышеотмеченный теплый период в этой же работе [42] подтверждается словами: "Каргинский межледниковый (мегастадиальный) горизонт имеет радиометрический возраст от 55 (50?) до 23 (22) тыс. лет назад и объединяет морские и аллювиальные, а также стадиальные лохподгортские ледниковые, озерно-ледниковые и субаэральные отложения …". Профессор Архипов А.С. известен обоснованностью своих выводов по палеоклимату. Поэтому приведенные выше его результаты заслуживают пристального внимания, несмотря на то, что они не согласуются с результатами прежней астрономической теории палеоклимата (см. 2 на рис. 12). Как уже отмечалось, на протяжении почти 100 лет палеоклиматологи вынуждены были ориентироваться на результаты прежней астрономической теории климата. Сейчас предстоит большая работа по переосмыслению всего массива данных по палеоклимату, в результате которой будет с большей достоверностью установлена история Земли.

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

       Астрономическая теория ледниковых периодов, как первоначально и предполагали ее создатели Л. Агасис, Ж. Адемар, Дж. Кроль, Р. Харгреаевес [43], М. Миланкович и другие исследователи, может объяснять наблюдавшиеся в прошлом наступления ледниковых периодов и смену их очень теплыми временами. Необходимо продолжать исследования в направлении ее дальнейшего развития. С одной стороны, должны быть продолжены решения задачи о вращении Земли за большие интервалы прошедшего времени. С другой стороны, необходимо в свете полученных колебаний инсоляции Земли проанализировать имеющиеся свидетельства палеоклиматических изменений на предмет их корреляции с колебаниями инсоляции Земли.

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

9. Благодарности

       Эта работа была начата в 1995 г. и на разных ее этапах мне помогали: П.А. Апасев, Н.Ю. Апохин, А.В. Байдин, И.В. Бинеев, В.С. Ботвина, Е.А. Коврижкина, О.И. Кротов, Е.Н. Невидимова, А.А. Павлова, М.Л. Панова, М.А. Пономарев, Е.Ф. Сафина, К.Е. Сеченов, Л.И. Смульский, Я.И. Смульский, И.К. Филатов, М.Э. Чикмарева, И.А. Шаболина. Задачи считались на суперкомпьютерах Института прикладной математики им. М.В. Келдыша РАН (г. Москва) и Сибирского суперкомпьютерного центра СО РАН (г. Новосибирск). Работа выполнялась при поддержке грантов губернатора Тюменской области 2003 и 2004 г.г. и интеграционной программы Президиума РАН № 13 и ОНЗ-11 в 2004 - 2011 гг. На всем протяжении настоящую работу поддерживал директор Института, академик В.П. Мельников. Перечисленным коллегам и коллективам выражаю свою признательность.

Литература

1. Джон Б., Дербишир Э., Янг Г., Фейрбрдж Р., Эндрюс Дж. Зимы нашей планеты: Пер. с англ. / Под ред. Б. Джона. – М.: Мир, 1982. – 336 с.

2. Имбри Дж., Имбри К.П. Тайны ледниковых эпох: Пер. с англ. / Под ред. Г.А. Авсюка .– М.: Прогресс, 1988. – 264 с.

3. Agassiz, L. Etudes sur les glaciers, Neuchatel, 1840.

4. Adhemar J.A. Revolutions de la mer: Deluges Periodiques, Carilian-Goeury et V. Dalmont. Paris, 1842.

5. Croll J. On the physical cause of the change of climate during geological epochs // Philosophical Magazine. 1864, V. 28. P. 121-137.

6. Миланкович М. Математическая климатология и астрономическая теория колебаний климата. – М.–Л.: ГОНТИ, 1939. – 207 с.

7. Большаков В.А., Капица А.П. Уроки развития орбитальной теории палеоклимата // Вестник Российской Академии Наук, 2011, т. 81, № 7, с. 603-612.

8. Brauwer D., Van Woerkom A.J.J. The secular variation of the orbital elements of the principal planets // Astron. Pap. – 1950, Vol. 13, P. 2.

9. Вурком А. Астрономическая теория изменения климата/В кн. Изменение климата. М.: ИИЛ. – 1958.

10. Шараф Ш. Г. и Будникова Н. А. Вековые изменения элементов орбиты Земли и астрономическая теория колебаний климата // Тр. Инст. теоретич. астрономии. – Вып. XIV. – Л.: Наука, 1969 г. – С. 48-109.

11. Berger A. and Loutre M. F. Insolation values for the climate of the last 10 million years // Quaternary Science Reviews. 1991. № 10. Р. 297-317.

12. Laskar J., Robutel P., Joutel F., Gastineau M., Correia A.C.M., and Levrard B. A Long-term numerical solution for the Earth // Icarus 170, 2004. Iss. 2: 343-364.

13. Edvardsson S., Karlsson K.G. and Engholm M. Accurate Spin Axes and Solar System Dynamics: Climatic Variations for the Earth and Mars // Astronomy &Astrophysics, 384, 2002, 689-701. http://dx.doi.org/10.1051/0004-6361:20020029

14. Adem J. Numerical Experiments on Ice Ages climates // Climatic change. 1981, 3. P. 155-171.

15. CLIMAP Project Members: "The Surface of the Ice-Age Earth" // Science, 1976, 191. P. 1131-1137.

16. Lisiecki L.E., Raymo M.E. A Pliocene-Pleistocene Stack of 57 globally distributed benthic δ¹⁸O records // Paleocenography, 2005, 20. P. 1-17.

17. Смульский И.И. Анализ уроков развития астрономической теории палеоклимата // Вестник Российской Академии Наук, 2013. Т. 83. № 1. С. 31-39. http://elibrary.ru/item.asp?id=18448265. doi:10.7868/S0869587313010118

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

19. Смульский И.И., Кротов О.И. Новый алгоритм расчета инсоляции Земли / Институт криосферы Земли СО РАН. – Тюмень, 2013. – 38 с. – Деп. в ВИНИТИ 08.04.2013 № 103-В2013. http://www.ikz.ru/~smulski/Papers/NwAlClI2c.pdf.

20. Smulsky J.J. and Krotov O.I. New Computing Algorithm of the Earth's Insolation // Applied Physics Research, Vol. 6, No. 4; 2014, p. 56-82. http://dx.doi.org/10.5539/apr.v6n4p56.

21. Смульский И.И. Математическая модель Солнечной системы / В сб. Теоретические и прикладные задачи нелинейного анализа. Российская Академия Наук: ВЦ им. А.А. Дородницына. М.: ВЦ РАН А.А. Дородницына, 2007. – С. 119-138. http://www.ikz.ru/~smulski/Papers/MatMdSS5.pdf.

22. Смульский И.И. Программа расчета инсоляции Земли в среде MathCad. / Институт криосферы Земли СО РАН. – Тюмень, 2013. http://www.ikz.ru/~smulski//Data/Insol/.

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

24. Смульский И.И. Система Galactica. Институт криосферы Земли СО РАН. – Тюмень, 2012. http://www.ikz.ru/~smulski/GalactcW/.

25. Смульский И.И. Оптимизация пассивной орбиты с помощью гравиманевра // Космические Исследования, 2008, том 46, No 5. – С. 484-492. http://www.ikz.ru/~smulski/Papers/KOS0484.pdf.

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

27. Smulsky J.J. New Components of the Mercury's Perihelion Precession // Natural Science. Vol. 3, No.4, 268-274 (2011). doi:10.4236/ns.2011.34034. http://www.scirp.org/journal/ns или http://www.ikz.ru/~smulski/Papers/NCmMePNSc.pdf.

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

29. Smulsky J.J. and Smulsky Ya.J. Dynamic Problems of the Planets and Asteroids, and Their Discussion // International Journal of Astronomy and Astrophysics, 2012, 2, 129-155. doi:10.4236/ijaa.2012.23018. http://www.scirp.org/journal/PaperDownload.aspx?paperID=23224.

30. Newcomb S. The elements of the four inner planets and the fundamental constants of astronomy. Washington: Government printing office. 1895. – 202 p.

31. Simon J.L., Bretagnon P., Chapront J. et. al. Numerical Expression for Precession Formulae and Mean Elements for the Moon and the Planets // Astron. and Astrophys, 1994, vol. 282. P. 663-683.

32. Мельников В.П., Смульский И.И. Астрономическая теория ледниковых периодов: Новые приближения. Решенные и нерешенные проблемы. – Новосибирск: Академическое изд-во "Гео", 2009. – 98 с. Книга на двух языках. С обратной стороны: Melnikov V.P., Smulsky J.J. Astronomical theory of ice ages: New approximations. Solutions and challenges. – Novosibirsk: Academic Publishing House "GEO", 2009. – 84 p. http://www.ikz.ru/~smulski/Papers/AsThAnR.pdf

33. Смульский И.И. Расчет взаимодействий в Солнечной системе за 50 млн. лет для изучения эволюции климата / "Большая Медведица". Журнал проблем защиты Земли. – Межрегиональный общественный фонд им. Ломоносова. Новосибирск, Россия, 2005 г., No. 1, с. 44-56. http://www.ikz.ru/~smulski/Papers/RasVSS2c.pdf.

34. Гребеников Е.А., Смульский И.И. Эволюция орбиты Марса на интервале времени в сто миллионов лет / Сообщения по прикладной математике. Российская Академия Наук: ВЦ им. А.А. Дородницына. М.: ВЦ РАН А.А. Дородницына. – 2007. – 63 с. http://www.ikz.ru/~smulski/Papers/EvMa100m4t2.pdf.

35. Morbidelli A. Modern Interactions of Solar System Dynamics // Annul Rev. Earth planet. Sci. 2002, 30. P. 89-112.

36. Laskar J., Correia A. C. M., Gastineau M., Joutel F., Levrard B. and Robutel P. Long term evolution and chaotic diffusion of the insolation quantities of Mars // Astron. and Astrophys, 2004, vol. 428. P. 261-285.

37. Laskar J. Large-scale chaos in the solar system // Astron. and Astrophys, 1994: L9 - L12.

38. Смульский И.И., Сеченов К.Е. Уравнения вращательного движения Земли и их решения при воздействии Солнца и планет / Институт криосферы Земли СО РАН. – Тюмень, 2007. – 35 с. –Деп. в ВИНИТИ 02.05.07 г. № 492-В2007. http://www.ikz.ru/~smulski/Papers/UVrVzSPc.pdf.

39. Smulsky J.J. The Influence of the Planets, Sun and Moon on the Evolution of the Earth's Axis // International Journal of Astronomy and Astrophysics, 2011, 1, 117-134. doi:10.4236/ijaa.2011.13017. http://www.SciRP.org/journal/ijaa.

40. Хайрер Э., Нёрсетт С.П., Ваннер Г. Решение обыкновенных дифференциальных уравнений. – М.: Мир, 1990. – 513 с.

41. Крутько П.Д., Максимов А.И., Скворцов Л.М. Алгоритмы и программы проектирования автоматических систем - М.: Радио и связь, 1988. – 298 с.

42. Архипов С.А. Хронология геологических событий позднего плейстоцена Западной Сибири. Геология и геофизика, 1997, т. 38, № 12. – С. 1863-1884.

43. Hargreaeves R. Distribution of Solar Radiation on the surface of the earth and its dependence on astronomical elements. Trans. Cambr. Phil. Soc. 16, 1896.

44. Комаров И.А., Исаев В.С. Криология Марса и других планет Солнечной системы. – М.: Научный мир, 2010. – 232 с.

       Печатается в соответствии с решением Ученого совета Института криосферы Земли Сибирского отделения Российской академии наук от 12 сентября 2014 г., протокол № 6.

       Приложение. Переписка с журналом "Вестник РАН"

Вестник РАН ______________________________ 119049, Москва, Крымский вал, Мароновский пер., 26 ______________________________ Тел.: 8 (499) 238-21-44 E-mail: Vestnik@Naukaran.ru

Уважаемая редакция!

       В Ваш журнал отправляю почтой свою статью: Смульский И.И. "Основные положения и новые результаты астрономической теории изменения климата".

       В разделе "Обозрение" была опубликована статья Большакова В.А. и Капицы А.П. Уроки развития орбитальной теории палеоклимата // Вестник Российской Академии Наук, 2011, т. 81, № 7, с. 603-612, а в разделе "Дискуссия" – статья Смульского И.И. Анализ уроков развития астрономической теории палеоклимата // Вестник Российской Академии Наук, 2013. Т. 83. № 1. С. 31-39., а также ответ Большакова В.А. на нее.

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

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

           Приложение

           1. Статья - на 26 стр. в 2 экз.

           Дополнительно на электронный адрес журнала в архивированном файле OsPoATLP2.zip отправляю следующие файлы:

           1) OsPoATLP2.doc – статья;

           2) OsPoATLPFg1b.gif, … 14b.gif – 14 файлов черно-белых рисунков.

           Автор     16.05.2014 г.                                                       И.И. Смульский

-------------------------------------------------------

От: "Joseph J. Smulsky" <jsmulsky@mail.ru>

Кому: "Журнал \"Вестник РАН\"" <vestnik@naukaran.ru>

Тема: Re: Из "Вестника РАН"

Дата: 1 августа 2014 г. 10:25

       В отношении предложения рецензента: Мне думается, автору следует предложить статью в журнал, где его работу могла бы оценить узкопрофессиональная аудитория" можно сказать следующее. Оно бесполезно для автора и не оправдывает решение рецензента.

С уважением               01.08.2014 г.                                                       И.И. Смульский (Иосиф Иосифович)

Дополнительная информация.

1. Скорректированную книгу (Смульский И.И. Электромагнитное и гравитационное воздействия (нерелятивистские трактаты). Новосибирск: "Наука", 1994. 225с.) Можно скопировать в djvu-формате здесь: http://www.ikz.ru/~smulski/ElGrVz.djvu и в pdf-формате здесь: http://www.ikz.ru/~smulski/ElGrVz2.pdf.

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

3. Курс лекций (Смульский И.И. Динамика (Конспект лекций по теоретической механике для строительных специальностей). Вторая редакция. – Тюмень: ТюмГАСА. – 2004 г. 88 с. 200 экз.) Можно скопировать в pdf-формате здесь: http://www.ikz.ru/~smulski/Dinmk3A5.pdf

4. Курс лекций (Смульский И.И. Курс лекций по теории взаимодействия (основания механики). – Тюмень: ТюмГАСУ. – 2005 г. 83 с. 50 экз.) Можно скопировать в pdf-формате здесь: http://www.ikz.ru/~smulski/LEKTV4A5.pdf.

5. Книгу (Smulsky J.J. The Theory of Interaction . – Ekaterinburg: Publishing house "Cultural Information Bank". – 2004, – 304 p.) Можно скопировать в pdf-формате здесь: http://www.ikz.ru/~smulski/TVEnA5_2.pdf

6. Скорректированную книгу (Смульский И.И. Аэродинамика и процессы в вихревых камерах. – Новосибирск: ВО "Наука". – 1992. – 301 с.) Можно скопировать в djvu-формате здесь: http://www.ikz.ru/~smulski/Aerpro/aerpro.djvu. Компьютерные программы к книге на Фортране можно скопировать здесь: http://www.ikz.ru/~smulski/Aerpro/PrgrmAVC.rar.

7. Информация по книге: Смульский И.И. Электродинамика движущихся тел. Определение сил и расчет движений. Saarbrucken, Germany: "Palmarium Academic Publishing", 2014. 324 с. ISBN 978-3-659-98421-1. – находится здесь http://www.ikz.ru/~smulski/Papers/InfElMvB.pdf.

8. Изменение параметров орбит за 100 млн. лет можно скопировать по адресу: http://www.ikz.ru/~smulski/Data/OrbtData/ ( см. файл ReadMeOR.pdf).

9. Система Galactica для решения задач взаимодействия n-тел по закону тяготения Ньютона или закону Кулона доступна по адресу: http://www.ikz.ru/~smulski/GalactcW/. Описание системы Galactica находится в файлах GalDiscrp.pdf и GalDiscrE.pdf на русском и английском языках, соответственно. Оно позволяет даже начинающему исследователю ставить и решать задачи с помощью программы Galactica.

10. Готовые к использованию другие Системы имеются по адресу: http://www.smul1.newmail.ru/Systems.html.

------------------------------------------------------------------------

От: ""Журнал \"Вестник РАН\""" <vestnik@naukaran.ru>