파이썬으로 지구와 달의 공전궤도 시뮬레이션하는 방법, 삼체 문제 (three-body problem), 만유인력법칙, 궤도역학

▶ 삼체 문제 (three-body problem)

- 고전역학에서 삼체 문제란, 서로 상호작용하는 세 물체의 운동을 다루는 문제이다. 

- N체 문제 (N-body problem)에서 N=3인 경우를 의미한다.

 

▶ 시뮬레이션 주제

- 태양-지구-달 삼체 시스템에서의 지구와 달의 공전궤도

 

▶ 시뮬레이션 설정

- 프로그래밍 언어: 파이썬3.9

- 편집기: Jupyter notebook

------------------------------------

- 물리량 단위: MKS units (미터(metre, m), 킬로그램(kilogram, kg), 초(second, s))

- 좌표계: 데카르트 좌표계 (2D, (x, y))

- time step: 1 day

- number of steps : 3 years (365x3 days)

 

▶ 시뮬레이션 결과

red dot: 태양, blue dot: 지구, green dashed line: 달의 공전궤도

 

태양과 지구 사이의 거리가 달의 공전궤도에 비해 매우 크기 때문에 달의 공전을 확인할 수 없다.

달의 공전을 확인하기 위해 지구를 frame의 중앙에 놓고 애니메이션을 실행했고, 그 결과는 아래와 같다.

 

time step: 1 hour, number of time steps: 24*27 (27 days)

▶ 코드

### 0. 라이브러리 불러오기 ###
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import matplotlib.image as mpimg

%matplotlib notebook

### 1. 가시화를 위해 두 천체의 위치 좌표를 시뮬레이션 스텝마다 저장
trajectory_earth = []
trajectory_sun = []
trajectory_moon = []

### 2. 물리량 설정 ###
G = 6.7e-11 # 중력상수 G, J*m/kg^2
m_sun = 2e30 # 태양 질량 설정, kg
m_earth = 6.0e24 # 지구 질량 설정, kg
m_moon = 7.4e22 # 달 질량 설정, kg

### 3. 초기 조건 설정 ###
dist_sun_to_earth = 1.5e11 # 초기 태양 - 지구 거리 설정, m
dist_earth_to_moon = 3.85e8 # 초기 지구 - 달 거리 설정, m
dist_sun_to_moon = dist_sun_to_earth + dist_earth_to_moon # 초기 지구 - 달 거리 설정, m
x_sun = np.array([0., 0.], dtype='float64') # 태양 초기 위치 설정, m
x_earth = np.array([dist_sun_to_earth, 0.], dtype='float64') # 지구 초기 위치 설정, m
x_moon = np.array([dist_sun_to_earth + dist_earth_to_moon, 0.], dtype='float64') # 지구 초기 위치 설정, m
vel_sun = np.array([0., 0.], dtype='float64') # 태양 초기 속도 설정, m/s
vel_earth = np.array([0., 2 * np.pi * dist_sun_to_earth / (365 * 24 * 3600)], dtype='float64') # 지구 초기 속도 설정, m/s
vel_moon = np.array([0., 2 * np.pi * dist_sun_to_earth / (365 * 24 * 3600) + 2 * np.pi * dist_earth_to_moon / (27 * 24 * 3600)], dtype='float64') # 달 초기 속도 설정, m/s

trajectory_sun.append(x_sun.copy())
trajectory_earth.append(x_earth.copy())
trajectory_moon.append(x_moon.copy())

### 4. 시뮬레이션 환경 설정 ###
num_steps = 27 * 24 # number of steps : 27 days
dt = 60 * 60 # time step : 1 hour

### 5. 시뮬레이션 작동 (시뮬레이션 결과 저장) ###
for i in range(num_steps):
    # 각 천체 간의 거리 계산
    dist_sun_to_earth = np.linalg.norm(trajectory_earth[i] - trajectory_sun[i])
    dist_earth_to_moon = np.linalg.norm(trajectory_moon[i] - trajectory_earth[i])
    dist_sun_to_moon = dist_sun_to_earth + dist_earth_to_moon

    # 중력에 의한 힘 계산
    force_sun_to_earth = (G * m_sun * m_earth / dist_sun_to_earth ** 2) * (trajectory_earth[i] - trajectory_sun[i]) / dist_sun_to_earth
    force_earth_to_moon = (G * m_earth * m_moon / dist_earth_to_moon ** 2) * (trajectory_moon[i] - trajectory_earth[i]) / dist_earth_to_moon
    force_sun_to_moon = (G * m_sun * m_moon / dist_sun_to_moon ** 2) * (trajectory_moon[i] - trajectory_sun[i]) / dist_sun_to_moon

    # 가속도 계산
    acc_sun = (force_sun_to_earth + force_sun_to_moon) / m_sun
    acc_earth = (-force_sun_to_earth + force_earth_to_moon) / m_earth
    acc_moon = -(force_earth_to_moon + force_sun_to_moon) / m_moon

    # 속도 업데이트
    vel_sun += acc_sun * dt
    vel_earth += acc_earth * dt
    vel_moon += acc_moon * dt

    # 위치 업데이트
    x_sun += vel_sun * dt
    x_earth += vel_earth * dt
    x_moon += vel_moon * dt
    
    # 결과 저장
    trajectory_sun.append(x_sun.copy())
    trajectory_earth.append(x_earth.copy())
    trajectory_moon.append(x_moon.copy())

# numpy 배열로 변환
trajectory_sun = np.array(trajectory_sun)
trajectory_earth = np.array(trajectory_earth)
trajectory_moon = np.array(trajectory_moon)

### 6. 시뮬레이션 결과 시각화 ###
fig, ax = plt.subplots(figsize=(7, 7))
ax.set_xlim(-6e8, 6e8)
ax.set_ylim(-6e8, 6e8)
ax.set_xlabel('x [m]', fontsize = 15, labelpad = 7)
ax.set_ylabel('y [m]', fontsize = 15, labelpad = 7)
ax.tick_params(axis = 'x', labelsize = 13)
ax.tick_params(axis = 'y', labelsize = 13)
ax.set_aspect('equal')

# 배경색을 검정색으로 설정
ax.set_facecolor('black')

# 이미지를 불러옵니다.
earth_img = mpimg.imread('')
moon_img = mpimg.imread('')

# 이미지 크기 설정
earth_size = 1.2e7  # 지구 이미지 크기
moon_size = 3.4e6   # 달 이미지 크기

# 태양과 지구 이미지 초기 위치 설정 (임의로 초기 위치 설정)
earth = ax.imshow(earth_img, extent=(dist_sun_to_earth - earth_size, dist_sun_to_earth + earth_size, -earth_size, earth_size))
moon = ax.imshow(moon_img, extent=(dist_sun_to_moon - moon_size, dist_sun_to_moon + moon_size, -moon_size, moon_size))


trail_earth, = ax.plot([], [], 'b--', lw=0.5)
trail_moon, = ax.plot([], [], 'g--', lw=0.5)

def update(frame):
    earth.set_extent((trajectory_earth[frame, 0] - earth_size, trajectory_earth[frame, 0] + earth_size,
                    trajectory_earth[frame, 1] - earth_size, trajectory_earth[frame, 1] + earth_size))
    moon.set_extent((trajectory_moon[frame, 0] - moon_size, trajectory_moon[frame, 0] + moon_size,
                      trajectory_moon[frame, 1] - moon_size, trajectory_moon[frame, 1] + moon_size))
    trail_earth.set_data(trajectory_earth[:frame, 0], trajectory_earth[:frame, 1])
    trail_moon.set_data(trajectory_moon[:frame, 0], trajectory_moon[:frame, 1])
    
    # 축의 위치를 현재 지구의 위치로 업데이트
    ax.set_xlim(-6e8 + trajectory_earth[frame, 0], 6e8 + trajectory_earth[frame, 0])
    ax.set_ylim(-6e8 + trajectory_earth[frame, 1], 6e8 + trajectory_earth[frame, 1])
    
    return earth, moon, trail_earth, trail_moon


ani = FuncAnimation(fig, update, frames=num_steps, interval=1, blit=True)

output_path = ''
ani.save(output_path, writer='pillow', fps=30)

plt.show()