Документ взят из кэша поисковой машины. Адрес оригинального документа : http://comet.sai.msu.ru/~dmbiz/prac/next/bincep/node3.html
Дата изменения: Thu Mar 15 16:42:25 2001
Дата индексирования: Tue Oct 2 04:19:57 2012
Кодировка: koi8-r
Поисковые слова: carl sagan

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

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

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

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

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

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

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

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

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

где - скорость центра масс системы; - орбитальная скорость цефеиды; - пульсационный вклад в лучевую скорость, а - невязка. Скорость орбитального движения выражается в стандартном виде через элементы орбиты (формула (6.16) на стр. 113 в $\cite{куто}$):

$$\( V_{orb} (i) = K \cdot [e \cdot \cos (\omega ) + \cos (v(i) + \omara)], \)$$ (4)

где - истинная аномалия. Для описания пульсационного вклада в лучевую скорость используется тригонометрический ряд, подобный ряду (1):

$$\( 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} }})]} \),$$ (5)

гoе - пульсационный период, а - эпоха максимума блеска (взятые из ОКПЗ $\cite{окпз}$; см. выше); - порядок разложения. Задача определения параметров орбиты и коэффициентов тригонометрического ряда (5) сводится, таким образом, к решению системы условных уравнений (3) с весами (2) с помощью одного из оптимизирующих методов, например, МНК (метода наименьших квадратов). Неизвестные параметры входят в эти уравнения линейно, тогда как элементы орбиты - нелинейно, в том числе и через истинную аномалию ; это легко видеть из выражения (4). Вспомним связь истинной аномалии с элементами орбиты и способ ее вычисления. Приводим вспомогательные формулы, которые должны последовательно использоваться для вычисления истинной аномалии:

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

$$\( m = \frac{{2\pi }} {{P_{orb}}}; \)$$ (6)

$\bullet$ средняя аномалия: equation

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

$$\( E(i) = M(i) + e \cdot \sin (E(i)); \)$$ (7)

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

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

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

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

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

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

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

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

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

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

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

Начальную оценку основных орбитальных элементов - периода и эксцентриситета - рекомендуется произвести следующим образом. Вначале, как и ранее для оценки весов , короткий и плотный ряд наиболее точных измерений представим в виде (1), сделав таким образом предварительную оценку формы пульсационной кривой. Используя найденные с помощью линейного МНК коэффициенты этого тригонометрического ряда , вычислим по формуле (1) для всех моментов времени наблюдений приближенный пульсационный вклад . Очевидно, что разность является грубым представлением орбитального движения. Остаточные уклонения можно проанализировать с целью поиска периодичности с помощью одной из популярных программ Фурье - анализа. Как показывает опыт, найденный период можно использовать в качестве первого приближения , а по форме "орбитальной кривой" несложно оценить начальный эксцентриситет орбиты (при малом эксцентриситете кривая симметрична по орбитальной фазе).

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

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

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

где компоненты скорости и вычисляются по формулам (4, 4). Тогда Функции

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

и

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

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

$$f_{orb} = Fraction\{(t_i-T_0)/P_{orb}\},$$

$$f_{pls} = Fraction\{(t_i-t_0)/P_{pls}\}.$$

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

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

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