파이썬을 이용하여 케플러 제 3법칙을 그래프로 나타내기

이번 포스팅은 아주 오래전에 제가 올렸었는데 너무 대충 포스팅한 것 같아 다시 만들려고 하는 포스팅입니다.(ㅠㅠ)

  1. 개요

케플러 제 3법칙은 태양계 행성 뿐 아니라 우주 전체의 궤도를 가지는 행성이나 쌍성 등에 적용되는 매우 유명한 법칙 중 하나입니다. 물리교과는 물론이고 지구과학2에서도 매우 중요하게 다루어지는 내용이기도 합니다. 여기서는 케플러 제3법칙에 대한 이야기를 하고, 다음 천문 포스팅을 할 때 케플러 법칙에 대한 포스팅을 해 볼 예정입니다.

 

케플러 제3법칙은 다들 아시는 것처럼 아래의 관계를 가집니다.

 

$$P^{2}=k×a^{3}$$

 

여기서 P는 행성의 공전궤도 주기를, a는 궤도 장반경을 의미하는데, k의 경우 태양계의 경우 k=1 이 되기 때문에 위 식은 간단히 아래와 같이 표현할 수 있습니다.

 

$$P^{2}=a^{3}$$

 

또한 위 식이 성립하기 위해서는 공전궤도 주기를 년 단위로, 장반경의 거리 단위를 AU(천문단위, astronomical unit)를 사용할 때 k=1로 성립됩니다. 태양계에서 k=1이 되는건 이번 포스팅의 주제를 벗어나기 때문에 설명하지 않고, 케플러 법칙에 대한 다음 포스팅에서 설명하겠습니다.

 

그래프를 표현하기 위해서는 데이터가 필요합니다. 행성의 공전궤도와 공전주기 데이터는 아래 표를 참고하세요.

	수성	금성	지구	화성	목성	토성	천왕성	해왕성
주기(year)	0.24	0.615	1	1.88	11.86	29.457	84.022	164.8
궤도반지름 (AU)	0.39	0.7	1	1.52	5	10	20	30

  2. 그래프 형식 살펴보기

우리에게 익숙한 케플러 제 3법칙은 그래프로 표현하면 아래와 같습니다.

반응형

그래프를 잘 보면 이상한게 있습니다. x축과 y축의 간격이 처음에는 10배, 100배, 1000배로 등간격이 아닙니다. 눈치채신분이 있겠지만, 케플러 제3법칙에서 x축과 y축은 로그스케일 입니다. 파이썬에서 로그스케일로 그래프를 표현하는 방법은 여러가지가 있는데, 여기서는 matplotlib.pyplot의 명령어 중 하나인 plt.xscale을 사용하여 아래와 같은 것을 쓸것입니다.

 

plt.xscale('log')

plt.yscale('log')

 

그리고 또 한가지 주의사항.

 

P에 행성의 공전궤도를 대입하여 계산된 궤도 반지름을 이용하면 문제가 없지만, 알려진 주기와 궤도반지름을 각각 대응하여 그래프를 그리면

위의 그래프처럼 나옵니다. 누가봐도 삐뚤빼뚤합니다. 이건 당연하겠지만, 실제 궤도 주기와 궤도반지름은 케플러 제3법칙을 적용하여 그대로 나타냈을 때 상호간에 약간의 오차가 있기 때문입니다. 그래서 실제 데이터를 이용해 주기에 제곱을, 궤도반지름에 세제곱을 하여 그래프로 대응하기 위해서는 추세선을 이용하는게 더 낫습니다.

 

  3. 전체 코드

전체 코드를 보며 코드에 대한 설명을 하나씩 하겠습니다.

복붙하여 사용하실 수 있게 아래에도 동일한 내용을 적겠습니다.

import matplotlib.pyplot as plt  
import numpy as np
from scipy.optimize import curve_fit



def objective(x, a,b):
    return a*np.array(x)+b
period=[0.2, 0.6, 1, 1.9, 11.9, 29.5, 84.0, 164.8]
radius=[0.4, 0.7, 1, 1.5, 5, 10, 20, 30]
period2=[]
radius2=[]
for x in period:
    y=x**2
    period2.append(y)
for z in radius:
    a=z**3
    radius2.append(a)
popt,pcov = curve_fit(objective, period2, radius2,method='dogbox', bounds=([0.9999999,-0.000001],[1.0000001,0.000001]))
c,d=popt
plt.figure(figsize=(8,6))
plt.plot(period2,objective(period2,*popt),color='red',linestyle='--')
plt.scatter(period2,radius2)
plt.xscale('log')
plt.yscale('log')
plt.grid()
plt.title("Kepler's 3rd law", size=15)
plt.xlabel('$P^{2}$($AU^{2}$)', size=15)
plt.ylabel('$a^{3}$($year^{3}$)', size=15)
plt.savefig("C:\\111\\kepler.jpg", dpi=300)

plt.show()

  4. 코드 분석

이제 이 내용을 하나씩 보면

import matplotlib.pyplot as plt  
import numpy as np
from scipy.optimize import curve_fit

필수 라이브러리를 불러옵니다. 추세선을 그리기 위해 가장 아래에 curve_fit을 호출하는 코드를 사용한 것에 주목하세요.

def objective(x, a,b):
    return a*np.array(x)+b

 

추세선을 그리기 위한 1차 함수를 def로 정의하였습니다.

그런데 여기서 보통 return 뒤에 나오는 반환값은 curve_fit을 이용하는 경우 a*x+b라고만 하는데, 뒤에 코드에서 계산하는 소수가 부동소수점이라서 오류가 뜹니다. 그래서 부동소수점을 계산할 수 있도록 x를 np.array(x)라고 행렬형으로 집어넣어버렸습니다. 부동소수점에 대한 소개는 본 글의 목적에 한참 어긋나기 때문에 생략합니다.

 

아무튼 다음 코드는

period=[0.2, 0.6, 1, 1.9, 11.9, 29.5, 84.0, 164.8]
radius=[0.4, 0.7, 1, 1.5, 5, 10, 20, 30]
period2=[]
radius2=[]
for x in period:
    y=x**2
    period2.append(y)
for z in radius:
    a=z**3
    radius2.append(a)

 

period라는 리스트에 행성의 주기를,

radius라는 리스트에는 행성의 반지름을 집어넣었고, 각 리스트의 값을 for문을 이용하여 제곱과 세제곱하여 계산한 결과를 집어넣을 빈 리스트 period2와 radius2를 만들었습니다. 이렇게 하여 P² 과 A³ 의 계산 결과를 모두 period2와 raduis2라는 리스트에 모조리 집어넣었습니다.

 

이제부터 curve_fit을 이용하여 본격적으로 추세선을 그릴 단계입니다. 어려우니까 주목!!

 

popt,pcov = curve_fit(objective, period2, radius2,method='dogbox', bounds=([0.9999999,-0.000001],[1.0000001,0.000001]))
c,d=popt
plt.figure(figsize=(8,6))
plt.plot(period2,objective(period2,*popt),color='red',linestyle='--')

 

지난번 허블 - 르메트르의 법칙에서도 포스팅을 했지만 다시한번 이야기 해 보겠습니다.

 

popt,pcov에 curve_fit으로 objective에 있는 맞춤(fitting)함수 a*x+b를 이용해 period2와 radius2에 대한 공분산 값을 계산하여, 가장 정확하게 맞춤하였을 때 1차함수의 a, b 값을 계산하라는 함수가 첫 번째 코드에 나와있습니다. 지굼의 경우 popt에 들어갈 값은 a, b에 해당하겠지만, 2차함수일 경우 3개, 3차함수일 경우 4개가 될 수도 있습니다. 어째되었든 popt에는 맞춤했을 때 공분산 값이 최소인 a,b의 값이 지정되고, pcov는 그 때 공분산 값이 지정됩니다.

 

첫 번째 코드 맨 뒤에 bounds라는 값을 주었습니다. 저렇게 하여 a,b의 최소값과 최대값을 지정해 줍니다. 첫번째 리스트형 데이터는 a, b의 최소값을, 두 번째 리스트형 데이터에는 최대값을 지정해 주었습니다.

 

다음에는 c,d=popt라고 하였는데, c,d는 각각 계산된 a,b의 값이 얼마인지 알고자 집어넣었습니다. 궂이 안넣어도 되는 코드이고

 

plt figure는 그냥 그래프 크기 지정이니 넘어가겠습니다.

 

plt.plot 다음에 나오는건, 첫 번째로 x축 값으로는 period2를, y축 값으로는 objective 함수에 x축 값으로 period2를 집어놓고, curve_fit으로 맞춤했을 때 나온 최적의 새로운 ax+b함수에 period2를 집어넣었을 때 계산된 y축 값에 해당합니다. 지난번에도 이야기 했지만, *popt라고 쓴 건 값이 2개 이상일 수 있기 때문입니다. color는 단순하게 그래프 색을, linestyle은 그냥 그래프 형식입니다.

 

이제 밑으로는 그냥 단순히 그래프 꾸미는 코드입니다. 다만 주목할 것은 그래프의 x, y축이 로그 스케일이 되도록 하는 코드  

 

plt.xscale('log')
plt.yscale('log')

 

를 집어넣었구요, x, y축 제목이 위첨자가 들어가는 수식이어서 다음과 같이 표현하였습니다.

 

plt.xlabel('$P^{2}$($AU^{2}$)', size=15)
plt.ylabel('$a^{3}$($year^{3}$)', size=15)

 

코드를 자세히 보시면 아마 코드의 형식이 눈에 들어올 것입니다.

 

이렇게 하여 만들어진 최종 그래프는

  5. 추체선을 그리지 않고 1차 함수 형태로 나오게 하는 방법은 없을까?

 

있습니다. y축 값으로 실제 행성 공전궤도를 사용하는게 아니고, x축 값을 바탕으로 계산한 값을 쓰면 됩니다. 그럼 리스트형 데이터는 period 하나만 있으면 됩니다.

 

코드는 훨씬 단순해 질 것입니다.

 

import matplotlib.pyplot as plt
period=[0.2, 0.6, 1, 1.9, 11.9, 29.5, 84.0, 164.8]
radius=[]
for x in period:
    a=x**(2/3)
    radius.append(a)
period2=[]
radius2=[]
for x2 in period:
    x3=x2**2
    period2.append(x3)
for y2 in radius:
    y3=y2**3
    radius2.append(y3)
plt.figure(figsize=(7,5))
plt.grid('--')
plt.xscale('log')
plt.yscale('log')
plt.plot(period2,radius2, '--o')
plt.xlabel('$P^{2}$($year^{2}$)',size=15)
plt.ylabel('$a^{3}$($AU^{3}$})', size=15)
plt.title("Kepler's 3rd law", size=17)
plt.savefig("C:\\111\\kepler 2nd.jpg")
plt.show()

포스팅을 하느라 방금 5분만에 급조한 코드라 머리를 안쓰고 후다닥 하다보니 코드가 좀 길어졌지만, curve_fit을 하는 것 보다 생각할 요소는 훨씬 적습니다. 이렇게 하면, 공전궤도 긴반지름 값은 실제 값과 다소 차이가 있겠지만, 그래프는 오차없는 1차함수 형태가 됩니다. 어쨌든 이렇게 계산된 공전궤도 반지름 값을 보면 실제 데이터와 다소 오차가 있습니다. 그래프 형태는 아래와 같습니다.

어떤 형식으로 그래프를 그릴지는 개인 취향일 것 같습니다.

 

  6. 마치며

curve_fit은 본 포스팅의 내용을 눈으로 꼼꼼히만 읽으면 이해하기가 매우 어렵습니다. 우선 복붙하여 코드를 돌려보신 다음에, 이것저것 바꿔가며 연구해 보면 금새 이해되실 것이라 생각합니다. 물론 굉장히 어려운 내용이긴 합니다만, 충분히 해 내실 수 있을 것이라 생각합니다.

 

다음번 천문 관련 포스팅에서는 케플러 제3법칙에 대한 것을 작성해 보겠습니다. 긴 글 읽느라 고생하셨어요~ㅎㅎ : )