Научная работа: Эволюция орбиты Марса на интервале времени в сто миллионов лет

2. Сравнение результатов вычислений с аппроксимацией данных наблюдений по Ньюкомбу

2.1 Приведение вычисленной долготы перигелия к неподвижной точке на орбите

При сравнении результатов интегрирования уравнений (3) с наблюдениями возникает ряд трудностей, особенно если мы увеличиваем точность сравнения. Можно, например, сравнивать координаты тел, периоды их обращения вокруг Солнца, величины вековых возмущений. Сравнение положений тел, полученных в результате интегрирования в интервале времени, равном 50-ти годам, с эфемеридами DE406/LE406, указывает на их хорошее совпадение (см. табл. 4 из Приложения 3В). Под вековыми возмущениями мы понимаем законы изменения, прежде всего, угловых параметров орбит на интервале времени порядка одного тысячелетия. Так как вековые возмущения С. Ньюкомба, представленные им в виде «законов изменения» средних элементов эллиптической орбиты, определены на основании аппроксимации известных в астрономии наблюдений, то сравнение решений уравнений (3) с последними, позволяет убедиться в том, что они, вообще говоря, пригодны и на большем интервале времени.

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

На рис. 1 показаны основные параметры орбиты Марса (i, _р и ), которые определяются в результате интегрирования уравнений (3). Так как долгота перигелия _р = DB отсчитывается от подвижного узла D, то при изменении углов, определяющих положение плоскости орбиты (то есть величин i и ), точка D будет перемещаться вдоль линии DAB и тем самым будет «вносить вклад» в величину _р. К сожалению, невозможно зафиксировать некоторую неподвижную точку на орбите, отсчет от которой позволил бы определить чистое перемещение перигелия, не искаженное колебаниями самой орбиты. Только при совершении перигелием нескольких оборотов можно определить, более или менее точно, среднюю скорость его перемещения. Чтобы точнее «зафиксировать» мгновенное перемещение перигелия, можно, например, выбрать «условно неподвижную» точку на орбите и пересчитать положение перигелия по отношению к ней. С этой целью, через восходящий узел орбиты Земли ₀ на неподвижном экваторе (см. рис. 1), проведем большой круг, перпендикулярный к орбите Марса. Обозначим точку их пересечения буквой G. Если положение плоскости орбиты планеты относительно вектора (момента количества движения всей Солнечной Системы) определять с помощью углов Эйлера (угол прецессии -? и угол нутации- ?), тогда эволюция плоскости орбиты планеты будет определяться прецессионным движением со скоростью вокруг вектора и колебаниями угла нутации ? [17, 18]. Так как вектор близок к перпендикуляру, проведенному к плоскости орбиты Марса, то, при колебаниях орбиты GDAB, точка G будет иметь пренебрежимо малые смещения по орбите, что дает право отсчитывать от нее угловое расстояние перигелия _р0. Величина _р0 состоит из двух дуг

_p0 = GD + DB = GD + _р. (7)

Найдем из прямоугольного треугольника ₀GD дугу GD, если учесть, что дуга ₀D =,, угол GD = /2, а угол D = i. Из сферической астрономии (см., например, Справочное руководство, [8]) имеем, что

sin₀G = sin₀D·sini, (8)

cosGD = cos₀D/cos₀G. (9)

Из этих равенств можно получить, что угловое расстояние перигелия относительно неподвижной точки G равно

. (10)

2.2 Приведение данных наблюдений к неподвижной системе координат

2.2.1 Параметры орбиты Марса в подвижной системе координат по данным наблюдений

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

В работах Ньюкомба [21], опубликованных в 1895-1898 годах, приведены приближенные формулы для осредненных параметров орбиты Марса (эксцентриситета, долготы перигелия, долготы узла и наклона) в подвижной эклиптической системе координат, которые он получил на основе анализа и аппроксимации, имеющихся в его распоряжении, наблюдений:

e_a = 0.09331290 + 0.000092064·T_j - 0.000000077·T_j ²; (11)

р_a = 334^_1305.53 + 6626.73·T_j + 0.4675·T_j ² - 0.0043·T_j ³; (12)

Щ_a = 48^_4711.19 + 2775.57·T_j - 0.005·T_j ² - 00192·T_j ³; (13)

i_еа = 1^_5101.20 - 2.430·T_j + 0.0454·T_j ²; (14)

где T_j -время, отсчитываемое в юлианских столетиях по 36525 «эфемеридных» суток от фундаментальной эпохи 1900, январь, 12^h UT с номером юлианского дня JD_b = 2415020.3134. Индекс a означает, что коэффициенты в формулах получены в результате обработки наблюдений.

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

Формулы (11) - (14) заимствованы из Справочного руководства [8], в которых учтены некоторых незначительных уточнений Росса. Долгота Щ_a (см. рис. 1) также отсчитывается в плоскости подвижной орбиты Земли от ее подвижного восходящего узла г. Отметим, что р_a = Щ_a + щ_a является некоторой «условной» долготой перигелия, так как величины , определяются в разных плоскостях (щ_a - угловое расстояние точки перигелия B в плоскости орбиты Марса от его восходящего узла A, а величина Щ_a вычисляется в другой плоскости.

2.2.2 Параметры орбиты Земли в неподвижной системе координат

Так как элементы орбиты Марса отнесены к подвижной орбите Земли, сначала найдем их значения в неподвижной системе координат. Известно [8], что плоскость эклиптики «поворачивается» в своей плоскости (см. рис. 2 a) относительно мгновенной оси вращения, расположенной под углом П к оси x_e,,

П = 1735703 + 3287·T_t + 0,6·T_t², (15)

с годичной скоростью

= 0,47107 - 0,00068·T_t, (16)

где T_t - время в тропических столетиях по 36524.22 «эфемеридных» суток от фундаментальной эпохи с номером юлианского дня JD_b.

Если мы фиксируем ось x_e неподвижной системы координат в эпоху с номером юлианского дня JD_S, то проекция угловой скорости ·cosП на эту ось будет означать угловую скорость поворота эклиптики в неподвижной системе координат. Так как в неподвижной системе координат положение экватора А₀А₀ (см. рис. 2б), к которому подвижная эклиптика наклонена под углом i_Ea, известно, отсюда вытекает, что скорость изменения этого угла будет равна

·cos . (17)

Обозначим угол между эклиптикой и неподвижным экватором в неподвижной системе координат в эпоху JD через = е_a(T_S). Эта величина является начальным условием при интегрировании уравнения (17). Учтем при этом, что угол наклона подвижной эклиптики к подвижному экватору Земли определяется выражением

е_a(T_j) = 23^_2708,26 - 46,845·T_j - 0,0059·T_j² + 0,00181·T_j³,

а величина T_S = (JD_S - JD_b)/36525 представляет собой промежуток времени в юлианских столетиях от фундаментальной эпохи JD_b до эпохи JD_S.

Так как равенства (15) и (16) содержат время T, то оно входит в оба сомножителя соотношения (17). Используя среднее значение угла П для среднего момента (T_S + T_j)/2, проинтегрируем уравнение (17), в результате чего получим:

i_Ea = ?_S + 100[0,47107·(T_t - T_St) - 0,00034·(T_t - T_St)²]*

*cos1735703 + 3287·0.5·(T_St + T_t) + 0,6·0.25·(T_St + T_t)², (18)

где ?_S - наклон эклиптики к экватору в эпоху T_S, а величины

T_St = T_S·k_j, k_j = 36525/36524.22.

Формула (18) определяет «закон» изменения наклона эклиптики по отношению к неподвижной, «фиксированной» на эпоху T_S, плоскости экватора плоскости экватора, согласно данным наблюдения.

В результате решения уравнений (3) мы получаем переменные значения параметров орбиты Земли и, в частности, закон изменения угла наклона i_E плоскости её орбиты к неподвижному экватору и движение восходящего узла _E по неподвижному экватору. Законы изменения наклонов i_E и i_Ea с высокой точностью совпадают, так как их относительная разность на интервале времени 1000 лет не превышает величины .

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

_о2 = aT_n² + b T_n + c при -25 T_n 25. (19)

Здесь T_n - используемое при численном интегрировании время в юлианских столетиях от начальной эпохи JD₀,

a = 2.397234912503438·10^-6; b = 5.535655301401231·10^-5; c = -2.089118512971404·10^-6.

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

a = 2.376361216684775·10^-6; b = 4.788549396485488·10^-5; c = -2.325022430798324·10^-5.

Следует отметить, что перемещение восходящего узла орбиты Земли по неподвижному экватору _о2a мы также получили, используя средние элементы орбиты Земли, приведенные к системе 2000.0 г., в теории движения планет, разработанной во Франции

Симоном Дж.Л. и др. [20]. Отметим при этом, что величины дуг _о2 и _о2a совпадают с высокой точностью в любой из моментов рассмотренного диапазона времени -25 T_n 25.

В результате прецессии земной оси восходящий узел (см. рис. 2 б) перемещается против движения Земли, по подвижной эклиптике EE', с угловой скоростью, отнесенной к одному году (см. например, формулу из [ 8 ]),

p = 50,25641 + 0,02223 T_t. (20)

Тогда в любой момент времени Т смещение точки весеннего равноденствия по подвижной эклиптике относительно момента Т_S будет определяться дугой ₂ (см. рис. 2 б), полученной в результате интегрирования равенства (20), в следующем виде:

₂ = 5025.641·(T_t - T_St) + 2.223·(T_t - T_St)²/2. (21)

Для сравнения с результатами интегрирования уравнений (3) необходимо «время» в тропических столетиях T_t в зависимостях (18) и (21) заменить формулой: T_t = k_j·[T_n + (JD₀ - JD_b)/36525]. Таким образом, наклон подвижной эклиптики к неподвижному экватору определяется выражением (18), перемещение ее восходящего узла - формулой (19), а положение подвижного экватора - равенством (21).

2.2.3 Параметры орбиты Марса в неподвижной системе координат

Долготы _а и _а восходящего узла (точки D на рис. 2 б) и перигелия (точка B) Марса, соответственно, отсчитываются от точки , поэтому на их величины влияет прецессия ₂ восходящего узла Земли. Аналогично, угол наклона плоскости орбиты Марса i_еа зависит от изменения угла наклона орбиты Земли i_Ea, поэтому целесообразно вычислить вековые возмущения орбиты Марса относительно неподвижной плоскости экватора 1950.0 г. A₀A₀' (см. рис. 2 б) как функции от следующих параметров: 1) углового расстояния восходящего узла ^a=г₀D, отсчитываемого в плоскости экватора 1950 г. от восходящего узла Земли ₀; 2) _p0^a = GB - углового положения перигелия, отсчитываемого в плоскости орбиты Марса от неподвижной точки G; 3) угла наклона i^a = DG плоскости орбиты к неподвижной плоскости экватора.