[천문 + 프로그래밍] Missing Star 생성기

8월 10일부터 20일까지 폴란드에서 진행된 2023 IOAA(국제 천문 및 천체물리 올림피아드)에 참가하였는데요, 이 대회를 준비하려고 쓴 Missing Star 문제를 생성해주는 프로그램을 만들었었습니다. 이 글에서는 그 프로그램을 소개하려고 합니다.

---

Missing Star

Missing Star는 성도에서 랜덤한 별이 몇 개가 빠지고, 그 별이 어떤 별인지, 원래 어디 있어야 하는지를 찾는 문제입니다. Missing Star 문제를 풀려면 밤하늘에서 어떤 별이 어디에 있는지를 알아야 하기 때문에, 성도를 잘 암기해야 하죠.

예를 들면

이 성도에서는 아래 별들이 빠졌습니다. 빠진 별들을 표시해보면 아래와 같습니다.

$\alpha$ Aql
$\alpha$ Aur
$\alpha$ Gem

---

프로그램 소개

1. Setting

우선, 필요한 패키지를 모두 불러오고,

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from matplotlib.backends.backend_pdf import PdfPages
import random

히파르코스 위성의 데이터를 받아와서 별들의 적경 적위 값 데이터를 얻고, 계산을 위해 도(degree)를 라디안(radian)으로 변환하습니다.

pdf = PdfPages('Star Chart.pdf')

df = pd.read_csv("히파르코스 위성 데이터 파일 경로")

df = pd.DataFrame(df, columns=['RAJ2000','DEJ2000','Vmag'])
ra = df['RAJ2000'].to_numpy()
dec = df['DEJ2000'].to_numpy()
mag = df['Vmag'].to_numpy()
ra = ra*np.pi/180
dec = dec*np.pi/180

hip1.csv

3.43MB

다음으로, 랜덤한 위치와 시간의 밤하늘을 얻기 위해, 항성시(남중한 천체의 적경)와 위도를 np.random.rand 함수로 무작위 값을 얻어 설정하였습니다. 또한 히파르코스 위성 데이터에는 눈에 보이지 않는 어두운 별의 데이터 또한 포함되어 있으므로 성도 한계 등급을 설정해줍니다.

lat = np.pi*np.random.rand(1) - 0.5*np.pi
lst = 2*np.pi*np.random.rand(1)
mag_lim = 5.5 #성도 한계 등급

 

2. 함수 정의

어떤 별의 적경/적위 값을 알고, 항성시와 위도를 안다면 그 천체의 방위각과 고도 또한 구면삼각법을 이용하여 계산할 수 있습니다. 적경($RA$)과 항성시($LST$) 사이 관계를 이용하면 시간각($HA$)은

$$HA=LST-RA$$

이고, 이를 통해 계산한 고도($h$)는

$$\sin h=\sin\phi\sin DEC+\cos\phi\cos DEC\cos HA$$

여기서 $DEC, \phi$는 적위와 위도입니다.

또한 방위각($A$)은 위 값들과 아래와 같은 관계를 가지고 있습니다.

$$\sin DEC=\sin\phi\sin h+\cos\phi\cos h\cos A$$

$$\sin A\cos h=-\cos DEC\sin HA$$

이를 통해 $A$값을 결정할 수 있습니다. 고도와 방위각을 계산했으니, 성도에 그 좌표대로 별을 찍으면 우리가 설정한 랜덤한 밤하늘에서의 성도를 얻게 되는 것입니다.

이 방법으로 고도와 방위각을 구하는 함수를 정의하여 사용합시다.

def a_func(raf, decf, latf, lstf):
    haf = lstf - raf
    sin_af = np.sin(latf)*np.sin(decf) + np.cos(latf)*np.cos(decf)*np.cos(haf)
    return np.arcsin(sin_af)

def az_func(af, raf, decf, latf, lstf):
    haf = lstf - raf
    cos_az = (np.sin(decf) - np.sin(latf)*np.sin(af))/(np.cos(latf)*np.cos(af))
    sin_az = -np.cos(decf)*np.sin(haf)/np.cos(af)
    az_0 = np.arcsin(sin_az)
    if sin_az >= 0 and cos_az >= 0:
        az = az_0
        k = 0
    elif sin_az >= 0 and cos_az < 0:
        az = np.pi - az_0
        k = 1
    elif sin_az < 0 and cos_az < 0:
        az = -np.pi - az_0
        k = 2
    elif sin_az < 0 and cos_az >=0:
        az = az_0
        k = 3
    return az

다음으로 성도를 그릴 plot을 설정하는 함수를 정의합니다. 극좌표 plot을 사용하였고, 눈금을 모두 없애주었습니다.

def plotsetting():
    ax = fig.add_subplot(projection='polar')
    plt.grid(False)
    ax.set_ylim(0, 1)
    title = 'Star Chart (latitude :' + str(lat*180/np.pi) + 'degree , LST :' + str(lst*12/np.pi) + 'hour)'
    ax.set_title(title , va='bottom')
    ax.set_yticks(np.arange(0,1,1), labels=[''])
    ax.set_xticks(np.arange(0,2.0*np.pi,np.pi/12.0), labels=['', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', ''])

 

3. 성도 생성

방위도 무작위 방향으로 돌려놓기 위해 방위각을 np.random.rand 함수로 랜덤 값을 뽑아주고, 한계 등급보다 밝은 별들의 데이터만 뽑아 리스트에 저장해줍니다.

plot 크기를 설정하고, 위에서 정의한 plotsetting 함수로 plot을 설정해줍니다.

az_random = 2*np.pi*np.random.rand(1)
star_num = [x for x in range(len(mag)) if mag[x] < mag_lim and a_func(ra[x], dec[x], lat, lst) >= 0]

fig = plt.figure(figsize=(20, 20))
plotsetting()

이제 성도에서 뺄 별을 설정해야 합니다. missing star의 개수와 한계 등급을 설정하고, 아까 뽑은 별의 리스트에서 missing star 조건을 만족하고, 고도가 10도 이상인 별들의 리스트를 만들어줍니다. 그리고 만든 리스트를 무작위로 섞어줍니다.

missing_cnt = 3 #미싱스타 개수
missing_lim = 3.5 #미싱스타 한계 등급
missing_star_num = [x for x in range(len(mag)) if mag[x] < missing_lim and a_func(ra[x], dec[x], lat, lst) >= 10*np.pi/180]
random.shuffle(missing_star_num)

별들을 찍을 차례입니다. 아까 정의한 함수들을 이용하여 고도와 방위각을 계산해주고, 평사도법을 사용하여 별을 극좌표 plot에 찍어줍니다. 이때 별의 크기는 실제 광도와 비례하게 등급에 따라 지수함수적으로 달라지게 설정해주었습니다. 이때 아까 만든 missing star 리스트에서 맨앞부터 missing star의 개수만큼의 별들은 표시하지 않습니다. 별들을 모두 표시하면 plot을 이미지로 저장합니다.

for i in star_num:
	a_val = a_func(ra[i], dec[i], lat, lst)
	az_val = az_func(a_val, ra[i], dec[i], lat, lst) + az_random
	z_val = np.pi/2 - a_val
	l_val = np.sin(z_val)/(1+np.cos(z_val)) 
	star_size = 10**(0.40*(6-mag[i])) - 0.7
	if not(i in missing_star_num[:missing_cnt]):
    	plt.scatter(az_val, l_val, s = star_size, color='k')

pdf.savefig(fig)

문제를 만들었으니 답지도 만들어야겠죠. 똑같은 방식으로 별을 찍는데, 이번에는 아까 표시하지 않았던 별들은 빨간색으로 표시해줍니다.

fig = plt.figure(figsize=(20, 20))
plotsetting()
for i in star_num:
    a_val = a_func(ra[i], dec[i], lat, lst)
	az_val = az_func(a_val, ra[i], dec[i], lat, lst) + az_random
    z_val = np.pi/2 - a_val
    l_val = np.sin(z_val)/(1+np.cos(z_val))
    star_size = 10**(0.40*(6-mag[i])) - 0.7
    if i in missing_star_num[:missing_cnt]:
        plt.scatter(az_val, l_val, s = star_size, color='r')    
    else:
        plt.scatter(az_val, l_val, s = star_size, color='k')

pdf.savefig(fig)

다음 문제를 생성하기 위해 위도와 항성시를 다시 랜덤 설정해주고, 이전에 그린 성도 plot을 닫아줍니다.

lat = np.pi*np.random.rand(1) - 0.5*np.pi
lst = 2*np.pi*np.random.rand(1)
plt.close()

이 3번 과정을 for문에 넣어 원하는 개수만큼 성도를 생성해줍니다. 모두 생성했으면 생성한 성도가 들어있는 pdf를 닫고 저장합니다.

 

따라서 전체 코드는 아래와 같습니다.

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from matplotlib.backends.backend_pdf import PdfPages
import random

pdf = PdfPages('Star Chart.pdf')

df = pd.read_csv("히파르코스 위성 데이터 파일 경로")

df = pd.DataFrame(df, columns=['RAJ2000','DEJ2000','Vmag'])
ra = df['RAJ2000'].to_numpy()
dec = df['DEJ2000'].to_numpy()
mag = df['Vmag'].to_numpy()
ra = ra*np.pi/180
dec = dec*np.pi/180

lat = np.pi*np.random.rand(1) - 0.5*np.pi
lst = 2*np.pi*np.random.rand(1)
mag_lim = 5.5 #성도 한계 등급

def a_func(raf, decf, latf, lstf):
    haf = lstf - raf
    sin_af = np.sin(latf)*np.sin(decf) + np.cos(latf)*np.cos(decf)*np.cos(haf)
    return np.arcsin(sin_af)

def az_func(af, raf, decf, latf, lstf):
    haf = lstf - raf
    cos_az = (np.sin(decf) - np.sin(latf)*np.sin(af))/(np.cos(latf)*np.cos(af))
    sin_az = -np.cos(decf)*np.sin(haf)/np.cos(af)
    az_0 = np.arcsin(sin_az)
    if sin_az >= 0 and cos_az >= 0:
        az = az_0
        k = 0
    elif sin_az >= 0 and cos_az < 0:
        az = np.pi - az_0
        k = 1
    elif sin_az < 0 and cos_az < 0:
        az = -np.pi - az_0
        k = 2
    elif sin_az < 0 and cos_az >=0:
        az = az_0
        k = 3
    return az    

def plotsetting():
    ax = fig.add_subplot(projection='polar')
    plt.grid(False)
    ax.set_ylim(0, 1)
    title = 'Star Chart (latitude :' + str(lat*180/np.pi) + 'degree , LST :' + str(lst*12/np.pi) + 'hour)'
    ax.set_title(title , va='bottom')
    ax.set_yticks(np.arange(0,1,1), labels=[''])
    ax.set_xticks(np.arange(0,2.0*np.pi,np.pi/12.0), labels=['', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', '', ''])

map_num = 20

for j in range(map_num):

    az_random = 2*np.pi*np.random.rand(1)
    star_num = [x for x in range(len(mag)) if mag[x] < mag_lim and a_func(ra[x], dec[x], lat, lst) >= 0]

    #missing star
    fig = plt.figure(figsize=(20, 20))
    plotsetting()

    missing_cnt = 3 #미싱스타 개수
    missing_lim = 3.5 #미싱스타 한계 등급
    missing_star_num = [x for x in range(len(mag)) if mag[x] < missing_lim and a_func(ra[x], dec[x], lat, lst) >= 10*np.pi/180]
    random.shuffle(missing_star_num)

    for i in star_num:
        a_val = a_func(ra[i], dec[i], lat, lst)
        az_val = az_func(a_val, ra[i], dec[i], lat, lst) + az_random
        z_val = np.pi/2 - a_val
        l_val = np.sin(z_val)/(1+np.cos(z_val))
        star_size = 10**(0.40*(6-mag[i])) - 0.7
        if not(i in missing_star_num[:missing_cnt]):
            plt.scatter(az_val, l_val, s = star_size, color='k')

    pdf.savefig(fig)

    #missing star solution
    fig = plt.figure(figsize=(20, 20))
    plotsetting()
    for i in star_num:
        a_val = a_func(ra[i], dec[i], lat, lst)
        az_val = az_func(a_val, ra[i], dec[i], lat, lst) + az_random
        z_val = np.pi/2 - a_val
        l_val = np.sin(z_val)/(1+np.cos(z_val))
        star_size = 10**(0.40*(6-mag[i])) - 0.7
        if i in missing_star_num[:missing_cnt]:
            plt.scatter(az_val, l_val, s = star_size, color='r')    
        else:
            plt.scatter(az_val, l_val, s = star_size, color='k')

    pdf.savefig(fig)

    lat = np.pi*np.random.rand(1) - 0.5*np.pi
    lst = 2*np.pi*np.random.rand(1)
    plt.close()

pdf.close()

---

성도 좋아요

Missing Star 시간 때우기 정말 좋습니다