234 Глава 8. Обнаружение далеких экзопланет
Я
В а а
К а а
Г а
В 
Рис. 8.3. Глубина представляет полное изменение яркости, наблюдаемое в кривой блеска
Проект #11. Симуляция транзита экзопланеты
Прежде чем полететь в Айдахо фотографировать Великое американское затмение в 2017 году, я подготовился. Было известно, что полная фаза, в течение которой Луна будет полностью закрывать Солнце, продлится всего 2 минуты 10 секунд. Так что времени на эксперименты, тесты и выяснение деталей на ходу не будет. Чтобы успешно заснять конус полутени, тень, протуберанцы и эффект бриллиантового кольца (рис. 8.4), мне нужно было знать, какое именно брать оборудование, какие настройки камеры использовать и когда эти явления должны произойти.
Схожим образом компьютерные симуляции помогают нам подготовиться к наблюдениям природных явлений. Они помогают понять, чего ожидать, когда ожидать и как откалибровать инструменты. В текущем проекте мы создадим симуляцию транзита экзопланеты. Эту симуляцию вы сможете выполнять с планетами различных размеров, чтобы понять влияние размера на кривую блеска. Позже с помощью данной симуляции мы оценим кривые блеска, относящиеся к полям астероидов и возможным инопланетным мега структурам.
ЗАДАЧА
Написать программу Python, которая будет симулировать транзит экзопланеты, отображать получающуюся кривую блеска и вычислять радиус этой планеты.
Проект #11. Симуляция транзита экзопланеты 235
Рис. 8.4. Эффект бриллиантового кольца в конце фазы полного солнечного затмения в 2017 году
Стратегия
Для генерации кривой блеска потребуется возможность измерять изменения яркости. Обратимся к OpenCV и выполним математические операции над пикселями, такими как поиск среднего, минимального и максимального значений.
Вместо использования изображения реального транзита и звезды мы нарисуем круги на черном прямоугольнике аналогично тому, как рисовали прямоугольники на карте Марса в предыдущей главе. Для отражения графика кривой блеска можно использовать matplotlib, основную библиотеку Python для создания графиков. Ее мы установили в разделе «Установка NumPy и других научных библиотек с помощью pip» на с. 33 и уже использовали для создания графиков в главе 2.
Код для транзита
Программа transit.py с помощью OpenCV генерирует визуальную симуляцию экзопланеты, проходящей на фоне звезды, с помощью matplotlib создает
236 Глава 8. Обнаружение далеких экзопланет
график итоговой кривой блеска и в конце оценивает размер планеты, используя уравнение планетарного радиуса со с. 233. Код можете ввести сами или скачать с https://nostarch.com/real-world-python/.
Импорт модулей и присваивание констант
В листинге 8.1 мы импортируем модули и присваиваем константы, представляющие вводимые пользователем данные.
Листинг 8.1. Импорт модулей и присваивание констант transit.py, часть 1
import math
import numpy as np import cv2 as cv
import matplotlib.pyplot as plt
IMG_HT = 400
IMG_WIDTH = 500
BLACK_IMG = np.zeros((IMG_HT, IMG_WIDTH, 1), dtype='uint8')
STAR_RADIUS = 165
EXO_RADIUS = 7
EXO_DX = 3
EXO_START_X = 40
EXO_START_Y = 230
NUM_FRAMES = 145
Импортируем модуль math для уравнения радиуса планеты, NumPy для вычисления яркости изображения, OpenCV для рисования симуляции и matplotlib для графического отображения кривой блеска. Затем присваиваем константы, которые представляют пользовательские значения.
Начинаем с высоты и ширины окна симуляции. Это окно будет черным прямоугольным изображением. Его создаем с помощью метода np.zeros(), который возвращает заполненный нулями массив установленной формы и типа.
Напомню, что изображения OpenCV — это массивы NumPy, элементы которых должны иметь одинаковый тип. Тип данных uint8 представляет беззнаковое целое между 0 и 255. Полезный список других типов данных и их описания можете найти на странице https://numpy.org/devdocs/user/basics.types.html.
Далее присваиваем значения радиуса в пикселях для звезды и экзопланеты. OpenCV использует эти константы, когда рисует изображающие их круги.
Экзопланета будет двигаться перед звездой, поэтому нужно задать скорость ее движения. Константа EXO_DX инкрементирует позицию x планеты на три пикселя с каждым циклом программы, в результате чего планета будет смещаться слева направо.
Проект #11. Симуляция транзита экзопланеты 237
Присваиваем две константы, устанавливая начальную позицию планеты. Затем присваиваем константу NUM_FRAMES для управления количеством обновлений симуляции. Несмотря на то что это число можно вычислить (IMG_WIDTH/EXO_DX), его присваивание позволяет точно подстроить продолжительность симуляции.
Определение функции main()
Код листинга 8.2 определяет функцию main() для выполнения программы. Можно определить main() в любом месте, но размещение ее в начале позволит ей предоставлять контекст для функций программы, определяемых далее. В составе main() можно также вычислить радиус экзопланеты, вложив соответствующее уравнение в вызов функции print().
Листинг 8.2. Определение функции main() transit.py, часть 2
def main():
intensity_samples = record_transit(EXO_START_X, EXO_START_Y) relative_brightness = calc_rel_brightness(intensity_samples) print('\nestimated exoplanet radius = {:.2f}\n'
.format(STAR_RADIUS * math.sqrt(max(relative_brightness)
- min(relative_brightness))))
plot_light_curve(relative_brightness)
После определения main() создаем переменную intensity_samples и вызываем функцию record_transit(). Интенсивность (intensity) означает количество света, представленное численным значением пикселя. Функция record_transit() отрисовывает симуляцию на экране, измеряет ее интенсивность, добавляет полученные измерения в список intensity_samples, после чего этот список возвращает. Данной функции требуется начальная точка координат (x, y) экзопланеты. Передаем ей начальные константы EXO_START_X и EXO_START_Y, помещая тем самым планету в позицию, аналогичную на рис. 8.2. Обратите внимание, если существенно увеличить радиус экзопланеты, то может потребоваться сдвинуть начальную точку левее (отрицательные значения допускаются).
Далее создаем переменную relative_brightness и вызываем функцию calc_rel_ brightness(). Эта функция вычисляет относительную яркость, которая равна измеренной интенсивности, поделенной на ее максимальное зарегистрированное значение. Она получает в качестве аргумента список измерений интенсивности, конвертирует их в относительную яркость и возвращает новый список.
Этот список значений относительной яркости мы используем для вычисления радиуса экзопланеты в пикселях с помощью уравнения со с. 233. Можно выполнить эти вычисления в рамках функции print(). Для сообщения ответа с точностью до двух десятичных знаков используем {:.2f}.
238 Глава 8. Обнаружение далеких экзопланет
Завершается main() вызовом функции для рисования графика кривой блеска, которой передается список relative_brightness.
Регистрация транзита
Код листинга 8.3 определяет функцию для симуляции и регистрирования транзита. Он рисует звезду и экзопланету на изображении черного прямоугольника, после чего выполняет смещение экзопланеты. Здесь также происходит вычисление и отображение средней интенсивности изображения с каждым движением, добавление этого значения интенсивности в список и возвращение итогового списка в конце.
Листинг 8.3. Рисование симуляции, вычисление интенсивности изображения и возвращение ее значений в виде списка
transit.py, часть 3
def record_transit(exo_x, exo_y):
"""Нарисовать планету, проходящую мимо звезды, и вернуть список изменений интенсивности"""
intensity_samples = []
for _ in range(NUM_FRAMES): temp_img = BLACK_IMG.copy()
cv.circle(temp_img, (int(IMG_WIDTH / 2), int(IMG_HT / 2)), STAR_RADIUS, 255, -1)
cv.circle(temp_img, (exo_x, exo_y), EXO_RADIUS, 0, -1) intensity = temp_img.mean()
cv.putText(temp_img, 'Mean Intensity = {}'.format(intensity), (5, 390), cv.FONT_HERSHEY_PLAIN, 1, 255)
cv.imshow('Transit', temp_img) cv.waitKey(30)
intensity_samples.append(intensity) exo_x += EXO_DX
return intensity_samples
Функция record_transit() получает в качестве аргументов пару координат (x, y). Они показывают начальную позицию экзопланеты или, говоря точнее, пиксель, используемый в качестве центра для первого нарисованного круга. При этом он не должен пересекаться с кругом, представляющим звезду и размещенным по центру изображения.
Далее создаем пустой список для хранения измерений интенсивности. После этого начинаем цикл for, который использует константу NUM_FRAMES для повторения симуляции определенное число раз. Симуляция должна длиться немного дольше, чем требуется для полного выхода экзопланеты за пределы диска звезды. Таким образом, мы получим кривую блеска, включающую измерения после транзита.
Картинки и текст, которые мы размещаем на изображении при помощи OpenCV, становятся частью этого изображения. Следовательно, с каждым циклом нужно
Проект #11. Симуляция транзита экзопланеты 239
заменять предыдущее изображение, копируя его оригинал BLACK_IMG в локальную переменную temp_img.
Теперь можно рисовать звезду, используя метод OpenCV circle(). Передаем ему временное изображение, координаты (x, y) для центра круга, соответствующего центру изображения, константу STAR_RADIUS, белый цвет для заполнения и толщину линии. Используя отрицательное значение для толщины, мы заполняем цветом весь круг.
Далее рисуем экзопланету. Берем координаты exo_x и exo_y в качестве начальной точки, константу EXO_RADIUS в качестве размера и черный цвет для заполнения .
В этот момент нужно начинать регистрировать интенсивность изображения. Поскольку пиксели уже представляют интенсивность, нужно просто взять среднее по изображению. Количество проводимых измерений зависит от константы EXO_DX. Чем больше ее значение, тем быстрее движется экзопланета и тем реже будет зарегистрирована средняя интенсивность.
Для отображения вычисленной интенсивности на рисунке используем метод OpenCV putText(). Ему передаются временное изображение, текстовая строка, включающая сами измерения, координаты (x, y) для нижнего левого угла текстовой строки, шрифт, его размер и цвет.
Теперь называем окно Transit и отображаем его с помощью метода OpenCV imshow(). На рис. 8.5 показана итерация цикла.
Рис. 8.5. Экзопланета, проходящая на фоне звезды
240 Глава 8. Обнаружение далеких экзопланет
После вывода изображения используем метод OpenCV waitKey() для его обновления каждые 30 мс. Чем меньше переданное этому методу значение, тем быстрее экзопланета будет перемещаться.
Добавляем измеренное значение интенсивности в список intensity_samples, после чего продвигаем круг экзопланеты, увеличивая ее значение exo_x на константу EXO_DX . Завершается функция возвращением списка значений средней интенсивности.
Вычисление относительной яркости и построение кривой блеска
Код листинга 8.4 определяет функцию для вычисления относительной яркости каждого значения интенсивности и построения графика кривой блеска. Затем следует вызов функции main(), если программа не используется в качестве модуля в другой программе.
Листинг 8.4. Вычисление относительной яркости, отображение графика кривой блеска и вызов main()
transit.py, часть 4
def calc_rel_brightness(intensity_samples):
"""Вернуть список относительной яркости из списка значений интенсивности"""
rel_brightness = []
max_brightness = max(intensity_samples) for intensity in intensity_samples:
rel_brightness.append(intensity / max_brightness) return rel_brightness
def plot_light_curve(rel_brightness):
"""Построить график изменения относительной яркости в зависимости от времени"""
plt.plot(rel_brightness, color='red', linestyle='dashed', linewidth=2, label='Relative Brightness')
plt.legend(loc='upper center') plt.title('Relative Brightness vs. Time') plt.show()
if __name__ == '__main__': main()
Кривые блеска отражают относительную яркость с течением времени так, что полностью незатемненная звезда имеет значение 1.0, а полностью затемненная — значение 0.0. Чтобы преобразовать измерения средней интенсивности в относительные значения, мы определяем функцию calc_rel_brightness(), которая получает в качестве аргумента список усредненных измерений интенсивности.
Проект #11. Симуляция транзита экзопланеты 241
Внутри этой функции создаем пустой список для хранения конвертированных значений, после чего используем встроенную в Python функцию max() для поиска максимального значения в списке intensity_samples(). Для получения относительной яркости перебираем элементы этого списка и делим их на максимальное значение. По ходу дела добавляем результаты в список rel_brightness. Завершаем функцию возвращением нового списка.
Определяем вторую функцию для нанесения графика кривой блеска и передаем ей список rel_brightness . Используем метод plot() библиотеки matplotlib(), передавая ему этот список, цвет линии, стиль линии, ее толщину и метку для легенды графика. Добавляем легенду и название графика, после чего его отображаем. Результат — на рис. 8.6.
Рис. 8.6. Пример графика кривой блеска из transit.py
На первый взгляд отклонение яркости на этом графике может показаться чрезвычайно большим, но если поближе рассмотреть ось y, то можно видеть, что экзопланета уменьшила яркость звезды всего на 0.175 процента. Чтобы понять,