2021.05.18 - [과학/천문학] - [우주과학과 학생들은 어떤 과제를 하나요?] 국제우주정거장 (ISS)을 스마트폰을 통해 관측하기
밸런스 있는 삶을 찾아주신 여러분들 모두 반갑습니다 :)
오늘은 비행/궤도 역학에 관련된 개념을 파이썬의 함수(def)를 통해 구현하는 과정을 보일 것입니다!!
저번 게시글에서는 스마트폰 어플리케이션을 통하여 국제우주정거장 관측을 진행하였습니다. 그리고 저번 게스글 마지막에 예고했듯이 이번 시간에는 ISS Detector 어플리케이션에서 제공해주는 국제우주정거장의 거리, 적경, 적위 값과 위성 궤도 계산 과정을 파이썬에서 함수에서 구현하여, 국제우주정거장의 궤도에 대한 정보를 구하는 과정을 보일 것입니다. 다만, 내용이 너무 길어 두 번의 게시글에 걸쳐서 보여드리겠습니다!
스마트폰 어플리케이션을 통하여 국제우주정거장을 관측하는 과정이 궁금하시다면 맨 위에 첨부되어있는 링크를 통해 ISS Detector 어플로 국제우주정거장 관측 방법에 대해 참고하시면 되겠습니다.
다시 본론으로 돌아와서 이번 게시글에서는 ISS Detector에서 제공하는 국제우주정거장의 적경, 적위, 거리 값만을 이용하여 주어진 값들을 파이썬에서 함수(def)로 제가 직접 새로운 함수를 구현할 것입니다.
제가 이번에 구현할 함수는 총 4가지 입니다.
1) ISS Detctor에서 제공한 적경의 시간단위를 도(degree)단위로 바꾸는 함수
2) ISS Detctor에서 제공한 적경과 적위, 그리고 거리를 통해 국제우주정거장의 적도좌표계 위치를 나타내는 함수
3) 2개의 시간대에서 구해낸 적도좌표계 위치들에서 좌표간 거리 공식을 구현하여 특정 시간동안 국제우주정거장의 이동거리를 구하는 함수
4) '국제우주정거장의 이동거리 / 관측 시간' 을 통해 국제우주정거장의 공전(비행)속도를 구하는 함수
이렇게 4가지 함수를 통해 ISS Detctor에서 제공한 국제우주정거장의 적경과 적위, 그리고 거리라는 한정적인 정보를 파이썬에 함수(def)로 구현하여 국제우주정거장의 적도좌표계 위치, 이동거리, 그리고 속력이라는 국제우주정거장의 궤도에 대한 정보를 도출해낼 것입니다. (아래 코드블럭의 1번~4번 과정)
아래는 제가 두 번의 게시물을 거쳐 소개해드릴 ISS Detector 어플리케이션에서 제공해주는 국제우주정거장의 거리, 적경, 적위 값을 통해 국제우주정거장의 적도좌표계 위치, 공전(비행) 속도, 궤도 장반경, 단위질량당 에너지, 단위질량당 각운동량, 이심률이라는 국제우주정거장 궤도에 대한 최대한의 정보를 구하는 과정입니다.
# 1 시간단위의 적경을 도(degree)단위의 각도로 고치기
# 2 두 시간대의 적경과 적위, 거리값을 적도좌표계인 x,y,z 3차원의 좌표로 출력한다.
# 3 그 두 좌표 사이의 거리를 구한다.
# 4 그러면 속도를 알 수 있다.
# 5 거리는 고정값으므로 장반경을 알 수 있다.
# 6 뮤와 장반경을 알기 때문에 에너지도 알 수 있다.
# 7 속도 벡터와 거리 벡터를 알고 있으므로 h를 구할 수 있다.
# 8 h를 알기에 이심률 또한 알 수 있다.(근/원일점 거리도)
그리고 다음은 제가 이번 과정동안 이용할 각종 상수들입니다. (출처 : 위키백과)
# ISS 기본 정보
m = 419725 # iss 질량
M = 5.972e24 # 지구 질량
MU = 398600 # 지구 질량 * 중력상수
r_p = 331 + 6378 # 근지점
r_a = 410 + 6378 # 원지점
i_angle = 51.6410 # 궤도 경사각
velocity = 7.7066 # 속도
period = 91.68 * 60 # 주기
그리고 이 상수들과 지구를 중심으로한 궤도역학적 개념을 사용하여 제가 구하고자 한 공전(비행) 속도, 궤도 장반경, 단위질량당 에너지, 단위질량당 각운동량, 이심률의 이론값을 구해줄 것입니다.
우선 사인, 코사인, 제곱근 등 사용해야 할 수학적 개념이 많으므로 넘파이를 import 해줍니다.
import numpy as np
그리고 적경과 적위, 그리고 거리를 유도하기 위한 국제우주정거장의 고도값은 ISS Detector 어플리케이션을 통해 다음 사진과 같은 위치에서 구할 수 있습니다.
이제 적경, 적위 그리고 고도를 찾으셨나요?
그럼 본격적으로 궤도역학적 개념을 함수로 표현하는 과정을 시작합니다!
저는 1번의 "시차단위의 적경을 도(degree)단위의 각도로 고치기" 과정을 구현하기 위해 아래의 세 단계를 거쳤습니다.
- 시, 분, 초로 나타나 있는 적경을 list형태로 나타내기
- list 안의 값을 시, 분, 초 별로 1, 60, 3600을 곱하여 도의 단위로 수정하기
- 수정된 도 단위의 적경과 마찬가지로 도 단위인 적위를 라디안 단위로 고쳐주기
시,분,초를 동시에 담는 적경 값을 한 번에 계산하기 위해 저는 적경을 인덱싱이 가능한 list의 형태로 두었습니다.
그 후 적경이 8h 42' 51''로 표현되어 있으면 42'에는 60'/1h을 나누고, 51''에는 3600''/1h를 나누어 8h 42' 51''를 모두 시(h)에 대한 값으로 고쳐주었습니다.
(∵ 1'' = 1º를 3600등분 한 것)
그리고 도(degree)단위의 각도는 시차단위의 각도의 15배이므로 이를 표현하는 deg_angle = (angle1 + angle2 + angle3) * 15라는 코드를 작성하여 시차 단위의 적경을 모두 도 단위의 적경으로 바꾸어주었습니다.
(∵ 24h가 지났을 때 한바퀴, 360º 돌았다고 표현한다. 그러므로 15배의 관계 성립)
ISS Detector에서 적위값은 도 단위로 제공하기에 이 값은 그대로 변수로 지정하였습니다.
※ Right Ascension 적경(= RA), Declination 적위(= DEC) ※
변수들을 지정한 후 def과정을 통해 다음과 같이 코드를 작성하고 print를 해주면 다음과 같은 결과가 나옵니다.
# 1 시간단위의 적경을 도(degree)단위의 각도로 고치기
RA_0 = [8,42,51] # [시,분,초] (처음 관측 시간 때 정보) , 적경을 나타냄
DEC_0 = 37.1 # '도' 단위 (처음 관측 시간 때 정보) , 적위을 나타냄
RA_t = [8,52,24] # [시,분,초] (나중 관측 시간 때 정보)
DEC_t = 35.5 # '도' 단위 (나중 관측 시간 때 정보)
def RA_time_to_deg(RA): # 시차에 해당하는 적경을 라디안으로 수정하는 함수
angle1 = RA[0] # h 적경의 시간에 해당
angle2 = RA[1] / 60. # min 적경의 분에 해당
angle3 = RA[2] / 3600. # sec 적경의 초에 해당
deg_angle = (angle1 + angle2 + angle3) * 15 # 시차 단위의 적경을 각도로 변경
radian_angle = deg_angle * (np.pi/180) # '도'로 표현된 각을 라디안으로 수정
return radian_angle
print('처음 관측 할 때 적경 = ',RA_time_to_deg(RA_0), 'rad')
print('처음 관측 할 때 적위 = ',DEC_0 * (np.pi/180), 'rad')
print('나중 관측 할 때 적경 = ', RA_time_to_deg(RA_t), 'rad')
print('처음 관측 할 때 적위 = ',DEC_t * (np.pi/180), 'rad')
참고로, 저는 함수명이나 변수들을 지정해줄 때 남들도 알아보기 쉽게 조금은 길지만 정확한 단어로 표기하려고 합니다. 여러분은 저처럼 이렇게 긴 단어를 사용할 필요는 없지만, 적어도 누군가가 봤을 때 어떤 의미의 코드이고 어떤 의미의 변수인지 알아볼 수 있게 작성하시면 더 좋을 것 같습니다 :)
그 후 지구의 반경과 어플리케이션에 나타난 고도를 더하여 거리를 구하였습니다. 그 후 적경, 적위, 거리를 통해 임의의 두 지점에서의 국제우주정거장의 적도좌표계 상 좌표로 표현하였습니다. 적도좌표계에서 거리, 적경, 적위로 위치(좌표)를 표현하는 방법은 다음과 같습니다.
위의 그림에서 빨간 상자 안에 있는 내용을 코드로 구현하여 1번째 항을 x에 대한 요소, 2번째 항을 y에 대한 요소, 그리고 3번째 항을 z에 대한 요소로 표현하여 각각을 list 형식에 담아주어 이를 함수로 정의할 것입니다.
코드와 그 출력값은 다음과 같습니다.
# 2 두 시간대의 적경과 적위, 거리값을 적도좌표계인 x,y,z 3차원의 좌표로 출력한다.
earth_radius = 6378 # 지구 반경
R_0 = 422 + earth_radius # km 단위 (어플에서 주는 처음고도 값 + 지구반경)
R_t = 421 + earth_radius # km 단위 (어플에서 주는 나중고도 값 + 지구반경)
def Equatorial_Coordinate_System(RA, DEC, R): # 적경,적위,거리를 x,y,z평면의 죄표로 표현하는 함수
x = R * np.cos(DEC* (np.pi/180)) * np.cos(RA_time_to_deg(RA))
y = R * np.cos(DEC* (np.pi/180)) * np.sin(RA_time_to_deg(RA))
z = R * np.sin(DEC* (np.pi/180))
location_list = [x,y,z]
return location_list # km
print('처음 관측할 때의 좌표 (x,y,z) = ', Equatorial_Coordinate_System(RA_0, DEC_0, R_0), 'km')
print('나중 관측할 때의 좌표 (x,y,z) = ', Equatorial_Coordinate_System(RA_t, DEC_t, R_t), 'km')
그리고 두 점 사이의 거리를 구하기 위해
아래 그림의 파란색 글씨로 보이는 3차원에서의 거리 공식을 코드로 구현할 것입니다.
이 때의 한계점은 ISS가 타원궤도로 움직였기에 그 경로가 곡선이지만, 현재 저의 프로그래밍에 대한 지식으로는 이를 함수로 구현하기에 한계가 있어 두 시각 사이에 국제우주정거장이 움직인 경로를 직선 경로로 근사하였습니다.
이 함수의 구현에서 제가 즐겨쓰는 방법을 알려드리도록 하겠습니다.
바로 빈 list 정의 -> for문을 통해 원하는 과정을 거침 -> 원하는 과정을 거친 새로운 값을 빈 list에 append하여 원하는 list 구하기
저는 주로 이 과정을 즐겨씁니다. 왜냐하면 for문을 통해 많은 개수의 데이터들을 가공하거나 if구문을 통해 선별하여 원하는 값들만을 추려낼 수 있고, 그것을 빈 list에 append함으로써 제가 진정으로 원하는 데이터를 도출해낼 수 있기 때문입니다. 아마 다른 분들도 이 방법을 흔히 사용하실 겁니다.
이 방법과 3차원에서의 거리공식을 코드로 구현하여 함수를 만들고 출력하면 결과는 다음과 같습니다.
# 3 그 둘 사이의 거리를 구한다.
def iss_distance(iss_0, iss_t): # 관측 시간 사이의 iss이동거리를 구하는 함수
n = len(iss_0)
location_list = []
for i in range(n):
location_list.append((iss_0[i] - iss_t[i])**2) # 두 좌표의 x-x1 , y-y1, z-z1 구하고 각각을 제곱
distance = np.sqrt(sum(location_list)) # x-x1 , y-y1, z-z1 각각을의 제곱값들을 모두 더한 후 제곱근 해주어 거리 유도
return distance
iss_0 = Equatorial_Coordinate_System(RA_0, DEC_0, R_0) # 처음 관측시간의 위치
iss_t = Equatorial_Coordinate_System(RA_t, DEC_t, R_t) # 나중 관측시간의 위치
print('관측시간 사이의 거리 변화 = ', iss_distance(iss_0, iss_t), 'km')
이 때 두 지점에서의 시각을 list에 시, 분, 초에 대한 요소로 정의한 후 세 요소를 모두 초 단위로 수정하는 과정을 거치면 km/s의 속력을 유도할 수 있는 time scale이 유도됩니다. 바로 위의 코드 유도한 거리(iss_distance(iss_0, iss_t))와 이 과정으로 유도한 시간을 나누어 국제우주정거장의 평균 속력의 근삿값을 유도하였습니다.
# 4 iss의 속력을 구한다
def iss_velocity(distance, t0, t): # 거리와 처음/나중 관측 시간을 통해 iss 속도 유도하는 함수
sec_t0 = (t0[0] * 3600) + (t0[1]*60) + (t0[2]) # 처음관측시간을 초단위로 나타내기
sec_t = (t[0] * 3600) + (t[1]*60) + (t[2]) # 나중관측시간을 초단위로 나타내기
time = sec_t - sec_t0
velocity = distance / time # km/s
return velocity
t0 = [20,51,52] # 처음관측시간 [시,분,초]
t = [20,52,31] # 나중관측시간 [시,분,초]
distance = iss_distance(iss_0, iss_t)
real_v = 7.664 # 국제우주정거장 속력 이론값
print('국제우주정거장(iss) 속력의 근삿값 = ', iss_velocity(distance, t0, t), 'km/s')
print('오차율 = ', 100.*(iss_velocity(distance, t0, t)-real_v)/real_v, '%')
저는 오차율이 -0.66% 정도 나왔습니다. 오차율을 보니 일련의 과정들이 어느 정도는 정확했음을 알 수 있습니다 ^^b
제가 코딩을 사용하여 국제우주정거장의 궤도 정보를 표현한 다른 친구와 비교한 결과 아이디어는 같았지만, 그 과정과 코드에는 엄청난 차이가 있었습니다.
파이썬을 이용한 다른 코드와 달리 이번 코드는 생각할 점이 많았습니다.
단순히 데이터를 가지고 노는 것이 아닌 한정적인 데이터들을 궤도 / 비행 역학적인 개념을 적용하고,
그 적용하는 과정을 하나하나 코드로 구현해야 했습니다. 이 점에서 창의성과 풍부한 아이디어, 궤도 / 비행 역학적인 개념에 대한 이해, 파이썬 프로그래밍 언어 코드에 대한 이해 이 3박자가 어우러져야해서 이번 과제가 특별히 인상깊지 않았나 싶습니다.
이제 끝이냐고요?
놀랍게도 지금까지 인트로였습니다 ~~~ /^0^/
다음 게시물에서는 지구를 중심으로하는 타원궤도에서 위성의 운동을 기술할 수 있는 궤도/비행 역학적인 공식들에 대해서 소개하고 아직 구하지 않은 궤도 장반경, 단위 질량 당 에너지, 단위 질량 당 각운동량, 이심률,을 구하는 과정에 대해서 보여드리도록 하겠습니다.
그럼 다음 시간에 뵙겠습니다!! 읽어주셔서 감사합니다~