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

▶ 이체 문제 (two-body problem)

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

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

 

↓ 태양-지구-달 삼체 시스템에서의 지구와 달의 공전궤도 시뮬레이션 프로젝트 

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

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

 

▶ 시뮬레이션 주제

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

 

▶ 시뮬레이션 설정

- 프로그래밍 언어: 파이썬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)

 

▶ 시뮬레이션 결과

time step: 1day, number of steps: 2 years

▶ 코드 설명

0. 라이브러리 불러오기

- 시뮬레이션에 필요한 라이브러리를 불러옵니다. 

- 시뮬레이션 결과를 애니메이션으로 확인하기 위해 matplotlib에서 animation에 대한 기능을 가져왔습니다.

 

1. 시뮬레이션 결과를 저장할 배열 생성

- 시뮬레이션 결과를 저장할 배열을 미리 정의합니다.

- 시뮬레이션 step 마다의 지구와 태양의 2차원 위치 정보가 이에 해당합니다.

 

2. 물리량 설정

- 중력 계산을 위한 중력 상수 G와 태양, 지구의 질량을 설정합니다.

 

$ G = 6.7\times 10^{-11} Jm/kg^{2} $

$ m_{sun}=2.0\times10^{30}kg $

$ m_{earth}=6.0\times10^{24}kg $

 

3. 초기 조건 설정

- 시뮬레이션의 초기 조건을 설정합니다.

- 태양과 지구의 초기 위치, 속도를 설정합니다.

 

$ (x,y)_{sun} = (0,0) $ [m]

$ (x,y)_{earth} = (1.5\times 10^{11},0) $ [m]

$ v_{sun} = [0, 0] $ [m/s]

$ v_{earth} = [0, 30000] $ [m/s]

 

※ 실제 태양은 질량 중심을 기준으로 운동하지만, 그 중심이 태양과 거의 일치하므로, 태양의 초기 속도를 (0,0) [m/s]으로 가정했습니다.

 

※ 지구의 초기 속도의 경우, 태양과 지구 사이의 거리를 기반으로 원운동을 가정하여 30000 m/s로 설정했습니다.

 

4. 시뮬레이션 환경 설정

- 시뮬레이션 환경을 설정합니다.

- 계산이 이뤄지는 시간 간격(time step)과 그 간격의 수(number of time steps)를 정의합니다.

 

time step = 1 day (60*60*24 seconds)

number of time steps = 2 years (365*2 days)

 

5. 시뮬레이션 실행

- 시뮬레이션을 작동시킵니다.

- 시뮬레이션은 4.에서 정의한 time step 수만큼 작동하며, 중력 계산을 통해 산출된 값을 1.에서 정의한 배열에 저장합니다.

 

6. 시뮬레이션 결과 시각화

- 시뮬레이션 결과를 확인합니다.

 

▶ 코드

### 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 = []
dist_sun_to_earth_arr = []

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

### 3. 초기 조건 설정 ###
dist_sun_to_earth = 1.5e11 # 초기 태양 - 지구 거리 설정, m
x_sun = np.array([0., 0.], dtype='float64') # 태양 초기 위치 설정, m
x_earth = np.array([dist_sun_to_earth, 0.], dtype='float64') # 지구 초기 위치 설정, m

vel_sun = np.array([0., 0.], dtype='float64') # 태양 초기 속도 설정, m/s
vel_earth = np.array([0., 30000.], dtype='float64') # 지구 초기 속도 설정, m/s

trajectory_sun.append(x_sun.copy())
trajectory_earth.append(x_earth.copy())
dist_sun_to_earth_arr.append(dist_sun_to_earth)

### 4. 시뮬레이션 환경 설정 ###
num_steps = 365 * 2
dt = 60 * 60 * 24

### 5. 시뮬레이션 실행 ###
for i in range(num_steps):

    dist_sun_to_earth = np.linalg.norm(trajectory_earth[i] - trajectory_sun[i])

    force = (G * m_sun * m_earth / dist_sun_to_earth ** 2) * (trajectory_earth[i] - trajectory_sun[i]) / dist_sun_to_earth

    acc_sun = force / m_sun
    acc_earth = -force / m_earth

    vel_sun += acc_sun * dt
    vel_earth += acc_earth * dt

    x_sun += vel_sun * dt
    x_earth += vel_earth * dt

    trajectory_sun.append(x_sun.copy())
    trajectory_earth.append(x_earth.copy())
    dist_sun_to_earth_arr.append(dist_sun_to_earth)

trajectory_sun = np.array(trajectory_sun)
trajectory_earth = np.array(trajectory_earth)
dist_sun_to_earth_arr = np.array(dist_sun_to_earth_arr)

### 6. 시각화 ###
fig, ax = plt.subplots(figsize=(7, 7))
ax.set_xlim(-2e11, 2e11)
ax.set_ylim(-2e11, 2e11)
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')

bg_img = mpimg.imread('')
ax.imshow(bg_img, extent=[-2e11, 2e11, -2e11, 2e11], aspect='auto')

# 태양과 지구 이미지 로드
sun_img = mpimg.imread('')
earth_img = mpimg.imread('')

# 이미지 크기 설정
sun_size = 1.4e10  # 태양 이미지 크기
earth_size = 3e9   # 지구 이미지 크기

# 태양과 지구 이미지 초기 위치 설정
sun = ax.imshow(sun_img, extent=(-sun_size, sun_size, -sun_size, sun_size))
earth = ax.imshow(earth_img, extent=(dist_sun_to_earth - earth_size, dist_sun_to_earth + earth_size, -earth_size, earth_size))

trail_earth, = ax.plot([], [], 'b--', lw=1)
trail_sun, = ax.plot([], [], 'r-', lw=1)

def update(frame):
    sun.set_extent((trajectory_sun[frame, 0] - sun_size, trajectory_sun[frame, 0] + sun_size,
                    trajectory_sun[frame, 1] - sun_size, trajectory_sun[frame, 1] + sun_size))
    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))
    trail_earth.set_data(trajectory_earth[:frame, 0], trajectory_earth[:frame, 1])
    trail_sun.set_data(trajectory_sun[:frame, 0], trajectory_sun[:frame, 1])
    return sun, earth, trail_earth, trail_sun

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

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

plt.show()