Построение кривой вращения галактики по данным длиннощелевой спектроскопии¶
Составители:¶
- к.ф.-м.н. Катков Иван Юрьевич
- профессор д.ф.-м.н. Засов Анатолий Владимирович
Методическое пособие к задаче специального астрофизического практикума для студентов 3-4го курсов астрономического отделения ФФ МГУ.
Данная задача знакомит студентов с:
- основами первичной редукции спектральных данных
- продвинутой методикой по-пиксельной аппроксимации спектров
pPXFдля восстановления информации о движении звезд и ионизованного газа в галактике.
Опубликовано: 10 сентября 2020; Последнее обновление: 24 октября 2020
Содержание:
- 1. Общее описание задачи
- 1.1 Введение
- 1.2 Наблюдательные данные
- 1.3 Первичная редукция данных
- 1.3.1 Какими данными располагаем?
- 1.3.2 Построение супербайеса
- 1.3.3 Вычитание байеса и кадр ошибок
- 1.3.4 Учет неравномерной засветки и вариаций чувствительности матрицы
- 1.3.5 Чистка следов космических частиц
- 1.3.6 Сложение экспозиций
- 1.3.7 Калибровка шкалы длин волн
- 1.3.8 Вычитание спектра ночного неба
- 1.3.9 Заключительные замечания к первичной редукции
- 1.4 Анализ редуцированного спектра
- 2. Выполнение задачи
- 2.1 Первичная редукция
- 2.2 Анализ спектра
- 2.3 Заключительные замечания
1. Общее описание задачи¶
1.1 Введение¶
Практически во всех современных астрономических наблюдениях для регистрации сигнала от исследуемого объекта используются ПЗС приемники. Сигнал с ПЗС матрицы считывается и записывается для хранения в формате FITS - Flexible Image Transport System, который стал стандартным форматом данных в астрономии. В FITS формате можно хранить изображения, многомерные массивы данных (например спектральные кубы, получаемые в результате анализа наблюдений с панорамным спектрографом) и табличные данные. Кроме самих данных, FITS формат позволяет хранить мета информацию в заголовке файла (FITS header), где может быть записана информация о конфигурации телескопа, времени наблюдения, качестве атмосферы в момент наблюдения и т.п.
Note: Важным аспектом анализа астрономических данных является подробно и полно сформированные FITS заголовки не только "сырых" данных (т.е. тех, что записываются непосредственно в ходе наблюдений), но и в редуцированных данных. Например, отредуцированное изображения звездного поля или галактики должны иметь корректные поля FITS заголовка, определяющие координатную сетку (WCS - World Coordinate System - fields). В FITS заголовке файла со спектром должны быть правильно прописаны поля, определяющие шкалу длин волн спектра. WCS стандарт для небесных и спектральных координат описан в статьях, приведенных на этой странице. Правильно оформленный FITS файл будет корректно отображаться в стандартных астрономических программах (ds9, Aladin и др.).
Основной целью этой задачи являются знакомство с процессом обработки и анализа данных, полученных при спектральных наблюдениях галактик.
Спектры галактик содержат в себе большое количество полезной физической информации о них. Во-первых, благодаря эффекту Доплера, по спектрам можно определить среднюю скорость движения звезд и ионизованного газа в галактике. Во-вторых, сам вид спектра позволяет сказать из совокупности каких звезд состоит галактика - какой средний возраст и металличность звездного населения. И, наконец, анализируя отношения эмиссионных линий ионизованного газа, можно получить информацию о его состоянии (температура, электронная плотность) и о механизме, ответственном за его ионизацию. В настоящей задаче мы остановимся на восстановлении кинематики (скорости $v$ и дисперсии скоростей $\sigma$, точнее их проекции на луч зрения) звезд и газа для одной из галактик, отнаблюденной на 6-м Большом Телескопе Азимутальном (БТА) Специальной Астрофизической Обсерватории.
В ходе работы над задачей ожидается, что из предоставленных "сырых" данных длинощелевых спектральных наблюдений будут:
- получен редуцированный и готовый к научному анализу спектр галактики;
- определена кривая вращения галактики по звездам - т.е. зависимость скорости вращения от расстояния до центра галактики;
- получена кривая вращения ионизованного газа по эмиссионным линиям.
Эти результаты могут быть использованы для определения массы галактик, динамического состояния звездных дисков, радиальных градиентов параметров звездного населения и состояния ионизованного газа, что требуется в широком круге научных задач исследования галактик. Методы и подходы к анализу данных, предлагаемые в этой задаче, могут применяться к любым спектральным данным.
1.2 Наблюдательные данные¶
Для выполнения задачи предлагается работать со спектрами одной галактики, полученными на 6-м телескопе Специальной Астрофизической Обсерватории (САО РАН) с помощью фокального редуктора светосилы SCORPIO (или SCORPIO-2), установленных в прямом фокусе телескопа.
Ниже приводится пример полного цикла обработки и анализа одной галактики, наблюдения которой проводились на SCORPIO в длиннощелевой моде, размер щели $1''\times 6'$. Спектральное разрешение спектрографа около $\approx2.5$ Å, что в рассматриваемом диапазоне 4800-5500 Å составляет в терминах скорости
$$\frac{\delta \lambda}{2.355 \lambda} \cdot c_0 \approx 60-65~\mathrm{км/c}$$Этот диапазон включает как контрастные абсорбционные линии звездного населения, так и сильные эмиссионные линии ионизованного газа $H\beta$ 4861.325 Å, дублет [OIII] 4958.911,5006.843 Å, что дает возможность восстанавливать из кинематику звезд и газа.
1.3 Первичная редукция данных¶
В этом параграфе обсудим основные этапы первичной редукции спектральных данных.
1.3.1 Какими данными располагаем?¶
Самым первым шагом является знакомство с данными. В первую очередь нужно понять какие данные имеются в руках, какие снимались калибровки, все ли научные и калибровочные экспозиции сняты без сбоев аппаратуры. Желательно просмотреть каждый кадр, например, в ds9.
Как правило вместе с данными предоставляется наблюдательный лог-файл. Если такого нет, его можно создать на основе нужных полей из FITS заголовков всех предоставленных файлов (подробнее см. во второй части описания - Выполнение задачи ).
Фрагмент лог-файла наблюдений на 6-м телескопе 04.05.2008:
|S5630306| HZ 44 | obj|01:44:59| 60 | Spectra |VPHG2300G | 4800-5570AA| | mask |13:24:00.49| +36:05:19.1|0.0 | 23| 81| 62.2|223.9|'
|S5630307| HZ 44 |flat|01:53:15| 20 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:16:31.20| +35:18:09.7|0.0 | 16| 65| 53.5|224.8|'
|S5630308| HZ 44 |flat|01:54:09| 20 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:16:31.20| +35:18:09.7|0.0 | 16| 65| 53.8|225.0|'
|S5630401| NGC 5533 | obj|02:02:37| 30 | Image | |5470A EW=790| | V |14:16:31.20| +35:18:09.7|0.0 | 17| 68| 55.8|160.3|' '
|S5630402| NGC 5533 | obj|02:06:09| 20 | Image | |5470A EW=790| | V |14:16:31.56| +35:18:07.3|0.0 | 18| 70| 56.5|161.0|' '
|S5630403| NGC 5533 |neon|02:15:25| 60 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:16:31.56| +35:18:07.3|0.0 | 19| 73| 58.1|162.6|' '
|S5630406| NGC 5533 | obj|02:24:32| 1200 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:16:31.68| +35:18:05.8|0.0 | 21| 76| 59.3|163.8|' '
|S5630407| NGC 5533 | obj|02:45:41| 1200 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:16:31.68| +35:18:05.8|0.0 | 25| 81| 61.2|165.7|'
|S5630408| NGC 5533 | obj|03:06:37| 1200 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:16:31.68| +35:18:05.8|0.0 | 29| 85| 62.2|166.6|'
|S5630409| NGC 5533 | obj|03:30:45| 1200 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:16:31.68| +35:18:05.8|0.0 | 33| 90| 62.4|166.9|'
|S5630410| NGC 5533 |neon|03:52:31| 60 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:16:31.68| +35:18:05.8|0.0 | 37| 93| 62.2|166.7|'
|S5630411| NGC 5533 | obj|03:55:09| 20 | Image | |5470A EW=790| | V |14:16:31.68| +35:18:05.8|0.0 | 37| 94| 62.1|166.6|'
|S5630412| NGC 5533 | obj|03:57:23| 20 | Image | |5470A EW=790| | V |14:16:31.31| +35:18:08.2|0.0 | 38| 94| 62.1|166.6|' '
|S5630413| NGC 5533 | obj|03:59:03| 20 | Image | |5470A EW=790| | V |14:16:31.36| +35:18:07.9|0.0 | 38| 94| 62.1|166.5|' '
|S5630414| NGC 5533 |neon|04:01:49| 60 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:16:31.36| +35:18:07.9|0.0 | 39| 95| 62.0|166.5|' '
|S5630415| NGC 5533 | obj|04:03:32| 1200 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:16:31.36| +35:18:07.9|0.0 | 39| 95| 61.9|166.4|'
|S5630416| NGC 5533 | obj|04:24:32| 1200 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:16:31.36| +35:18:07.9|0.0 | 43| 98| 61.2|165.7|'
|S5630417| NGC 5533 |neon|04:45:55| 60 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:16:31.36| +35:18:07.9|0.0 | 47| 101| 60.3|164.7|'
|S5630418| NGC 5533 |flat|04:47:52| 20 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:16:31.36| +35:18:07.9|0.0 | 47| 101| 60.2|164.6|'
|S5630419| NGC 5533 |flat|04:48:40| 20 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:16:31.36| +35:18:07.9|0.0 | 47| 102| 60.1|164.6|'
|S5630501| HD135722 | obj|04:56:53| 10 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |15:15:52.81| +33:16:53.5|0.0 | 39| 91| 59.9|164.5|'
|S5630502| HD135722 | obj|04:58:40| 5 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |15:15:52.81| +33:16:53.5|0.0 | 39| 92| 59.9|164.5|'
|S5630503| HD135722 | obj|05:00:04| 5 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |15:15:52.81| +33:16:53.5|0.0 | 39| 92| 59.9|164.5|'
|S5630504| HD135722 |neon|05:01:08| 60 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |15:15:52.81| +33:16:53.5|0.0 | 40| 92| 59.8|164.5|'
|S5630601| HD147677 | obj|05:06:46| 5 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |16:22:26.50| +30:51:40.0|0.0 | 30| 77| 55.2|164.8|'
|S5630602| HD147677 | obj|05:08:16| 5 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |16:22:26.50| +30:51:40.0|0.0 | 30| 77| 55.3|164.9|'
|S5630603| HD147677 |neon|05:09:20| 60 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |16:22:26.50| +30:51:40.0|0.0 | 30| 77| 55.4|165.0|'
|S5630701| Flat field | obj|05:14:49| 20 | Image | |5470A EW=790| | V |14:09:39.84| +59:37:28.9|0.0 | 44| 136| 78.5|164.7|'
|S5630702| Flat field | obj|05:15:44| 15 | Image | |5470A EW=790| | V |14:09:39.84| +59:37:28.9|0.0 | 44| 136| 78.3|164.5|'
|S5630703| Flat field | obj|05:17:00| 15 | Image | |5470A EW=790| | V |14:09:39.84| +59:37:28.9|0.0 | 44| 136| 78.1|164.3|'
|S5630704| sunsky | obj|05:19:14| 60 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:09:39.84| +59:37:28.9|0.0 | 44| 136| 77.7|163.8|'
|S5630705| sunsky | obj|05:21:52| 180 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:09:39.84| +59:37:28.9|0.0 | 44| 137| 77.2|163.3|'
|S5630706| sunsky | obj|05:27:22| 120 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:09:39.84| +59:37:28.9|0.0 | 45| 137| 76.1|162.3|'
|S5630707| sunsky | obj|05:32:15| 60 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |14:09:39.84| +59:37:28.9|0.0 | 46| 137| 75.2|161.4|'
|S5630708| sunsky |neon|05:37:37| 60 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |22:33:00.60| -22:03:03.2|0.0 | 79| -44|325.8|160.8|'
|S5630801| bias 1x2 |bias|05:39:38| 0 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |22:35:04.48| -22:04:33.5|0.0 | 80| -44|325.8|160.8|'
|S5630802| bias 1x2 |bias|05:40:22| 0 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |22:35:51.73| -22:05:15.7|0.0 | 80| -44|325.8|160.8|'
|S5630803| bias 1x2 |bias|05:41:07| 0 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |22:36:37.00| -22:05:17.6|0.0 | 80| -44|325.8|160.8|'
|S5630804| bias 1x2 |bias|05:41:53| 0 | Spectra |VPHG2300G | 4800-5570AA| | slit_1.0 |22:37:22.88| -22:05:17.6|0.0 | 80| -44|325.8|160.8|'
1.3.2 Построение супербайеса¶
Кадры байеса являются калибровочными кадрами, которые снимаются с нулевой экспозицией. Несмотря на то, что на матрицу не попадает фотонов, полученный кадр имеет ненулевые отсчеты - смещение (bias) относительно нуля. Такое же смещение имеется на всех кадрах, включая научные и другие калибровочные. Чтобы учесть это смещение отсчетов, необходимо построить средний кадр байеса (супербайес) по всем имеющимся bias кадрам.
Внимание: В ходе наблюдательной ночи могут использовать разные режимы считывания ПЗС матрицы (разное бинирование и быстрое/медленное считывание, которые имеют разные характеристики, например, разный шум считывания). Для построения корректного супербайеса нужно использовать bias кадры снятые в том же режиме считывания, что и научные кадры.
ПЗС матрицы обладают темновым током, который приводит к росту отсчетов со временем без реального попадания фотонов на матрицу. Чтобы учесть этот эффект снимают так называемые dark кадры, т.е. делают экспозицию без засветки матрицы. Потом эти кадры необходимо усреднить, отмасштабировать на экспозицию текущего кадра и вычесть. Для высококачественных ПЗС матриц характерен очень низкий уровень темнового тока. Например, используемая в SCORPIO матрица EEV42-40 обладает уровнем темнового тока $\sim 0.03e^-$, поэтому калибровочные накопления dark не снимаются.
1.3.3 Вычитание байеса и кадр ошибок¶
Из всех кадров, которые будут использоваться в редукции (спектр объекта, флэт, арки), необходимо вычесть усредненный байес (супербайес). На этом же этапе рекомендуется из исходных изображений вырезать области, которые остались незасвеченными (overscan).
Note: В случаях, когда по какой-то причине bias кадры не снимались - можно оценить уровень bias по overscan области (если темновой ток и рассеянный свет пренебрежимо малы). В этой ситуации bias будет просто константой вместо попиксельного кадра, который может учитывать неравномерность смещения по матрице.
Очень важным элементом анализа данных является оценка ошибок определяемых параметров. Для корректного расчета ошибок при анализе спектра галактики, нужно помимо спектра иметь кадр ошибки. Эти кадры рассчитываются в самом начале редукции данных и пропагируются через все этапы редукции и анализа спектра.
Кадры ошибок вычисляются в предположении, что ошибки определяются пуассоновским фотонным шумом и шумом считывания матрицы. Величину шума считывания можно оценить из распределения ошибок относительно среднего уровня на отдельном кадре байеса (НЕ супербайеса!). Иногда эта величина записана в FITS заголовок, например, под полем READNOIS. Типичное значение шума считывания для используемой в SCORPIO матрицы $2e^--4e^-$ Кадр ошибки формируются в соответствии со следующей формулой:
$$
E_I = \frac{\sqrt{ I\cdot G + R^2}}{G},
$$
где $R$ - шум считывания матрицы, $G$ - коэффициент усиления матрицы (gain), $I$ - текущий кадр накопления в отсчетах (ADU). Коэффициент усиления как правило прописан в FITS заголовке - GAIN, характерные значения $\sim 0.5 - 2 e^-/ADU$.
1.3.4 Учет неравномерной засветки и вариаций чувствительности матрицы¶
Сложность оптической схемы спектрографа приводит к тому, что матрица засвечивается неравномерно. Кроме того, матрица может иметь неодинаковую чувствительность от пикселя к пикселю, а так же дефекты. Для учета обоих эффектов снимают специальные калибровочные кадры - плоское поле (flat). Для этого используется источник с плавным непрерывным спектром без явных абсорбционных или эмиссионных особенностей. Для этого хорошо подходит спектр лампы накаливания.
Щель спектрографа засвечивается равномерным светом от лампы, но в результате особенностей устройства спектрографа, матрица (вдоль щели, т.е. по Y) будет засвечена неравномерно. Чтобы это учесть, сначала кадры плоского поля усредняются и вычитаются следы космических частиц (об этом подробнее в следующем параграфе). Затем следует рассчитать относительные коэффициенты путем деления каждой строки кадра плоского поля усредненный спектр в центральной области кадра. Желательно использовать медианное усреднение, а центральную область выбрать шириной 30-50 px. После этого, все научные кадры делятся на кадр нормированных коэффициентов, что убирает как попиксельные вариации, так и неравномерную засветку. Кадры с пуасоновскими ошибками вычисляются по правилами вычислений ошибок.
1.3.5 Чистка следов космических частиц¶
Каждая экспозиция объекта длится 10-20 минут. За это время через матрицу проходит достаточно много высокоэнергетичных частиц, которые оставляют следы на матрице (космики, cosmic hits, CH). Эти следы необходимо убрать. Для этого есть несколько подходов:
- Самым простым вариантом является попарное сравнение двух кадров. Такой подход применим когда имеется более одной научной экспозиции объекта. Пиксели, в которых испорчены космиками, определяются из условия, что разница двух кадров больше наперед заданного предела. Это очень простая и эффективная методика. Но ее сложно применять, если научные экспозиции заметно отличаются друг от друга, например, в из-за растущего фона неба, либо переменной прозрачности в ходе наблюдений.
- Более продвинутый подход - это методика L.A.Cosmic (Laplacian Cosmic Ray Identification), в котором космики обнаруживаются с помощью лапласовского фильтра. Реализация этой методики есть для Python, IDL, IRAF и С++.
При использовании любого алгоритма чистки следов космических частиц следует следить за получаемой маской. При плохо подобранных параметрах ошибочно за космики могут приниматься яркие эмиссионные линии в галактике или яркий центр галактики. Рекомендуется каждую маску сохранять в файл или выводить в виде картинки и контролировать качество идентификации космиков.
1.3.6 Сложение экспозиций¶
После чистки космиков можно складывать отдельные кадры объекта и неона. Соответствующие кадры ошибок так же суммируются, но по правилам вычисления ошибок, т.е. суммируются квадратично.
Внимание: При сложении необходимо обращать внимание на относительные сдвиги между кадрами. Спектр объекта может смещаться от кадра к кадру как по X, т.е. вдоль направления дисперсии, так и по Y (вдоль щели). Эти смещения могут быть связаны с внутренними гнутиями прибора, которые могут изменяться в течении длительных наблюдений, поскольку спектрограф находится не в стационарном положении, а закреплен в прямом фокусе телескопа (в стакане). Кроме того, из-за проблем с гидированием объект может смещаться вдоль щели. Если замечены существенные сдвиги необходимо перед сложением произвести соответствующие коррекции кадров. Грубо на глаз смещения можно проверить с помощью программы ds9 в режиме блинкования. Более точно сдвиги можно определить по положению пика кросс-корреляционной функции посчитанной между реперным и текущим кадрами. Точность кросс-корреляции для определения сдвига очень высока и может достигать $\sim0.1$ px.
1.3.7 Калибровка шкалы длин волн¶
Для определения шкалы длин волн в ходе наблюдений снимаются специальные калибровочные наблюдения, которые называются арками (arc frames). В случае прибора SCORPIO используется лампа с гелий-неон-аргоновым (He-Ne-Ar) наполнением (обычно такие спектры называют просто спектрами неона). Спектр такой лампы состоит из набора узких эмиссионных линий, которые в результате особенностей оптической схемы на матрице имеют вид дуг (отсюда название arcs). Длины волн эмиссионных линий известны, поэтому отождествив наиболее яркие линии и измерив из точные положения в пиксельных координатах можно определить дисперсионное калибровочное соотношение $\lambda = f(x)$. Полиномиальная зависимость хорошо подходит для аппроксимации этой зависимости:
$$\lambda = \sum_{i=0}^{i=N} k_i x^i$$Для спектров SCORPIO c решеткой 2300G вполне достаточно использовать полином 4й-5й степени ($m=4..5$) по длине волны.
Спектр неона с отождествленными линиями, снятый с решеткой 2300G (из Руководства пользователя SCORPIO):
Линии равных длин волн сильно искривлены на кадрах неона, поэтому необходимо идентифицировать реперные линии в разных положениях по Y - $y_j$. После чего определять калибровочное соотношение $\lambda = f(x)$ для каждого положения на щели (для данной "строки" на кадре). Таким образом, на выходе этой процедуры будут определены длины волн всех пикселей изображения.
Можно использовать несколько иной подход: сначала идентифицировать положения известных линий спектра неона по всему кадру, что даст набор $\lambda_{k}$ в разных положениях на кадре ($x_k$, $y_k$). После чего следует искать дисперсионное решение в виде двумерного полинома:
$$\lambda = \sum_{i=0}^{m} \sum_{j=0}^{n} k_{ij} x^i y^j,$$где $m$ и $n$ - степени полиномов вдоль направления дисперсии (по X) и вдоль щели (по Y), соответственно.
Заключительным шагом калибровки шкалы длин волн является линеаризация спектра, т.е. процедура приведения исходного кадра к равномерной шкале длин волн по X. Это осуществляется путем интерполяции исходного спектра $F(x,y)$ к $F(\lambda, y)$, где шкала длин волн уже линейная.
Note: Удобно интерполировать исходный кадр к логарифмической шкале длин волн. Это дает преимущество для определения скоростей и дисперсии скоростей, в том числе при по-пиксельном моделировании спектра. Эффект Доплера приводит к смещению спектра на величину $\delta \lambda = \lambda_0 \cdot v/c$, который зависит от длины волны. Если спектр приведен к логарифмической шкале (т.е. логбинирован), тогда сдвиг $\delta \ln \lambda$ зависит только от скорости (Покажите это.). Иными словами, изменение спектра за счет скорости объекта приводит к линейному сдвигу логбинированного спектра. Это факт используется при моделировании спектров галактик, определении скорости и дисперсии скоростей звезд и газа.
Внимание: Чтобы спектры корректно отображались в ds9 и других приложениях нужно правильно заполнить ключевые слова FITS файла:
CDELT1,CD1_1,CRVAL1,CRPIX1и правильно указатьCTYPE1. Для линейной шкалы длин волн в воздухеCTYPE1='AWAV-LIN', для логбинированной шкалыCTYPE1='AWAV-LOG'.
Для контроля линеаризации рекомендуется линеаризовать кадры неона, чтобы убедиться что калибровочное соотношение уравнение $\lambda = f(x)$ и линеаризация выполнены правильно. После линеаризации линии неона должны быть прямыми.
1.3.8 Вычитание спектра ночного неба¶
При экспозиции объекта на щель попадает свет не только от объекта, но и от ночного неба. В итоге результирующий спектр есть сумма спектра объекта и спектра ночного неба.
В галактиках мы интересуемся не только центральными областями, где поверхностная яркость как правило в разы превышает поверхностную яркость ночного неба (если наблюдения проводятся в темную безлунную ночь), но и далекой периферией дисков, где все наоборот - ночное небо доминирует в спектре. Поэтому для изучения, необходимо вычитать вклад от спектра ночного неба.
В первом приближении мы можем предполагать, что спектр неба вдоль щели не меняется. Поэтому достаточно рассчитать спектр неба в области на щели, где нет объекта, и вычесть этот средний спектр со всей щели. Главная проблема здесь в том, что в действительности спектр неба в разных положения на щели немного отличается. Главным источником этого отличия являются вариации инструментального контура (LSF, line spread function) спектрографа. Для SCORPIO (в меньшей степени для SCOORPIO-2) характерны сильные вариации LSF вдоль щели, что проявляется как изменение формы ярких эмиссионных линий, так и, хоть менее заметно, абсорбционных линий спектра ночного континуума (отраженный и рассеянный солнечный спектр). Также этот эффект хорошо видно на спектре неона - видно как профиль сильных линий меняется вдоль щели. Поэтому простое вычитание усредненного спектра приводит к перевычету или недовычету спектра неба.
Есть продвинутые методики для учета вариаций контура (Katkov & Chilingarian 2011), но можно воспользоваться простым подходом, который позволяет нивелировать этот эффект. Для этого следует аппроксимировать спектр неба в каждом столбце изображения вне области объекта полиномом невысокой степени (2-4) и экстраполировать на положение объекта. Построенная таким образом модель ночного неба вычитается из кадра.
Спектр галактики и спектр ночного неба. Типичное положение областей для определения спектра ночного неба показано полупрозрачными полосами.
1.3.9 Заключительные замечания к первичной редукции¶
Линеаризованный суммарный спектр галактики с вычтенным спектром неба можно анализировать и извлекать различные параметры. Предлагаемая последовательность этапов первичной редукции может быть расширена, дополнена. Во-первых, при наличии спектра звезды спектрофотометрического стандарта рекомендуется для него проделать все этапы первичной редукции, сделать экстракцию одномерного спектра, сравнить с эталонным спектром и получить кривую спектральной чувствительности, которую следует учесть в спектре галактики. Этот шаг особенно важен для анализа отношений потоков далеко расположенных эмиссионных линий и для анализа свойств звездного населения. На определение кинематических параметров, спектральная чувствительность влияет слабо.
При редукции "сырых" данных следует стремиться к минимальному числу интерполяции кадров. Каждая интерполяция вносит артефакты в кадры, в частности нарушает статистику отсчетов, что приводит к систематическому смещению уровня ошибки в кадре. Пропагированные кадры ошибок через все этапы редукции будут лишь приближенно соответствовать реальному уровню ошибок в научном кадре. Оценить это смещение можно по разбросу уровня отсчетов (RMS) в области где нет линий на краю галактики, сопоставив его с рассчитанным кадром ошибки.
Note: Для лучшего понимания используемых в задаче данных и их редукции очень рекомендуется ознакомиться с замечательным циклом лекций А.В.Моисеева (САО РАН) "Многорежимный фокальный редуктор телескопа БТА" (2016/2020 гг.)
1.4 Анализ редуцированного спектра¶
1.4.1 Адаптивное бинирование спектра¶
Наша главная цель анализа - это получить скорости звезд и газа и дисперсию скоростей звезд в зависимости от расстояния от центра галактики (получить радиальные профили этих параметров). Мы можем анализировать спектр галактики в каждой строке, но, во-первых, это несколько сотен спектров, а главное - сигнал-шум в данных стремительно падает с расстоянием от центра, поскольку падает поверхностная яркость галактики. Адаптивный бининг - это простая методика и она позволяет протянуть профили параметров немного дальше. Суть в том, что спектр бинируется в бины разной длины - в центре, где SNR высокий бины могут быть размером 1-2 пикселя, а дальше размер бина (то есть области суммирования или усреднения спектра) увеличивается. Адаптивный означает, что размер бинов адаптивно увеличивается, чтобы обеспечить наперед заданный SNR в бинированном спектре.
1.4.2 По-пиксельная моделирование спектра галактики¶
Для анализа редуцированного спектра галактики мы будем использовать Penalized PiXel-Fitting (pPXF) метод (Cappellari & Emsellem (2004), Cappellari (2017)), который является одним из наиболее популярных публично доступных методик по-пиксельного моделирования галактик.
Его суть сводится к моделированию спектра галактики путем $\chi^2$ минимизации, в ходе которой определяются оптимальные параметры модели. Модельный спектр галактики представляется линейной комбинацией темплейтов, свернутых с распределением звезд по лучевым скоростям (LOSVD - Line-Of-Sight Velocity Distribution).
- LOSVD в
ppxfзадается в виде Гаусс-Эрмитового профиля: $$y=(v-v_{sys})/\sigma$$ $$L = \frac{1}{\sqrt{2\pi\sigma^2}} \exp{(-0.5 y^2)} \cdot (1 + h_3 H_3(y) + h_4 H_4(y))$$ Здесь $v$, $\sigma$ - это скорость и дисперсия скоростей на луче зрения. $H_3(y)$, $H_4(y)$ - полиномы Эрмита, а $h3$, $h4$ - коэффициенты Эрмита. Эрмитова добавка позволяет описывать несиметричные отклонения (в виде крыльев) от чисто Гауссовой формы LOSVD за счет коэффициента $h3$, тогда как $h4$ позволяет получить широкие симметричные крылья и острый пик в центре или, наоборот, профиль с резко убывающими крыльями и полкой в центре ("ящикоподобный" П-образный профиль). - В качестве темплейтов в
pPXFопционально могут использовать набор спектров звезд или набор моделей звездных населений. ВpPXFможно добавить произвольные темплейты (например описывающий спектр ночного неба), но в Python пакете уже идут модели простых звездных населений на основе звездной библиотеки MILES. Также в качестве темплейтов можно добавить набор эмиссионных линий, тогда скорости звезд и газа (по разным линиям) будут определяться одновременно. - Для наилучшего соответствия между линейной комбинацией темплетов, свернутых с параметризованной LOSVD, в модели есть мультипликативный и аддитивный континуумы, которые задаются в виде полиномов Лежандра.
- Суммарную модель можно записать так $$T^{full}=\sum_k c_k T_k(\tilde{\lambda})$$ $$M_i(\tilde{\lambda}) = m(\tilde{\lambda}) \left( \sum_q T_q^{full}(\tilde{\lambda}) \ast L(v_q, \sigma_q) \right) + a(\tilde{\lambda}),$$ здесь $T_k(\tilde{\lambda})$ индивидуальный темплейт - в нашем случае это спектры звездного населения MILES различного возраста и металличности, суммарный темплейт звездного населения $T_{star}^{full}$ сворачивается со звездной LOSVD $L(v_{star}, \sigma_{star})$, кроме того в нашей модели будет темплейт, описывающий эмиссионные линии в виде Гауссиан с $T_{gas}^{full} = T_{H\beta} + T_{[OIII]}$. $m(\tilde{\lambda})$, $a(\tilde{\lambda})$ - мультипликативный и аддитивный континуумы. Ниже мы будем использовать только мультипликативный континуум, потому что аддитивный следует использовать с особой осторожностью (только там, где действительно ожидается аддитивная добавка в наблюдаемом спектре), потому что он может систематически влиять на дисперсию скоростей звездного континуума $\tilde{\lambda}$ означает, что мы используем спектры в логарифмической шкале длин волн $\ln \lambda$. $\ast$ - означает операцию свертки. Полный $\chi^2$ будет выглядеть следующим образом: $$\chi^2 = \sum_i \frac{S_i - M_i}{\delta S_i},$$ где $S_i$ - наблюдаемый спектр галактики, $\delta S_i$ - ошибки спектра.
2. Выполнение задачи¶
Для уверенного выполнения задачи необходимы базовые навыки работы c python, а выполнение задачи удобно представить в виде jupyter notebook. Ниже приводится примерная реализация этапов редукции и анализа данных на примере спектров галактики NGC 5533, полученных с прибором SCORPIO на 6м телескопе БТА в рамках программы А.В.Засова по исследованию динамического состояния дисков галактик.
Композитное изображение галактики NGC 5533 из обзора Legacy Survey (https://www.legacysurvey.org/).
2.1 Первичная редукция¶
Далее предполагается, что у Вас уже настроено python окружение и установлены следующие пакеты:
numpyastropymatplotliblacosmicscikit-imageтребуется дляlacosmictqdm- красивый прогресс-барvorbin- адаптивный бининг (опционально)ppxf- попиксельное моделирования спектра галактики
Если какого-то из пакетов не хватает, его можно установить прямо из jupyter notebook выполнив команду !pip install lacosmic.
2.1.0 Где брать данные¶
Необходимые данные можно найти тут http://gal-04.voxastro.org/~katkov/data/.
| Image (Legacy Survey) | Galaxy | Obs. Date | Positional Angle (deg) | Folder |
|---|---|---|---|---|
 |
NGC 338 | 01.10.2006 | 109 | s061001/ |
 |
NGC 524 | 14.10.2009 | 37 | s091014/ |
 |
NGC 5440 | 12.05.2007 | 50 | s070512/ |
 |
NGC 5533 | 12.05.2007 | 28 | s080505/ |
2.1.1 Ревизия данных¶
!ls -lash data/n5533/raw/s080505/
Директория s080505 содержит все данные, наблюдавшиеся ночь 5-6 мая 2008 года. Лог файл s080505.txt содержит информацию о всех снятых кадрах:
# for some unknown reason converted HTML file has an empty cell
# !cat 'data/n5533/raw/s080505/s080505.txt'
with open('data/n5533/raw/s080505/s080505.txt', 'r') as f:
print(f.read())
Note: Через восклицательный знак можно использовать доступные команды терминала (
!ls,!cat,!unzipи др.).
Не всегда вместе с данными предоставляется лог файл, поэтому полезно уметь составить логфайл на основе информации в FITS заголовках. Сделаем это для тренировки этого навыка:
!unzip 'data/n5533/raw/s080505/*.zip' -d 'data/n5533/raw/s080505/'
В astropy есть удобные command-line tools, в том числе команда fitsheader, которая будет доступна в терминале, если активировано виртуальное python окружение с установленным astropy. fitsheader показывает FITS заголовок файла.
!fitsheader data/n5533/raw/s080505/S5640306.FTS
Также fitsheader может открывать много файлов и, более того, формировать таблицу по заданным ключевым словам из FITS заголовков:
!fitsheader -f -k OBJECT -k IMAGETYP -k MODE -k EXPTIME -k GAIN -k BINNING -k SEEING -k PARANGLE -k ROTANGLE data/n5533/raw/s080505/*.FTS
Мы будем работать со спектром, полученным вдоль большой оси галактики NGC 5533 $PA=28^\circ$. Если в названии объекта не указан позиционный уголь щели, его можно посчитать используя ключевые слова PARANGLE и ROTANGLE: $PA_{slit} = PARANGLE - ROTANGLE + 132.5^\circ$.
Рекомендуется на каждом этапе редукции просматривать все используемые файлы, например, в ds9. Если он установлен и корректно добавлен в переменную окружения PATH (те команда ds9 доступна в терминале), тогда открывать файлы можно прямо из jupyter ноутбука командами типа !ds9 data/n5533/raw/s080505/S5640304.FTS, !ds9 -multiframe data/n5533/raw/s080505/S56403*.FTS.
2.1.2 Построение супербайеса¶
Спектры считывались с матрицы в режиме бинирования 1x2, поэтому bias файлы следует брать с аналогичным бинингом.
import numpy as np
from astropy.io import fits
from astropy import stats
import glob
raw = 'data/n5533/raw/s080505'
wd = 'data/n5533/processed/'
idxs_bias = ['01', '02', '03', '04', '05', '06', '07', '08', '09']
# считываем один кадр, чтобы понять размер кадра
bias0 = fits.getdata(f"{raw}/S5640801.FTS")
shp = bias0.shape
# заполняем стэк
stack_bises = np.zeros((shp[0], shp[1], len(idxs_bias)))
readnoises = np.zeros(len(idxs_bias))
for i, idx in enumerate(idxs_bias):
current_file = f"{raw}/S56408{idx}.FTS"
data = fits.getdata(current_file)
header = fits.getheader(current_file)
# заодно посчитаем статистику по кадрам mean, median, std, std/GAIN и, тем самим, оценим шум считывани
mean, median, std = stats.sigma_clipped_stats(data)
readnoises[i] = std
print(f"Working on: {current_file}: {mean:.3f} {median:.3f} {readnoises[i]:.3f} {readnoises[i]*header['GAIN']:.3f}")
stack_bises[:, :, i] = data
readnoise = np.mean(readnoises)
readnoise_e = np.mean(readnoises*header['GAIN'])
print(f"Mean Read Noise (ADU): {readnoise:.4f}")
print(f"Mean Read Noise (electrons): {readnoise_e:.4f}")
# Считаем супербайес
# используем сигма-клипастое (sigma-clipping) среднее, чтобы убрать редкие отскакивающие отсчеты
superbias, _, _ = stats.sigma_clipped_stats(stack_bises, axis=2)
Для тестов, очень полезно рисовать интерактивные картинки. Для этого следует использовать магическую команду магическая команда %matplotlib notebook. Полезную информацию на эту тему можно найти тут.
После подбора параметров отображения, можно вернуть настройку для рисования статичной картинки прямо в ноутбуке. Она же останется в сохраненном файле ноутбука.
from matplotlib import pyplot as plt
plt.rcParams['figure.figsize'] = [14, 6]
im = plt.imshow(superbias, origin='lower')
plt.colorbar(im)
plt.show()
Чтобы было лучше видно структуру пикселей - добавим нормализацию (подробнее о simple_norm()).
В интерактивном режиме отображения графиков в ноутбуке %matplotlib notebook каждая следующая команда отрисовки будет действовать на первую интерактивную картинку. Чтобы избежать такого поведения нужно нажать на кнопку выключения интерактивного режима в верхнем правом угле.
from astropy.visualization import simple_norm
norm = simple_norm(superbias, 'linear', percent=60.0)
im = plt.imshow(superbias, origin='lower', norm=norm)
plt.colorbar(im)
plt.show()
Видно, что средний байес кадр (супербайес) имеет неравномерную структуру, которая не видна на одиночных кадрах байеса. Это несущественная вариация уровня отсчетом, поскольку оцененный шум считывания $\approx4.7$ ADU.
Наконец сохраним superbias в отдельный файл в директории ./data/n5533/processed/, где будут жить все продукты редукции.
Но перед этим мы зададим область overscan, которую будем "обрезать" на всех кадрах. Чтобы не прогадать с границами overscan лучше всего использовать хорошо освещенные кадры, например, изображение плоского поля.
!ds9 -scale 99.5 data/n5533/raw/s080505/S5640312.FTS
# задаем область оверскана, которая будет использоваться ниже
region_ovescan = [21, 2066, 0, 1025] # X0, X1, Y0, Y1
superbias_cutted = superbias[region_ovescan[2]:region_ovescan[3],
region_ovescan[0]:region_ovescan[1]] # Weird Python indexing (respect to IDL)
shp_cutted = superbias_cutted.shape
fits.PrimaryHDU(data=superbias_cutted, header=header).writeto('data/n5533/processed/superbias.fits', overwrite=True) # use last opened header
!ls -lsrth data/n5533/processed/
!ds9 data/n5533/processed/superbias.fits
2.1.3 Расчет коэффициентов для приведения к плоскому полю¶
!ds9 -scale mode zscale data/n5533/raw/s080505/S564031[2-3].FTS
Усредняем FLAT по двум кадрам.
hdr_flat = fits.getheader('data/n5533/raw/s080505/S5640312.FTS')
flat1 = fits.getdata('data/n5533/raw/s080505/S5640312.FTS')
flat2 = fits.getdata('data/n5533/raw/s080505/S5640313.FTS')
flat1 = flat1.astype(np.float) - superbias
flat2 = flat2.astype(np.float) - superbias
flat = 0.5 * (flat1 + flat2)
norm = simple_norm(flat, 'linear', percent=95.0)
im = plt.imshow(flat, origin='lower', norm=norm)
plt.colorbar(im)
plt.show()
Усредняем по центральной части кадра, чтобы рассчитать кадр нормировочный коэффициентов.
win = 30
y_cent = flat.shape[0]//2
print(f"Central position (px): {y_cent}")
flat_mean_1d = np.median(flat[y_cent-win:y_cent+win, :], axis=0)
plt.plot(flat_mean_1d)
plt.show()
flat_mean_2d = np.resize(flat_mean_1d, flat.shape)
norm = simple_norm(flat_mean_2d, 'linear', percent=95.0)
im = plt.imshow(flat_mean_2d, origin='lower', norm=norm)
plt.colorbar(im)
plt.show()
flat_norm = flat / flat_mean_2d
fits.PrimaryHDU(data=flat_norm, header=header).writeto('data/n5533/processed/flat.fits', overwrite=True)
norm = simple_norm(flat_norm, 'linear', percent=95.0)
im = plt.imshow(flat_norm, origin='lower', norm=norm)
plt.colorbar(im)
plt.show()
2.1.4 Чистка следов космических частиц¶
from lacosmic import lacosmic
files_object = ['S5640306', 'S5640307','S5640309','S5640310']
for file in files_object:
file_path = f"data/n5533/raw/s080505/{file}.FTS"
data, hdr = fits.getdata(file_path, header=True)
data = data.astype(np.float) - superbias # вычитаем супербайес
# создаем кадр ошибок (пуассоновский шум + шум считывания)
error = np.sqrt(np.clip(data,0.0,None)) + readnoise
print(f"Cleaning Cosmic Hits in {file}")
data_clean, mask_ch = lacosmic(data, 2.0, 10, 5.0, effective_gain=hdr['GAIN'], readnoise=readnoise_e)
# flat-fielding
data_out = data_clean / flat_norm
# здесь предполагаем, что нормировочные коэффициенты приведения к плоскому полю не вносят значимого вклада в ошибку
# (В действительности - это может быть не совсем так.)
error_out = error / flat_norm
# remove overscan
data_cutted = data_out[region_ovescan[2]:region_ovescan[3],
region_ovescan[0]:region_ovescan[1]] # Weird Python indexing (respect to IDL)
error_cutted = error_out[region_ovescan[2]:region_ovescan[3],
region_ovescan[0]:region_ovescan[1]]
# write to file
ofile = f"data/n5533/processed/{file}_proc.fits"
ofile_e = f"data/n5533/processed/{file}_proc_err.fits"
print(f"write file: {ofile}")
print(f"write err file: {ofile_e}")
fits.PrimaryHDU(data=data_cutted, header=hdr).writeto(ofile, overwrite=True)
fits.PrimaryHDU(data=error_cutted, header=hdr).writeto(ofile_e, overwrite=True)
# plotting
norm = simple_norm(data, 'log', percent=90.0)
plt.subplot(121)
plt.imshow(data, origin='lower', norm=norm)
plt.subplot(122)
plt.imshow(data_clean, origin='lower', norm=norm)
plt.show()
2.1.5 Калибровка длин волн¶
Подготовка спектра неона¶
Первый шаг - это усреднить кадры неона, по которому будем рассчитывать wavelength solution.
Но сначала проверяем проверяем все кадры неона и смотрим в режиме блинкирования нет ли заметных глазом сдвигов неона (позже рассчитаем этот сдвиг точно).
!ds9 data/n5533/raw/s080505/S56403{05,08,11}.FTS
Сильных сдвигов не видно, но это нужно точнее проверить. Сдвиг неонов, которые снимались в начале, в середине и в конце времени наблюдений объекта, будет свидетельствовать о возможных сдвигах между экспозициями объекта, что тоже следует учитывать.
Вычтем байес, нормируем к плоскому полю, уберем overscan и сформируем стэк из кадров.
files_neon = ['S5640305', 'S5640308','S5640311']
stack_neons = np.zeros((shp_cutted[0], shp_cutted[1], len(files_neon)))
for i, file in enumerate(files_neon):
file_path = f"data/n5533/raw/s080505/{file}.FTS"
print(f"Working on {file_path}")
# read data
data, hdr = fits.getdata(file_path, header=True)
# Subtract superbias and flat-fielding
data = (data.astype(np.float) - superbias) / flat_norm
# remove overscan
data_cutted = data[region_ovescan[2]:region_ovescan[3],
region_ovescan[0]:region_ovescan[1]] # Weird Python indexing (respect to IDL)
stack_neons[:, :, i] = data_cutted
norm = simple_norm(data_cutted, 'log', percent=99.0)
plt.imshow(data_cutted, origin='lower', norm=norm)
plt.show()
Теперь аккуратное разбираемся со сдвигами. TBC!!!
neon = np.mean(stack_neons, axis=-1)
fits.PrimaryHDU(data=neon, header=header).writeto('data/n5533/processed/neon.fits', overwrite=True)
norm = simple_norm(neon, 'log', percent=99.0)
plt.imshow(neon, origin='lower', norm=norm)
plt.show()
Кросс-идентификация линий неона¶
Это важный шаг. Если кросс-идентификация будет неправильная - wavelength solution будет с ошибкой, тогда все измерения скоростей будут неправильными.
Рисуем спектр неона из нижней части кадра и сопоставляем с реперной картинкой, приведенной в первой части описания. Записываем положения известных линий неона. Здесь удобно использовать интерактивный режим отображения matplotlib графиков.
%matplotlib notebook
neon_1d = np.mean(neon[0:50, :], axis=0)
plt.plot(neon_1d)
plt.show()
line_wav0 = np.array([4837.48, 4863.96, 4885.10, 4921.93 , 4956.89 , 5015.68 , 5037.75, 5047.74, 5080.38, 5116.5, 5145.00, 5152.10, 5162.28, 5188.3, 5203.9, 5221.91, 5330.78, 5341.6, 5400.56, 5433.65, 5495.52])
line_pix0 = np.array([114. , 193. , 257. , 361. , 463. , 627. , 688. , 715. , 804. , 901. , 976.4 , 995 , 1022. , 1090. , 1132. , 1178. , 1459. , 1486. , 1636. , 1718 , 1874. ])
# маскируем некоторые линии, которые после проверки в следующих пунктах оказались выпадающими, как правило из-за бленд
line_use = np.array([ 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1])
msk = np.array(line_use) == 1
line_wav0 = line_wav0[msk]
line_pix0 = line_pix0[msk]
nlines = line_wav0.size
%matplotlib inline
plt.plot(line_pix0, line_wav0, '+')
plt.show()
Сделаем начальное приближение еще более точным.
coef = np.polyfit(line_wav0, line_pix0, deg=4)
xx = np.arange(shp_cutted[1])
yy = np.polyval(coef, xx)
# обновленные значения положений линий
line_pix = np.polyval(coef, line_wav0)
# посмотрим насоклько промазали
plt.plot(line_wav0, line_pix0 - line_pix, '+')
plt.ylabel('$\Delta X, \AA$')
plt.axhline(0)
plt.show()
Имея в руках положение линий нужно определить положение линий с большей точностью и рассчитать положение этих линий по всему кадру.
from scipy import optimize
def gaus_cont(x, cent, sigma, int0, c0, c1):
"Simple gaussian + continuum"
yy = (x - cent) / sigma
return c0 + c1*x + int0*np.exp(-0.5*yy**2)
def find_peaks(array, guess_positions, window=15.0, sigma=3.0, plot=False):
"Accurately find centers of lines fitting them by Gaussian line + continuum"
xx = np.arange(array.size)
pos_out = np.array(guess_positions)
for i, pos0 in enumerate(guess_positions):
lims = np.array([pos0-window, pos0+window]).astype(np.int)
y = array[lims[0]:lims[1]]
x = xx[lims[0]:lims[1]]
popt, pcov = optimize.curve_fit(gaus_cont, x, y, p0=[pos0, sigma, np.nanmax(y), 0.0, 0.0])
pos_out[i] = popt[0]
if plot:
print(popt)
fit = gaus_cont(x, *popt)
plt.plot(x, y)
plt.plot(x, fit, color='red')
plt.axvline(popt[0], linestyle='--')
plt.show()
return pos_out
Посмотрим как фиты линий выглядят в одной "строке" неона.
position0 = find_peaks(neon[0, :], line_pix, window=15.0, sigma=2.0, plot=True)
Теперь пробегаемся по всем строкам и считаем положения линий неона. При этом в каждой строке берем начальное приближении с предыдущей.
from tqdm import tqdm
lines_2d_X = np.zeros((nlines, shp_cutted[0]))
lines_2d_Y = np.zeros((nlines, shp_cutted[0]))
lines_2d_wav = np.zeros((nlines, shp_cutted[0]))
for j in tqdm(range(shp_cutted[0])):
lines_2d_X[:, j] = find_peaks(neon[j, :], position0, window=15.0, sigma=2.0, plot=False)
position0 = lines_2d_X[:, j]
lines_2d_wav[:, j] = line_wav0
lines_2d_Y[:, j] = j
Чтобы уменьшить размер данных, можно предварительно пробинировать спектр неона по 3, 5, 10 пикселей. Но для простоты оставим анализ в каждом пикселей неона, благо неон как правило снимается с хорошим "сигнал-шум" и можно работать с построчным спектром.
Теперь нужно проверить корректность определения положения линий по кадру неона.
xx = np.arange(shp_cutted[0])
plt.figure(figsize=(16, 6))
norm = simple_norm(neon, 'log', percent=99.0)
plt.imshow(neon, origin='lower', norm=norm)
for k in range(nlines):
plt.plot(lines_2d_X[k, :], xx, '.', ms=0.5, color='red')
plt.show()
Определение wavelength solution¶
Выглядит адекватно! Теперь нужно найти желанное wavelength solution. Будем искать зависимость $\lambda = f(x, y)$ в виде двумерного полинома по оси X (вдоль дисперсии) степени 5, а вдоль щели 4: $\lambda(x, y) = \sum_{i=0}^{5} \sum_{j=0}^{3} a_{i,j} x^i y^j$.
import itertools
# Функция рассчета коэффициентов полиномиального ряда.
def polyfit2d(x, y, z, deg=[4, 5]):
"Least squares for the polynomial parametrizations of the surface."
ncols = (deg[0] + 1) * (deg[1] + 1)
matrix = np.zeros((x.size, ncols))
ij = itertools.product(range(deg[0] + 1),
range(deg[1] + 1))
for k, (i,j) in enumerate(ij):
matrix[:, k] = x**j * y**i
c = np.linalg.lstsq(matrix, z, rcond=None)[0]
return c.reshape((deg[0]+1, deg[1]+1))
# Функция для рассчета суммы полиномов в произвольной точке, используя рассчитанные коэффициенты
def polyval2d_norm(x, y, coef):
"Evaluate polynomial surface"
ij = itertools.product(range(coef.shape[0]),
range(coef.shape[1]))
z = np.zeros_like(x)
for i,j in ij:
z += coef[i, j] * x**j * y**i
return z
# Функция обертка для работы с пиксельными координатами и ангстремами, вместо нормированных массивов,
# которые используются в рассчете
def polyval2d(x, y, coef, x_minmax, y_minmax, z_minmax):
x_norm = (x - x_minmax[0]) / (x_minmax[1] - x_minmax[0])
y_norm = (y - y_minmax[0]) / (y_minmax[1] - y_minmax[0])
z_norm = polyval2d_norm(x_norm, y_norm, coef)
z = z_norm * (z_minmax[1] - z_minmax[0]) + z_minmax[0]
return z
deg = [4, 5]
xx1d = np.arange(shp_cutted[1], dtype=np.float)
yy1d = np.arange(shp_cutted[0], dtype=np.float)
xx2d, yy2d = np.meshgrid(xx1d, yy1d)
# Немного упрощаем рассчет полиномиальной зависимости, используя лишь каждую 20ю строчку измерений положений линий по неону.
# И разворачиваем массивы в 1d.
xp = lines_2d_X[:, ::20].flatten() # измеренные положения линий
yp = lines_2d_Y[:, ::20].flatten() # положение спектра неона вдоль щели
zp = lines_2d_wav[:, ::20].flatten() # значение длины волны идентифицированной линии
# Чтобы избежать связанных с численной точностью проблем - нормализуем данные.
# Если этого не сделать, коэффициенты при МНК фите будут неправильными.
xp_minmax = [np.min(xp), np.max(xp)]
yp_minmax = [np.min(yp), np.max(yp)]
zp_minmax = [np.min(zp), np.max(zp)]
xp_norm = (xp - xp_minmax[0]) / (xp_minmax[1] - xp_minmax[0])
yp_norm = (yp - yp_minmax[0]) / (yp_minmax[1] - yp_minmax[0])
zp_norm = (zp - zp_minmax[0]) / (zp_minmax[1] - zp_minmax[0])
# Делаем фит и рассчитываем коэффициенты
coef = polyfit2d(xp_norm, yp_norm, zp_norm, deg=deg)
print(coef)
# Для проверки смотрим на значения chi2 и средний разброс
fit = polyval2d_norm(xp_norm, yp_norm, coef)
fit_a = polyval2d(xp, yp, coef, xp_minmax, yp_minmax, zp_minmax)
print(f"chi2: {np.sum((zp_norm - fit)**2)}")
print(f"chi2/dof: {np.sum((zp_norm - fit)**2)/x.size}")
print(f"STD: {np.std(zp - fit_a):.4f} Å")
# Рассчитываем длины волн в каждом пикселе кадра
full_wave_solution = polyval2d(xx2d, yy2d, coef, xp_minmax, yp_minmax, zp_minmax)
Для проверки смотрим на несколько картинок.
Оранжевые крестики - показывают положения идентифицированных линий, на которые мы аппроксимируем. Белые линии отображают линии равных длин волн полученной аппроксимации. Кривизна линий примерно совпадает, что говорит об отсутствии грубой ошибки.
plt.figure(figsize=(16, 8))
im = plt.imshow(full_wave_2d, origin='lower')
plt.colorbar(im)
cs = plt.contour(xx1d, yy1d, full_wave_solution, colors='white', levels=25)
plt.clabel(cs, inline=1, fontsize=10, fmt='%.0f')
# plt.scatter(xp, yp, c='white', marker='+', s=40, edgecolors='k', linewidths=0.3)
plt.scatter(xp, yp, c='C1', marker='+', s=40)
plt.show()
Валидация найденного решения¶
Теперь более детально посмотрим на качество полученной аппроксимации и оценим точность калибровки длин волн.
for j in range(nlines):
fit = polyval2d(lines_2d_X[j, ::20], lines_2d_Y[j, ::20], coef, xp_minmax, yp_minmax, zp_minmax)
resid = lines_2d_wav[j, ::20] - fit
resid_std = np.std(resid)
resid_mean = np.mean(resid)
resid_median = np.median(resid)
print(f"Line {line_wav0[j]:.2f} | STD: {resid_std:5.3f} Å | Mean: {resid_mean:6.3f} Å | Median: {resid_median:6.3f} Å")
plt.figure(figsize=(16, 3))
plt.subplot(121)
plt.plot(lines_2d_Y[j, ::20], resid, '+')
plt.axhline(0.0, linestyle='--', color='k')
plt.ylabel("$\Delta \lambda, \AA$")
plt.xlabel("Y position, px")
plt.subplot(122)
line_position_X = lines_2d_X[j, 500]
lambda_position = polyval2d(line_position_X, 500, coef, xp_minmax, yp_minmax, zp_minmax)
x0, x1 = int(line_position_X)-30, int(line_position_X)+30
plt.plot(xx1d[x0:x1], neon[500, x0:x1])
plt.axvline(line_position_X, linestyle=':')
plt.title(f"{line_wav0[j]:3f}; solution: {lambda_position:3f}")
plt.show()
Ошибка идентификации линий и длинноволнового решения составляет $\approx0.02$ Å, что вполне хорошо для решетки 2300G прибора SCORPIO.
2.1.6 Линеаризация спектров и сложение экспозиций¶
files_object = ['S5640306', 'S5640307','S5640309','S5640310']
stack_lin = np.zeros((neon.shape[0], neon.shape[1], len(files_object)))
stack_err = np.zeros_like(stack_lin)
stack_exptime = np.zeros(len(files_object))
for ifile, file in enumerate(files_object):
file_path = f"data/n5533/processed/{file}_proc.fits"
e_file_path = f"data/n5533/processed/{file}_proc_err.fits"
data, hdr = fits.getdata(file_path, header=True)
err = fits.getdata(e_file_path)
data_new = data.copy()
err_new = err.copy()
wave_new = np.linspace(full_wave_2d[500, 0], full_wave_2d[500,-1], num=full_wave_2d.shape[1])
for j in tqdm(range(data.shape[0]), ):
data_new[j, :] = np.interp(wave_new, full_wave_2d[j, :], data[j, :])
err_new[j, :] = np.interp(wave_new, full_wave_2d[j, :], err[j, :])
stack_lin[:, :, ifile] = data_new
stack_err[:, :, ifile] = err_new
stack_exptime[ifile] = hdr['EXPTIME']
ofile = f"data/n5533/processed/{file}_lin.fits"
ofile_e = f"data/n5533/processed/{file}_lin_err.fits"
print(f"write file: {ofile}")
print(f"write file: {ofile_e}")
hdr['CRPIX1'] = 1
hdr['CRVAL1'] = wave_new[0]
hdr['CDELT1'] = wave_new[1] - wave_new[0]
hdr['CTYPE1'] = 'AWAVE-LIN'
fits.PrimaryHDU(data=data_new, header=hdr).writeto(ofile, overwrite=True)
fits.PrimaryHDU(data=err_new, header=hdr).writeto(ofile_e, overwrite=True)
# plotting
norm = simple_norm(data, 'log', percent=99.0)
plt.figure(figsize=(16,4))
plt.subplot(121)
plt.imshow(data, origin='lower', norm=norm)
plt.subplot(122)
plt.imshow(data_new, origin='lower', norm=norm)
plt.show()
Из-за технических особенностей гризма в приборе может быть установлена под небольшим углом по отношению к матрице, что приводит к небольшому наклону направления дисперсии по отношению к X оси матрица. Этот эффект можно проверить по положению объекта в ярким центром или по звезде в кадре. В приборе SCORPIO2 используется 13 точечный тест - это кадр флэта сделанный с маской (расположенной на щели). По искажениям этих 13 полос (спектров флэта) на спектре можно восстановить геометрические искажения в оптической системе и их исправить (подробнее см. в Лекция III. Многорежимный фокальный редуктор телескопа БТА (Моисеев А.В.). Но в нашем случае проверим сдвиг по самим кадрам.
Сначала простая графическая проверка.
plt.figure(figsize=(16, 4))
yrange = [460, 580]
for i in range(4):
plt.subplot(1,4,i+1)
spec_sub = stack_lin[yrange[0]:yrange[1],:, i]
norm = simple_norm(spec_sub, 'linear', percent=99.0)
plt.imshow(spec_sub, origin='lower', norm=norm, aspect=20)
plt.show()
Несмотря на сильное сжатие картинок (aspect=20), подписи к осям соответствуют реальным пикселям на кадрах. На каждом кадре четко видно, как центр галактики "едет" по диапазону на ~10 px. Этот эффект обязательно нужно исправлять.
Для этого нужно определить положение центра галактики в каждом X и потом сдвинуть спектр. Не забываем сдвигать кадры ошибок.
from matplotlib import gridspec
plt.figure(figsize=(16, 8))
gs = gridspec.GridSpec(2,4)
yrange = [460, 580]
xx = np.arange(stack_lin.shape[1])
yy = np.arange(stack_lin.shape[0])
# Это будет положение центра галактики во всех кадрах после корреции.
# Таким образом мы сразу учитываем возможный сдвиг положения центра галактики от кадра к кадру.
ycenter = stack_lin.shape[0] // 2
stack_lin_centered = np.zeros_like(stack_lin)
stack_err_centered = np.zeros_like(stack_err)
for i, file in enumerate(files_object):
spec_sub = stack_lin[yrange[0]:yrange[1],:, i]
ymax = np.argmax(spec_sub, axis=0)
coef = np.polyfit(xx, ymax, deg=1)
print(f"Poly Coef.: {coef}")
ymax_poly = np.polyval(coef, xx)
# рисуем положение центров в каждой колонке
plt.subplot(gs[0, i])
norm = simple_norm(spec_sub, 'log', percent=99.0)
plt.imshow(spec_sub, origin='lower', norm=norm, aspect=20)
plt.plot(ymax, ',')
plt.plot(ymax_poly)
# Исправлляем сдвиг всего кажра и рисуем кусочек кадра, чтобы лучше было видно
for q in range(stack_lin.shape[1]):
delta_y = ycenter - (ymax_poly[q] + yrange[0])
stack_lin_centered[:, q, i] = np.interp(yy, yy + delta_y, stack_lin[:, q, i])
stack_err_centered[:, q, i] = np.interp(yy, yy + delta_y, stack_err[:, q, i])
spec_sub = stack_lin_centered[ycenter-50:ycenter+50,:, i]
plt.subplot(gs[1, i])
plt.imshow(spec_sub, origin='lower', norm=norm, aspect=20)
# write fits files
hdr['CRPIX2'] = ycenter
hdr['CRVAL2'] = 0.0
hdr['CDELT2'] = np.float(hdr['IMSCALE'].split('x')[1])
ofile = f"data/n5533/processed/{file}_centered.fits"
e_ofile = f"data/n5533/processed/{file}_centered_err.fits"
print(f"write file: {ofile}")
print(f"write file: {e_ofile}")
fits.PrimaryHDU(data=stack_lin_centered[:, :, i], header=hdr).writeto(ofile, overwrite=True)
fits.PrimaryHDU(data=stack_err_centered[:, :, i], header=hdr).writeto(e_ofile, overwrite=True)
plt.show()
Отлично! Теперь все индивидуальный кадры отцентрованы по Y и можно их складывать в одну суммарную экспозицию.
objlin = np.sum(stack_lin_centered, axis=2)
objerr = np.sqrt(np.sum(stack_err_centered**2, axis=2))
exptime = np.sum(stack_exptime)
hdr['EXPTIME'] = exptime
objfile = f"data/n5533/processed/obj_lin.fits"
e_objfile = f"data/n5533/processed/obj_err.fits"
print(f"write file: {objfile}")
print(f"write file: {e_objfile}")
fits.PrimaryHDU(data=objlin, header=hdr).writeto(objfile, overwrite=True)
fits.PrimaryHDU(data=objerr, header=hdr).writeto(e_objfile, overwrite=True)
2.1.7 Вычитание спектра неба¶
Как видно на изображениях спектров выше, в анализируемом спектре - небо очень яркое - видны как эмиссионная линия неба $\approx5460$ Å, так и континуум, на котором видны линии триплета магния в рассеянном солнечном спектре. Короче говоря, для этой галактики вычитание небо особенно важный этап.
Оценим область на щели, где нет галактики.
prof = np.nanmedian(objlin[:, 700:1500], axis=1)
plt.plot(prof)
plt.yscale('log')
plt.xlabel('Position along the slit, px')
plt.show()
Видно, что галактика на самом деле простирается далеко. Поэтому область для вычитания неба придется брать далеко. Принимая во внимание, что SCORPIO имеет сильные вариации инструментального контура (что хорошо видно на спектре неона как линии становятся асимметричными, у них появляется сильное "голубое крыло" на краю кадра) - можно ожидать, что вычитание неба будет неидеальным.
sky_regions = [[10, 180], [800, 1000]]
yy = np.arange(objlin.shape[0])
sky_msk_y = np.zeros(objlin.shape[0], dtype=np.bool)
for sky_reg in sky_regions:
sky_msk_y[(yy >= sky_reg[0]) & (yy <= sky_reg[1])] = True
plt.plot(yy, prof)
plt.plot(yy[sky_msk_y], prof[sky_msk_y], 'o')
plt.yscale('log')
plt.xlabel('Position along the slit, px')
plt.show()
sky = np.zeros_like(objlin)
sky_deg = 2
for j in tqdm(range(objlin.shape[1])):
sky_coef = np.polyfit(yy[sky_msk_y], objlin[sky_msk_y, j], deg=sky_deg)
sky[:, j] = np.polyval(sky_coef, yy)
objsky = objlin - sky
# plotting
plt.figure(figsize=(16,4))
plt.subplot(121)
norm = simple_norm(objsky, 'log', percent=80.0)
plt.imshow(objsky, origin='lower', norm=norm)
plt.subplot(122)
norm = simple_norm(sky, 'log', percent=99.05)
plt.imshow(sky, origin='lower', norm=norm)
plt.show()
# writing to file
f_obj_sky = f"data/n5533/processed/obj_sky.fits"
print(f"write file: {f_obj_sky}")
fits.PrimaryHDU(data=objsky, header=hdr).writeto(f_obj_sky, overwrite=True)
f_sky = f"data/n5533/processed/sky.fits"
print(f"write file: {f_sky}")
fits.PrimaryHDU(data=sky, header=hdr).writeto(f_sky, overwrite=True)
Вертикальные полоски во внешней части галактики (лучше видны в ds9 при удачном подборе яркости) как раз являются артефактами вычитания ночного неба. В определении параметров галактики - они в первую очередь влияют на металличность звездного населения, но влияние на скорости незначительно. Эти артефакты влияют на определение дисперсии скоростей звезд, особенно при низком сигнал-шум, когда вырождения металличность-дисперсия существеннее себя проявляет. Но в целом, это влияние на конечную оценку не столь критична, как, например, точное измерение инструментального контура.
Итак, у нас в руках отредуцированный суммарный спектр галактики, в котором отлично видна (лучше в ds9) структура абсорбционных линий, но кроме того видна протяженная линия $H\beta$ и чуть слабее [OIII] 5007 Å. Последний кламп $H\beta$ виден в $\approx300\mathrm{px}\approx110''$ от центра галактики, т.е. можно протянуть кривую вращения по газу вплоть до этого расстояния!
Но для начала оценим SNR по области спектра, где нет эмиссионных линий или артефактов вычитания неба, скажем в области $5190\pm20$ Å.
spec, hdr = fits.getdata('data/n5533/processed/obj_sky.fits', header=True)
error = fits.getdata('data/n5533/processed/obj_err.fits')
wave = hdr['CRVAL1'] + hdr['CDELT1']*np.arange(spec.shape[1])
msk_snr = np.abs(wave - 5190) <= 20
signal = np.sum(spec[:, msk_snr], axis=1)
noise = np.sum(error[:, msk_snr], axis=1)
plt.plot(signal/noise)
plt.xlabel("Y position, px")
plt.ylabel("SNR")
plt.show()
from vorbin.voronoi_2d_binning import voronoi_2d_binning
arcsec_per_pixel = 0.357 # according IMSCALE = '0.178x0.357'
def binning_linear(signal, noise, minbin=1, maxbinsize=10, maxdist=200):
"Helper function to make linearly increased bins"
binNum = np.zeros_like(signal, dtype=np.int) - 1
i0 = np.argmax(signal/noise)
grad = maxbinsize/maxdist
icurrent = i0
ibin = 0
while icurrent <= i0+maxdist:
binsize = int(1 + grad*(icurrent-i0))
binNum[icurrent:icurrent+binsize] = ibin
ibin += 1
icurrent += binsize
icurrent = i0
while icurrent >= i0-maxdist:
binsize = int(1 + grad*(i0-icurrent))
binNum[icurrent-binsize:icurrent] = ibin
ibin += 1
icurrent -= binsize
nbins = np.max(binNum)+1
xBin = np.zeros(nbins)
sn = np.zeros(nbins)
xx = np.arange(signal.size)
for ibin in range(nbins):
idx_bin = binNum == ibin
xBin[ibin] = np.nansum(xx[idx_bin]*signal[idx_bin]) / np.nansum(signal[idx_bin])
sn[ibin] = np.nansum(signal[idx_bin]) / np.sqrt(np.nansum(noise[idx_bin]**2))
return binNum, xBin, sn
binNum, xBin, sn = binning_linear(signal, noise, minbin=2, maxbinsize=70, maxdist=300)
imax = np.argwhere(binNum == 0)[0][0]
print(f"Nbins: {np.max(binNum)+1}")
# plotting
plt.figure(figsize=(16,4))
plt.subplot(121)
plt.plot(signal/noise)
plt.xlabel("Y position, px")
plt.ylabel("SNR")
# каждый bin отрисуем отдельным цветом
for i in range(np.max(binNum)+1):
bin_idx = (binNum == i)
color = (np.random.uniform(size=3)*255).astype(np.int)
color = np.random.uniform(size=3)
plt.axvspan(np.min(x[bin_idx]), np.max(x[bin_idx]), alpha=0.3, color=color)
ax2 = plt.subplot(122)
plt.plot(xBin, sn, 'o')
secax = ax2.secondary_xaxis('top', functions=(lambda x: (x-imax)*arcsec_per_pixel, lambda x: x/arcsec_per_pixel+imax))
ax2.set_xlabel('Y position, px')
secax.set_xlabel('Raius, arcsec')
ax2.axvline(imax, linestyle=':')
plt.show()
2.2.2 Моделирование спектра¶
Для моделирования спектра мы будем использовать пакет pPXF.
После установки пакета pPXF в директории site-packages/ppxf/examples будут находиться полезные примеры использования pPXF. В качестве опоры для нашего анализа мы будем использовать пример site-packages/ppxf/examples/ppxf_example_population_gas_sdss.py.
from os import path
from time import perf_counter as clock
import ppxf as ppxf_package
from ppxf.ppxf import ppxf
import ppxf.ppxf_util as util
import ppxf.miles_util as lib
ppxf_dir = path.dirname(path.realpath(ppxf_package.__file__))
print(ppxf_dir)
Подготовка спектра галактики¶
В редукции мы линеаризовали спектр и привели его к линейной шкале длин волн. Для анализа нам нужно использовать логбинированный спектр (т.е. у которого шкала длин волн логарифмическая). Чтобы избежать дополнительной интерполяции данных мы могли бы на шаге линеаризации сразу интерполировать спектр в логарифмическую шкалу $\lambda$. Но в практике спектры чаще встречаются в линейной шкале, поэтому разумно сделать логбинирование средствами доступными в pPXF.
galaxy_lin, hdr = fits.getdata('data/n5533/processed/obj_sky.fits', header=True)
# здесь по умолчанию предполагается, что реперный пиксель CRPIX=1. Но так бывает не всегда.
lamRange1 = hdr['CRVAL1'] + np.array([0., hdr['CDELT1']*(hdr['NAXIS1'] - 1)])
galaxy = np.zeros_like(galaxy_lin)
noise = np.zeros_like(galaxy_lin)
for j in range(galaxy_lin.shape[0]):
gal, logLam1, velscale = util.log_rebin(lamRange1, galaxy_lin[j, :])
err, _, _ = util.log_rebin(lamRange1, error[j, :])
galaxy[j, :] = gal
noise[j, :] = err
wave = np.exp(logLam1)
print(f"Velscale (km/s): {velscale[0]:.3f}")
В логбинированном спектре сдвиг по X в пикселях эквивалентен сдвигу по скорости. velscale - это размер пикселя в км/с. Опыт показывает, что при достойном SNR можно уверенно задетектировать сдвиги $\approx 0.1$ px, а в редких случаях $\approx 0.03-0.05$ px. Отсюда характерная точность измерения скорости при приличном SNR будет $\approx 0.1\mathrm{px} \approx 2 \mathrm{km/s}$.
Подготовка темплейтов¶
Так же pPXF идет с библиотекой звезд MILES и библиотекой SSP моделей MILES на основе звездной библиотеки. Спектральное разрешение MILES $2.51$ Å, что близко к спектральному разрешению наших данных.
Настраиваем темплейты по-примеру ppxf_example_population_gas_sdss.py. В рассматриваемой галактике кроме звездного спектра галактики видны эмиссионные линии $H\beta$ и [OIII]. pPXF позволяет одновременно фитировать и спектр галактики и эмиссии, но для этого нужно определенным образом сформировать стэк темплейтов, в который кроме SSP моделей будут входит гауссовы линии для фита эмиссий.
c = 299792.458 # speed of light in km/s
FWHM_gal = 2.5 # Примерное разрешение данных SCORPIO
vsys = 3863. # взято из HyperLeda http://leda.univ-lyon1.fr/ledacat.cgi?o=NGC5533
pathname = ppxf_dir + '/miles_models/Mun1.30*.fits'
# The templates are not normalized.
# In this way the weights and mean values are mass-weighted quantities.
# Use the keyword 'norm_range' for light-weighted quantities.
miles = lib.miles(pathname, velscale, FWHM_gal)
# The stellar templates are reshaped below into a 2-dim array with each
# spectrum as a column, however we save the original array dimensions,
# which are needed to specify the regularization dimensions
#
stars_templates = miles.templates.reshape(miles.templates.shape[0], -1)
# Estimate the wavelength fitted range in the rest frame.
lam_range_gal = np.array([np.min(wave), np.max(wave)])/(1 + vsys/c)
# Construct a set of Gaussian emission line templates.
# The `emission_lines` function defines the most common lines, but additional
# lines can be included by editing the function in the file ppxf_util.py.
gas_templates, gas_names, line_wave = \
util.emission_lines(miles.log_lam_temp, lam_range_gal, FWHM_gal,
limit_doublets=True)
# Combines the stellar and gaseous templates into a single array.
# During the PPXF fit they will be assigned a different kinematic
# COMPONENT value
#
templates = np.column_stack([stars_templates, gas_templates])
Посмотрим как выглядят используемые темплейты звездных SSP моделей библиотеки MILES. Отрисуем модели солнечной металличности и разных возрастов.
plt.figure(figsize=(16, 7))
miles_wave = np.exp(miles.log_lam_temp)
for iage in [7, 12, 17, 21, 24]:
spec = miles.templates[:, iage, 4]
age = miles.age_grid[iage, 0]
spec /= np.median(spec[5500:5600])
plt.plot(miles_wave, spec, label=f"{age:.1f} Gyr", lw=0.5)
plt.xlabel("Wavelength, Å")
plt.legend()
plt.show()
Моделирование спектра¶
Последние приготовления и можно запустить моделирование спектра галактики.
Следуя примеру ppxf_example_population_gas_sdss.py наша модель будет состоять из 3-х компонент со своей индивидуальной кинематикой: звездное население, бальмеровская линия $H\beta$ и запрещенная линия [OIII].
Перед запуском на всех бинах, стоит сделать тестовый запуск на одном, чтобы подобрать нужные параметры. Ниже приводится пример запуска уже для все бинов.
# dv - это искуственный сдвиг скорости за счет разного размера спектра галактики и темплейтов
dv = c*(miles.log_lam_temp[0] - np.log(wave[0]))
vel = c*np.log(1 + vsys/c) # eq.(8) of Cappellari (2017)
start = [vel, 100.] # (km/s) начальные приближения для [V, sigma]
n_temps = stars_templates.shape[1]
n_forbidden = np.sum(["[" in a for a in gas_names]) # forbidden lines contain "[*]"
n_balmer = len(gas_names) - n_forbidden
# Assign component=0 to the stellar templates, component=1 to the Balmer
# gas emission lines templates and component=2 to the forbidden lines.
component = [0]*n_temps + [1]*n_balmer + [2]*n_forbidden
gas_component = np.array(component) > 0 # gas_component=True for gas templates
# Fit (V, sig) moments=2 for the stars
# and (V, sig) moments=2 for the two gas kinematic components
moments = [2, 2, 2]
# Adopt the same starting value for the stars and the two gas components
start = [start, start, start]
# See the pPXF documentation for the keyword REGUL,
# regul_err = 0.013 # Desired regularization error
# regul = 1./regul_err
reg_dim = miles.templates.shape[1:]
regul = 0.0 # выключаем регуляризацию
# Here the actual fit starts.
#
# IMPORTANT: Ideally one would like not to use any polynomial in the fit
# as the continuum shape contains important information on the population.
# Unfortunately this is often not feasible, due to small calibration
# uncertainties in the spectral shape. To avoid affecting the line strength of
# the spectral features, we exclude additive polynomials (DEGREE=-1) and only use
# multiplicative ones (MDEGREE=10). This is only recommended for population, not
# for kinematic extraction, where additive polynomials are always recommended.
#
vel_stars = np.zeros(np.max(binNum)+1)
vel_Hbeta = np.zeros_like(vel_stars)
vel_OIII = np.zeros_like(vel_stars)
e_vel_stars = np.zeros_like(vel_stars)
e_vel_Hbeta = np.zeros_like(vel_stars)
e_vel_OIII = np.zeros_like(vel_stars)
sig_stars = np.zeros_like(vel_stars)
sig_Hbeta = np.zeros_like(vel_stars)
sig_OIII = np.zeros_like(vel_stars)
e_sig_stars = np.zeros_like(vel_stars)
e_sig_Hbeta = np.zeros_like(vel_stars)
e_sig_OIII = np.zeros_like(vel_stars)
age_stars = np.zeros_like(vel_stars)
met_stars = np.zeros_like(vel_stars)
ml_stars = np.zeros_like(vel_stars)
chi2 = np.zeros_like(vel_stars)
t0 = clock()
for i, ibin in enumerate(range(np.max(binNum)+1)):
# for i, ibin in enumerate(range(3)):
idx_bin = (binNum == ibin)
if np.sum(idx_bin) == 1:
galaxy_bin = galaxy[idx_bin, :].flatten()
noise_bin = noise[idx_bin, :].flatten()
else:
galaxy_bin = np.sum(galaxy[idx_bin, :], axis=0)
noise_bin = np.sqrt(np.sum(noise[idx_bin, :]**2, axis=0))
t = clock()
pp = ppxf(templates, galaxy_bin, noise_bin, velscale, start,
moments=moments, degree=-1, mdegree=10, vsyst=dv,
lam=wave, clean=False, reg_dim=reg_dim, regul=regul,
component=component, gas_component=gas_component,
gas_names=gas_names)
print('Elapsed time in given bin: %.2f s' % (clock() - t))
weights = pp.weights[~gas_component] # Exclude weights of the gas templates
weights = weights.reshape(reg_dim)/weights.sum() # Normalized
logage, met = miles.mean_age_metal(weights)
ml_r = miles.mass_to_light(weights, band="r")
vel_stars[i] = pp.sol[0][0]
vel_Hbeta[i] = pp.sol[1][0]
vel_OIII[i] = pp.sol[2][0]
e_vel_stars[i] = pp.error[0][0]
e_vel_Hbeta[i] = pp.error[1][0]
e_vel_OIII[i] = pp.error[2][0]
sig_stars[i] = pp.sol[0][1]
sig_Hbeta[i] = pp.sol[1][1]
sig_OIII[i] = pp.sol[2][1]
e_sig_stars[i] = pp.error[0][1]
e_sig_Hbeta[i] = pp.error[1][1]
e_sig_OIII[i] = pp.error[2][1]
age_stars[i] = 10**logage/1e9 # in Gyr
met_stars[i] = met
ml_stars[i] = ml_r
chi2[i] = pp.chi2
# Plot fit results for stars and gas.
from matplotlib import gridspec
gs = gridspec.GridSpec(1, 2, width_ratios=[2.7, 1])
plt.figure(figsize=(16, 4))
plt.subplot(gs[0])
pp.plot()
# Plot stellar population mass-fraction distribution
plt.subplot(gs[1])
miles.plot(weights)
plt.show()
print('Elapsed time in PPXF: %.2f s' % (clock() - t))
# коррекция ошибок на фактор sqrt(chi2)
for prof in (e_vel_stars, e_vel_Hbeta, e_vel_OIII,
e_sig_stars, e_sig_Hbeta, e_sig_OIII):
prof *= np.sqrt(chi2)
В идеальном случае, когда входные данные имеют адекватные ошибки, а модель действительно описывает эти данные (нет так называемой проблемы mismatch template), нормированный $\chi^2$ фита должен быть $\approx1$. При большом количестве ре-интерполяций кадров, статистика ошибок на кадре меняется и можно ожидать чуть меньшего значения $\chi^2\approx0.7...0.9$. Если $\chi^2$ сильно больше единицы - это означает либо модель плохая, или в данных есть выбросы (outliers), либо ошибки недооценены, а при $\chi^2<<1$ ошибки скорее переоценены. Таким образом, если мы ожидаем, что расхождения $\chi^2$ от единицы связаны с недо- или переоценкой ошибок, можно использовать корректирующий фактор $f_{err}=\sqrt{\chi^2}$ (здесь предполагается нормированный $\chi^2$). Поскольку форма вектора ошибок не меняется, т.е. не изменяется относительные веса разных точек спектра - фит будет идентичен, изменятся лишь ошибки выходных параметров фита на фактор $f_{err}$. Поэтому в последних строках предыдущего блока кода мы корректировали ошибки полученных параметров фактор $f_{err}$.
Представим профили полученных параметров графически.
# Plotting
arcsec_per_pixel = 0.357 # according IMSCALE = '0.178x0.357'
rsec = (xBin - xBin[0])*arcsec_per_pixel
v0 = vel_stars[0]
xlim = [-100, 100]
xlabel = 'R, arcsec'
plt.figure(figsize=(16, 5))
# sorted indexes
srt = np.argsort(rsec)
ax = plt.subplot(121)
plt.fill_between(rsec[srt], vel_OIII[srt]-e_vel_OIII[srt], vel_OIII[srt]+e_vel_OIII[srt], color='C0', alpha=0.1)
plt.plot(rsec[srt], vel_OIII[srt], marker='.', linestyle='-', label=r'[OIII]', color='C0')
plt.fill_between(rsec[srt], vel_Hbeta[srt]-e_vel_Hbeta[srt], vel_Hbeta[srt]+e_vel_Hbeta[srt], color='C2', alpha=0.1)
plt.plot(rsec[srt], vel_Hbeta[srt], marker='.', linestyle='-', label=r'$H\beta$', color='C2')
plt.fill_between(rsec[srt], vel_stars[srt]-e_vel_stars[srt], vel_stars[srt]+e_vel_stars[srt], color='k', alpha=0.1)
plt.plot(rsec[srt], vel_stars[srt], marker='.', linestyle='-', label=r'Stars', color='k')
plt.legend()
ax.secondary_yaxis('right', functions=(lambda x: x-v0, lambda x: x+v0))
plt.axhline(v0, linestyle=':')
plt.axvline(0.0, linestyle=':')
plt.ylim(v0-350, v0+350)
plt.xlim(xlim)
plt.xlabel(xlabel)
plt.title('$v_\mathrm{LOS}$, km/s')
plt.subplot(122)
plt.plot(rsec, sig_stars, 'k.')
plt.fill_between(rsec[srt], sig_OIII[srt]-e_sig_OIII[srt], sig_OIII[srt]+e_sig_OIII[srt], color='C0', alpha=0.1)
plt.plot(rsec[srt], sig_OIII[srt], marker='.', linestyle='-', label=r'[OIII]', color='C0')
plt.fill_between(rsec[srt], sig_Hbeta[srt]-e_sig_Hbeta[srt], sig_Hbeta[srt]+e_sig_Hbeta[srt], color='C2', alpha=0.1)
plt.plot(rsec[srt], sig_Hbeta[srt], marker='.', linestyle='-', label=r'$H\beta$', color='C2')
plt.fill_between(rsec[srt], sig_stars[srt]-e_sig_stars[srt], sig_stars[srt]+e_sig_stars[srt], color='k', alpha=0.1)
plt.plot(rsec[srt], sig_stars[srt], marker='.', linestyle='-', label=r'Stars', color='k')
plt.ylim(0, 200)
plt.xlim(xlim)
plt.axvline(0.0, linestyle=':')
plt.xlabel(xlabel)
plt.title('$\sigma_\mathrm{LOS}$, km/s')
plt.figure(figsize=(16, 4))
ax = plt.subplot(131)
plt.plot(rsec, age_stars, 'k.')
plt.axvline(0.0, linestyle=':')
plt.ylim(3, 14)
plt.yscale('log')
plt.xlim(xlim)
plt.xlabel(xlabel)
plt.title('$<Age>_\mathrm{mass}$, Gyr')
ax = plt.subplot(132)
plt.plot(rsec, met_stars, 'k.')
plt.axvline(0.0, linestyle=':')
plt.ylim(-1.0, 0.5)
plt.xlim(xlim)
plt.xlabel(xlabel)
plt.title('$<[Fe/H]>_\mathrm{mass}$, dex')
ax = plt.subplot(133)
plt.plot(rsec, ml_stars, 'k.')
plt.axvline(0.0, linestyle=':')
plt.ylim(0, 8)
plt.xlim(xlim)
plt.xlabel(xlabel)
plt.title('$M/L_\mathrm{r-band}$')
plt.show()
2.3 Заключительные замечания¶
В работе показан полный цикл анализа данных от первичной редукции до моделирования спектра. Отметим слабые места представленной реализации анализа:
- Инструментальный контур Мы предполагаем, что разрешение спектра SCORPIO очень близко к разрешению моделей MILES. Для точных измерений, в первую очередь, дисперсии скоростей - требуется точное измерение инструментального разрешения, например, по анализу эмиссий в кадрах неона или по анализу спектра рассветного неба (sunsky). Внимательно рассматривая линии неона можно заметить, что к краю кадра линии не просто меняют свою ширину (увеличиваются), но также меняют свою форму - у линий появляется "голубое крыло". Для особо точных измерений форму контура (LSF - line spread function) можно восстанавливать из спектра неона или sunsky и учитывать их при моделировании спектра. Это достаточно сложная процедура, особенно в случае использования моделей населений сопоставимого разрешения. Но существуют модели звездных населений более высокого разрешения, где для учета LSF достаточно свернуть эти модели с восстановленным контуром LSF перед попиксельным моделированием спектра.
- Калибровка спектральной чувствительности Как правило каждую наблюдательную ночь снимается спектро-фотометрический стандарт. Его следует анализировать как и кадры объекта, а после извлечь одномерный спектр и сопоставив с реперным спектром данного стандарта определить кривую спектральной чувствительности. Каждый спектр нужно делить на эту кривую, чтобы спектр получался правильно формы. В фите с помощью
pPXFиспользуется мультипликативным континуум, которые нивелирует расхождения между формой спектра моделей зв. населений и спектра галактики. Однако, чем меньше это расхождение, тем стабильнее будет фит. - Флэт Калибровка флэта - это тоже спектр калибровочной лампы. Поэтому его лучше сначала линеризовать, а только потом использовать для расчета коррекции плоского поля.
- Биннинг Сейчас мы анализируем звездный континуум и эмиссии использую одинаковую схему бинирования (binNum). Это не совсем корректно, потому что эмиссии и звезды имеют существенно разное распределения яркости вдоль щели. Поэтому правильная стратегия - это после первого прогона, рассчитать чисто эмиссионный спектр вычитая из наблюдаемого спектра спектр звездного континуума (в каждом небинированном спектра), а после этого рассчитать схему бинирования (адаптивно или вручную) используя сигнал-шум в эмиссионных линиях. После этого провести анализ эмиссионного спектра.