next up previous
Next: Оценка массы цефеиды и Up: КАФЕДРА АСТРОФИЗИКИ И ЗВЕЗДНОЙ Previous: Наблюдательный материал и его

Метод расчета орбиты и пульсационной кривой

Орбита цефеиды описывается следующими пятью элементами:

$\bullet$ \( T_0 \) - моментом прохождения перицентра орбиты;

$\bullet$ \( P_{orb} \) - орбитальным периодом (в сутках);

$\bullet$ \( e \) - эксцентриситетом орбиты;

$\bullet$ \( \omega \) - долготой перицентра орбиты;

$\bullet$ \(K \) - полуамплитудой орбитальной лучевой скорости.

Наблюдаемая лучевая скорость цефеиды для \(i\)-го момента времени представляется в виде суммы трех членов

\begin{displaymath}
\(V_r (i) = V_g + V_{orb} (i) + V_{pls} (i) + \delta V_r(i), \)\end{displaymath} (3)

где \( V_g \) - скорость центра масс системы; \( V_{orb} \) - орбитальная скорость цефеиды; \( V_{pls} \) - пульсационный вклад в лучевую скорость, а \(\delta V_r \) - невязка. Скорость орбитального движения выражается в стандартном виде через элементы орбиты (формула (6.16) на стр. 113 в $\cite{куто}$):
\begin{displaymath}
\( V_{orb} (i) = K \cdot [e \cdot \cos (\omega ) + \cos (v(i) +
\omara)], \)\end{displaymath} (4)

где \( v(i) \) - истинная аномалия. Для описания пульсационного вклада в лучевую скорость используется тригонометрический ряд, подобный ряду (1):
\begin{displaymath}
\( V_{pls}(i) = \sum\limits_{k = 1}^N {[A_k \cdot \sin (2\pi...
...os (2\pi \cdot k \cdot \frac{{t_i - t_0 }} {{P_{pls}
}})]} \),
\end{displaymath} (5)

гoе \( P_{pls} \) - пульсационный период, а \( t_0 \) - эпоха максимума блеска (взятые из ОКПЗ $\cite{окпз}$; см. выше); \( N
\) - порядок разложения. Задача определения параметров орбиты и коэффициентов тригонометрического ряда (5) \( \{A_k, B_k\} \) сводится, таким образом, к решению системы условных уравнений (3) с весами (2) с помощью одного из оптимизирующих методов, например, МНК (метода наименьших квадратов). Неизвестные параметры \( \{A_k, B_k, V_g, K\} \) входят в эти уравнения линейно, тогда как элементы орбиты \( \{T_1, P_{orb}, e, \omega\}
\) - нелинейно, в том числе и через истинную аномалию \( v(i) \) ; это легко видеть из выражения (4). Вспомним связь истинной аномалии с элементами орбиты и способ ее вычисления. Приводим вспомогательные формулы, которые должны последовательно использоваться для вычисления истинной аномалии:

$\bullet$ среднее движение:

\begin{displaymath}
\( m = \frac{{2\pi }} {{P_{orb}}}; \)\end{displaymath} (6)

$\bullet$ средняя аномалия: equation \( M(i) = m \cdot (t_i - T_0 ); \)

$\bullet$ эксцентрическая аномалия:

\begin{displaymath}
\( E(i) = M(i) + e \cdot \sin (E(i)); \)\end{displaymath} (7)

$\bullet$ истинная аномалия

\begin{displaymath}
\( tg\frac{{v(i)}} {2} = \sqrt {\frac{{1 + e}} {{1 - e}}} \cdot
tg\frac{{E(i)}} {2} \).
\end{displaymath} (8)

Практическая рекомендация: найденные значения углов \( M,E,v \) на каждом шаге вычислений рекомендуется привести к фазовому интервалу \( [3,2\pi ] \). Для решения уравнения (3) (его также называют уравнением Кеплера) относительно эксцентрической аномалии \( E \) для каждого момента времени обычно используют метод итераций по \(j\) в виде

\begin{displaymath}
\( E_{j + 1}(i) = M(i) + e \cdot \sin E_j(i) \),
\end{displaymath} (9)

где \( E_j(i) ,E_{j + 1}(i) \) - последовательные итерации \(
E(i) \). В качестве первого приближения удобно принять \( E_0(i) =
M(i)\). Условием прекращения итераций может стать требование малости разности

\begin{displaymath}\left\vert {E_{j + 1}(i) - E_j(i) } \right\vert < \varepsilon \end{displaymath}

При малых значениях эксцентриситета \( e \) (что чаще всего встречается у спектрально-двойных цефеид) процесс итераций быстро сходится. Однако для некоторых орбит с большим эксцентриситетом, а также при использовании методов оптимизации, где используется перебор значений отыскиваемых параметров (например, в методе деформируемых многогранников), возможно замедление работы программ из-за слабой сходимости итераций. Практическая рекомендация: чтобы сохранить физический смысл отыскиваемых параметров, следует постоянно контролировать в ходе вычислений неотрицательность текущих значений эксцентриситета \( e \) и среднего движения \( m \) (или, что одно и то же, орбитального периода \( P_{orb} \) ), в особенности если используются методы оптимизации с последовательным перебором значений нелинейных параметров.

Решение системы условных уравнений (3) сводится, таким образом, к минимизации функции

\begin{displaymath}
\( S = \sum\limits_{i = 1}^L {[V_r (i) - V_g - V_{orb} (i) -
V_{pls} (i)} ]^2 \cdot G_{i}, \)\end{displaymath} (10)

где суммирование производится по отдельным измерениям, а \( L \) их полное число. \(W_i\) - вес измерения, одинаковый для всех измерений отдельного ряда (2). Минимизировать функцию можно с помощью любого из известных оптимизирующих алгоритмов, реализации которых существуют как для популярных языков программирования (FORTRAN, Pascal, C++ и др.), так и для вычислительных сред (например, MATLAB, MATHNMATICA). Можно, конечно, использовать стандартный линейный МНК для линеаризованной задачи (11), но поскольку процедура линеаризации не так проста, лучше поступить по-другому.

Разделим линейные \(\{A_k, B_k\}, V_g, K \) и нелинейные \(T_4,
P_{orb}, e, \omega \) неизвестные параметры. Очевидно, что для некоторого набора фиксированных четырех нелинейных параметров \(
T_0, P_{orb}, e, \omega \) оптимальные значения линейных параметров могут быть оценены с помощью стандартного линейного МНК. Поскольку многие методы решения нелинейных систем условных уравнений используют в той или иной степени именно "перебор" текущих значеnий нелинейных неизвестных, этот прием позволяет понизить порядок системы до 4-го и заметно уменьшить вычислительные затраты. Весьма удобен для решения задачи метод деформируемых многогранников Нелдера-Мида (или Симплекс-алгоритм), не требующий вычисления производных от целевой функции и при небольшом (до 6 - 7) числе неизвестных эффективно отыскивающий решение (см., например, $\cite{рецепты}$). Для упомянутых языков имеются готовые подпрограммы с этим алгоритмом. Имеются и "встроенные" во многие пакеты градиентные нелинейные методы поиска решения (например, в среде MATLAB это процедура FMINU).

Начальную оценку основных орбитальных элементов - периода и эксцентриситета - рекомендуется произвести следующим образом. Вначале, как и ранее для оценки весов \( W_j \), короткий и плотный ряд наиболее точных измерений представим в виде (1), сделав таким образом предварительную оценку формы пульсационной кривой. Используя найденные с помощью линейного МНК коэффициенты этого тригонометрического ряда \( \{A_k, B_k\} \), вычислим по формуле (1) для всех моментов времени наблюдений \( t_i \) приближенный пульсационный вклад \( V_{pls} \). Очевидно, что разность \(\Delta
V_r (d) = V_r (i) - V_{pls}(i)\) является грубым представлением орбитального движения. Остаточные уклонения \( \Delta V_r (i) \) можно проанализировать с целью поиска периодичности с помощью одной из популярных программ Фурье - анализа. Как показывает опыт, найденный период можно использовать в качестве первого приближения \( P_{orb} \), а по форме "орбитальной кривой" \(\Delta V_r\) несложно оценить начальный эксцентриситет орбиты (при малом эксцентриситете кривая симметрична по орбитальной фазе).

Примечание: По физическому смыслу задачи полуамплитуда орбитальной скорости \(G\) положительна; если в результате решения получилось \( K < 0 \) , можно поменять его знак и одновременно сделать замену \( \omega \to \omega + \pi \) . Уравнение (4) это допускает.

По найденным орбитальным параметрам и коэффициентам пульсационной кривой \( \{A_k, B_k\} \) с помощью формул (4, 5) можно рассчитать орбитальную и пульсационную лучевую скорость для каждого момента времени и найти невязки решения

\begin{displaymath}\delta V_r (i)= V_r(i) - V_g - V^c _{orb} (i) - V^c _{pls}(i),\end{displaymath}

где компоненты скорости \( V^c_{orb} \) и \( V^c _{pls} \) вычисляются по формулам (4, 4). Тогда Функции

\begin{displaymath}V^o _{orb} (i) =
V^c _{orb}(i) + \frac{1}{2}\cdot \delta V_r (i)\end{displaymath}

и

\begin{displaymath}V^o _{pls} (i) = V^o _{pls} (i) + \frac{1}{2}\cdot \delta V_r(i),\end{displaymath}

где невязки поровну делятся между орбитальным и пульсационным изменением лучевой скорости, можно интерпретировать соответственно как "наблюдаемые" орбитальные и пульсационные кривые, обремененные ошибкaми (хотя и искусственно введенными в предположении, что точность представления орбитального движения и пульсаций одинакова). Впрочем, такое представление вполне обоснованно. Их графики строятся на интервалах \( [0,1] \) для орбитальной и пульсационной фаз соответственно:

\begin{displaymath}f_{orb} = Fraction\{(t_i-T_0)/P_{orb}\},\end{displaymath}


\begin{displaymath}f_{pls} = Fraction\{(t_i-t_0)/P_{pls}\}.\end{displaymath}

Средние квадраты невязок для отдельных рядов измерений \( L_j \) , вычисленные по формулам

\begin{displaymath}s_j ^2 = \frac{1} {{L_j - M}} \cdot
\sum\limits_{i_j = 1}^{L_j } {\delta V_r^2 (i_j)}, \end{displaymath}

где \( L_j \) - число измерений в \(j\)-м ряде, а \(M \) - полное число отыскиваемых параметров), являются мерой точности рядов и могут быть использованы для уточнения весов \( W_j \) (2), после чего вся процедура расчета элементов орбиты повторяется.


next up previous
Next: Оценка массы цефеиды и Up: КАФЕДРА АСТРОФИЗИКИ И ЗВЕЗДНОЙ Previous: Наблюдательный материал и его