КАФЕДРА АСТРОФИЗИКИ И ЗВЕЗДНОЙ АСТРОНОМИИ

КАФЕДРА ЭКСПЕРИМЕНТАЛЬНОЙ АСТРОНОМИИ

А.С.РАСТОРГУЕВ

ИЗУЧЕНИЕ ВРАЩЕНИЯ ГАЛАКТИКИ ПО ЛУЧЕВЫМ СКОРОСТЯМ И СОБСТВЕННЫМ ДВИЖЕНИЯМ ЗВЕЗД

  1. ПОСТАНОВКА ЗАДАЧИ.

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

Изучение кривой вращения галактического диска непосредственно связано со следующими актуальными направлениями исследований кинематики и динамики Галактики и астрофизики вообще:

  1. Как известно, дифференциальное вращение является основным крупномасштабным движением в диске Галактики. Как показывают исследования, уклонения от него сравнительно невелики (не более 10% от величины круговой скорости в окрестности Солнца). Обычно их объясняют возмущениями со стороны спиральных рукавов, и правильное ⌠снятие■ вращения облегчает задачу восстановления картины спирального узора по кинематическим данным и анализ особенностей поля остаточных скоростей.
  2. Изучение дифференциального вращения различных подсистем диска, различающихся возрастом и пространственным распределением, дает информацию об их происхождении и динамической эволюции.
  3. В радиоастрономии кривая вращения часто является единственным средством оценки расстояния до объекта с известной лучевой скоростью, поскольку других прямых методов локализации источников радиоизлучения (космические мазеры, молекулярные облака и др.), как правило, нет.
  4. Наблюдаемое поле скоростей весьма чувствительно к принимаемой величине расстояния до центра Галактики. Оптимизация кривой вращения позволяет оценить и этот важнейший параметр. Сама же кривая вращения плоских подсистем обычно используется как основной критерий адекватности моделей распределения массы в Галактике.

Хорошо известно, что различные подсистемы диска Галактики вращаются с разной скоростью, зависящей от дисперсии остаточных скоростей [1]. Это ⌠отставание■ центроидов от локального стандарта покоя находит полное объяснение в рамках гидродинамического описания осесимметричной галактики [2]. Наибольшей скоростью вращения обладают самые кинематически ⌠спокойные■ подсистемы с наименьшей дисперсией скоростей - газ (нейтральный и молекулярный водород), а из звездных подсистем - наиболее молодые звезды и звездные ассоциации, а также рассеянные скопления с возрастом менее лет. Именно к таким объектам относятся и классические цефеиды. Эти пульсирующие переменные звезды играют в современной астрономии исключительно важную роль благодаря наличию у них четкой связи периода пульсаций со светимостью (зависимости ⌠период - светимость■). По данным многоцветной BVRIJHK-фотометрии 9 цефеид - членов рассеянных скоплений нашей Галактики эта зависимость в лучах B и V имеет вид [3]:

,

, (1)

где - период фундаментального тона пульсаций, а нижний индекс I обозначает абсолютную величину, соответствующую потоку излучения, усредненному по периоду пульсаций. Из (1) можно вывести зависимость ⌠период - цвет■

, (2)

позволяющую учесть поглощение света и определить достаточно надежные в случайном смысле расстояния до цефеид.

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

  1. КИНЕМАТИЧЕСКАЯ МОДЕЛЬ И ОСНОВНЫЕ ОБОЗНАЧЕНИЯ.
  2. Картина движений в диске Галактике усложнена влиянием возмущений от спирального узора и другими потоковыми движениями. Однако дифференциальное вращение выделяется очень легко (см. рис. 1 и 2, на которых показана зависимость лучевой скорости и собственного движения по галактической долготе для классических цефеид по последним данным. На них легко выделяется хорошо известная ⌠двойная волна■ по долготе.) Поэтому мы рассмотрим простейшую модель чисто кругового движения центроидов. Кроме того, будем считать, что вектор вращательной скорости параллелен плоскости симметрии Галактики, а угловая скорость вращения не зависит от z-координаты. В рамках принятой модели влияние дифференциального вращения Галактики на лучевые скорости и собственные движения звезд описывается формулами Боттлингера [1]. Введем некоторые обозначения. Пусть - лучевая скорость, и - компоненты собственного движения по прямому восхождению и склонению (перед началом вычислений следует убедиться, что Вы используете локальные собственные движения ), и - компоненты собственного движения по галактической долготе и широте соответственно; и - угловая скорость исследуемого центроида на расстоянии и на расстоянии Солнца соответственно; - гелиоцентрическое расстояние объекта; - компоненты скорости Солнца относительно апекса в галактической прямоугольной системе координат (ось x которой направлена в центр Галактики, ось y - в направлении галактического вращения, z - к северному полюсу Галактики). Будем выражать собственные движения в ² /год, линейные (в том числе лучевые) скорости в км/с, расстояния в кпк, а угловые скорости в единицах км/с/кпк. Для удобства введем коэффициент перевода расстояний и собственных движений в линейные скорости, входящий в формулу , где - компонент тангенциальной скорости.

    2.1. Учет движения Солнца.

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

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

    2.2. Преобразование координат.

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

    ;

    где тот же единичный вектор в галактических координатах будет равен ; а - галактическая широта, - галактическая долгота (необходимо принять во внимание, что этот угол с учетом знаков аргументов может принимать значения от 0 до 360° ).

    Галактические компоненты собственного движения можно рассчитать с использованием следующего приема: вначале зададим вспомогательные углы ; затем последовательно вычислим:

    2.3. Формулы Боттлингера.

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

    , (3)

    где величины выражают вклад движения Солнца в лучевые и тангенциальные скорости; - остаточные значения соответственно лучевой скорости и компонентов собственного движения (при решении условных уравнений метода наименьших квадратов они рассматриваются как невязки, т.е. случайные величины. Для решения МНК два последних уравнения следует разделить на , поскольку случайными величинами являются именно остаточные уклонения собственных движений. Еще раз напомним, что если мы ставим задачу вычисления , в приведенных выше формулах следует положить ; если же мы задаем движение Солнца к некоторому апексу , то также будет одним из неизвестных параметров задачи.

    1. Метод раздельного решения уравнений Боттлингера.

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

, (4)

ограничиваясь, как правило, вторым порядком. Такое разложение дает хорошие результаты даже для расстояний порядка 5 √ 6 кпк от Солнца. Рассматривая треугольник ⌠Солнце - центр Галактики √ центроид S■ и применив к расстояниям теорему косинусов, легко вычислим по известным гелиоцентрическим расстояниям объектов и их галактическим координатам по формуле

(5)

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

Итак, неизвестными величинами, подлежащими определению, являются либо , либо . Они входят в условные уравнения линейно. Следовательно, для решения задачи могут быть использованы любые модификации линейного МНК; в пособиях приводятся также способы оценки ошибок коэффициентов регрессии [9].

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

  1. Для решения используются только лучевые скорости. Определяются значения следующих параметров: или и производные . Как правило, решение для этих параметров по собственным движениям значительно менее точное; особенно в тех случаях, когда собственные движения, как в большинстве каталогов, обременены не только случайными, но и систематическими ошибками. Однако следует отметить, что вертикальный компонент скорости по лучевым скоростям определяется ненадежно (т.к. основное движение в диске Галактики происходит в ее плоскости). Лучше всего вертикальная скорость определяется из уравнения для (3).
  2. Из уравнения для в принципе можно получить все параметры (что самое существенное, угловую скорость ), но их ошибки будут, как правило, значительно больше, чем для тех же параметров, определяемых по лучевым скоростям. Поэтому предлагается вначале получить решение для лучевых скоростей, а затем подставить найденные значения (или ) и в уравнения МНК для собственных движений по долготеи определить остальные параметры, т.е. и . И хотя в идеальном случае найденные по лучевым скоростям и по собственным движениям параметры должны быть одинаковыми, такой прием обосновывается тем, что из-за систематических ошибок в собственных движениях вторая производная от угловой скорости определяется крайне ненадежно. Сравнение значений первой производной от угловой скорости , полученных раздельно по лучевым скоростям и собственным движениям, служит одним из способов выявления систематических ошибок шкалы расстояний. Поскольку поле лучевых скоростей намного более чувствительно к используемой шкале расстояний, чем собственные движения, при удлинении шкалы расстояний величина , определяемая по лучевым скоростям, растет значительно быстрее производной, найденной по собственным движениям. Строго говоря, это один из вариантов метода статистических параллаксов [5].
  3. Уравнение для собственного движения по галактической широте лучше всего использовать только для вычисления одной неизвестной величины - вертикальной скорости Солнца , после подстановки в него всех остальных найденных параметров. Среднеквадратичное остаточное уклонение собственных движений
является хорошей независимой оценкой точности используемых собственных движений.

В качестве весов условных уравнений МНК (см., например, [9]) для системы (3) можно использовать обратные значения средних квадратов ошибок измерений .

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

  1. ОБЩЕЕ РЕШЕНИЕ УРАВНЕНИЙ БОТТЛИНГЕРА ПО ЛУЧЕВЫМ СКОРОСТЯМ И СОБСТВЕННЫМ ДВИЖЕНИЯМ.
  2. Этот метод применяется к выборке объектов, для которых одновременно известны как лучевые скорости, так и собственные движения (т.е. могут быть вычислены их пространственные скорости). Предполагается, что остаточные пространственные скорости распределены по трехосному эллипсоидальному закону Шварцшильда [1]. Для решения используется метод максимального правдоподобия, являющийся обобщением МНК.

    3.1. Обозначения.

    Предположим для простоты, что: а) форма и размеры эллипсоида не зависят от координат исследуемого центроида (в том числе и от z-координаты); б) малая ось эллипсоида параллельна оси вращения Галактики и в) большая ось эллипсоида направлена к оси вращения Галактики. Эти предположения с хорошей точностью для многих подсистем выполняются в реальной Галактике. Треугольник ⌠Солнце √ центроид S √ центр Галактики (GC)■ в проекции на галактическую плоскость схематически изображен на следующем рисунке:

    Здесь - расстояние до центроида S в проекции; - расстояние от Солнца до центра Галактики и - расстояние от оси вращения до исследуемого центроида. Оси направлены соответственно на центр Галактики и параллельно вектору линейной скорости ее вращения; вместе с осью они образуют правую тройку; - компоненты скорости в галактической прямоугольной системе координат; - лучевая и долготная тангенциальная компоненты скорости; - галактическая долгота. Вспомогательный угол , определяющий ориентацию эллипсоида относительно луча зрения, вычисляется по формуле

    , (6)

    где - гелиоцентрическое расстояние центроида S. Расстояние центроида S от оси Галактики вычисляется, как и в разделе 2.4, по формуле (5).

    3.2. Локальная система координат и тензор ковариации.

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

    . (7)

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

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

    , (8)

    так что

    .

    Пусть также форма и размеры эллипсоида скоростей задаются тензором дисперсий

    (9)

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

    , (10)

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

    , (11)

    где тензор ошибок можно записать в виде

    , (12)

    где - ошибки лучевой скорости и собственных движений. При этом мы пренебрегаем ошибками расстояния.

    Наблюдаемая скорость отражает также движение Солнца относительно исследуемой выборки объектов и влияние дифференциального вращения диска. Пусть вектор скорости местного центроида изучаемой выборки объектов относительно Солнца равен

    , (13)

    причем его вклад в локальную скорость, очевидно, равен

    ,

    где, по аналогии с (8) роль угла играет галактическая долгота , и матрица поворота для местного центроида равна

    . (14)

    Вклад дифференциального вращения Галактики описывается формулами Боттлингера (3), которые могут быть переписаны в векторном виде

    , (15)

    где, как и ранее, следует положить , если движение Солнца не задается заранее, а вычисляется.

    3.3. Распределение остаточных скоростей и функция правдоподобия.

    С учетом всего сказанного выше остаточная скорость объекта относительно центроида S с координатами в локальной системе координат равна

    , (16)

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

    , (17)

    где и - соответственно определитель и обратная матрица наблюдаемого тензора ковариации (11), а знак ⌠T■ обозначает операцию транспонирования.

    Поскольку объекты распределены в пространстве скоростей независимо друг от друга, их N-частичная (полная) функция распределения равна произведению функций (17) для всех объектов выборки:

    , (18)

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

    (19)

    и сводят задачу к поиску минимума функции правдоподобия с использованием одного из алгоритмов оптимизации. Легко видеть, что

    , (20)

    где индекс относится к текущему объекту выборки.

    В результате процедуры минимизации функции правдоподобия мы определяем , либо . Если шкала расстояний объектов надежна, можно оценить и расстояние до центра Галактики .

    3.4. Вычисление ошибок.

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

    . (21)

    Ясно, что

  3. НАБЛЮДАТЕЛЬНЫЙ МАТЕРИАЛ И ПРОВЕДЕНИЕ ВЫЧИСЛЕНИЙ.

Для выполнения работы предоставляется каталог цефеид, содержащий координаты звезд, их гелиоцентрические расстояния в кпк, вычисленные по зависимостям ⌠период - светимость■ и ⌠период √ цвет■ [3], лучевые скорости и их ошибки, взятые из компилятивной сводки в работе [7], и компоненты собственных движений и с их ошибками из каталога HIPPARCOS [8].

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

4.1. Раздельное решение уравнений Боттлингера.

После составления системы условных уравнений (см. п. 2.4) следует произвести поиск объектов, сильно уклоняющихся от общего решения для кривой вращения. Проще всего это сделать по лучевым скоростям. Проведя предварительное решение для неизвестных параметров, найдите среднеквадратичное остаточное уклонение лучевых скоростей (или среднеквадратичную невязку) от этого решения. Следует удалить из первоначальной выборки те объекты, уклонения которых от общего решения по модулю превышают . При нормальном законе распределения невязок таких объектов должно быть менее 1%, а в реальных ситуациях может быть заметно больше, например, из-за неправильной классификации объектов или грубых ошибок в вычислении расстояний или определении лучевых скоростей и собственных движений. Каждый из этих объектов может привести к смещению оценок кинематических параметров, поэтому они, безусловно, не дают правильного представления модели вращения. Рекомендуется повторить эту операцию несколько раз, пока не будут выявлены все сильно уклоняющиеся объекты. Обычно бывает достаточно 2 √ 3 итераций. Эти объекты (если их ошибки лучевых скоростей не превышают ), не следует использовать и в условных уравнениях для собственных движений, т.к. большие уклонения лучевых скоростей от выбранной модели вращения могут свидетельствовать, например, об ошибочности принятых расстояний или классификации объектов.

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

Рекомендуется решать уравнения Боттлингера (3) в нескольких вариантах: А) Задавая стандартную скорость Солнца (+10, +12, +7) км/с и вычисляя , разность между угловыми скоростями исследуемого центроида и центроида ярких звезд.

Б) Вычисляя скорость Солнца и полагая .

В) Ограничивая выборку диапазоном z-координат объектов, чтобы исключить попадание в выборку объектов, далеких от плоскости Галактики (хорошим выбором моно считать .

Г) Ограничивая выборку диапазоном расстояний от Солнца. Как известно, близкие объекты (до расстояний порядка 0.5 √ 0.8 кпк) могут принадлежать местной системе, локальному звездному комплексу и показывать общее вращение, отличное от вращения Галактики. Чтобы прояснить их влияние на общее решение, в некоторых вариантах вычислений их допустимо удалять.

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

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

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

4.2. Общее решение.

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

Существует целый ряд алгоритмов решения нелинейных задач оптимизации (минимизации функции правдоподобия). Отметим среди наиболее эффективных метод Левенберга - Маркардта [10, стр. 523] и метод деформируемых многогранников Нелдера √ Мида [10, стр. 289] (или симплекс-алгоритм). Первый из них обеспечивает быструю сходимость к решению, но требует явного (аналитического) задания производных от целевой функции по неизвестным параметрам, тогда как второй значительно медленнее, но не требует вычисления производных. В хорошо известной интегрированной среде MATLAB имеются алгоритмы минимизации, где производные вычисляются численно.

Как показывает система (3) уравнений Боттлингера, поле скоростей чувствительно к принятому значению расстояния Солнца от центра Галактики. Исходя из этого можно попытаться определить и этот параметр для задачи (20). Ясно, что ⌠оптимальное■ значение расстояния тесно связано с используемой шкалой расстояний: чем длиннее шкала, тем больше вычисленное расстояние до центра Галактики. Поэтому не имеет смысла одновременно уточнять шкалу расстояний и определять расстояние до центра Галактики по кривой вращения.

ЛИТЕРАТУРА

  1. Куликовский П.Г. ⌠Звездная астрономия■. М.: ⌠Наука■. 1985.
  2. Binney J., Tremain S. ⌠Galactic Dynamics■. P.198. Princeton: Princeton University Press. 1987.
  3. Бердников Л.Н., Возякова О.В., Дамбис А.К. ⌠Зависимость период √ светимость для классических цефеид Галактики в полосах BVRIJHK.■ Письма в астрон. журн. Т.22. N.12. С.936-944. 1996.
  4. Дамбис А.К., Мельник А.М., Расторгуев А.С. ⌠Кривая вращения системы классических цефеид и расстояние Солнца от центра Галактики■. Письма в астрон. журн. Т.21. N.5. С.331-347. 1995.
  5. Маррей К.Э. ⌠Векторная астрометрия■. Киев: ⌠Наукова думка■. 1986.
  6. Hawley S.L., Jeffreys W.H., Barnes T.G., III, Wan Lai. Absolute magnitudes and kinematic properties of RR Lyrae stars. Astrophys.J. V.302. P.626-631. 1986.
  7. Mishurov Yu.N., Zenina I.A., Dambis A.K., Mel▓nik A.M., Rastorguev A.S. ⌠Is the Sun located near the corotation circle ?■ Astron. Astrophys. V.323. P.775-780. 1997.
  8. The HIPPARCOS catalogue. ESA SP-1200. 1997.
  9. Щиголев Б.М. ⌠Математическая обработка наблюдений■. М.: ⌠Наука■. 1969.
  10. Press W.H., Flannery B.P., Teukolski S.A., Vetterling W.T. ⌠Numerical Recipes. The Art of Scientific Computing■. Cambridge: Cambridge University Press. 1987.