ChatGPTで系外惑星Gliese 12 bの軌道をシミュレーションしてみた（AI宇宙開発＃1）

はじめまして！大学で宇宙について学んでいる最中の大学生です。
『徒然草』ではないですが、心の赴くままになんとなくQiitaを書いてみようと思いました。

系外惑星Gliese 12 b

2024年5月、国立天文台が『宇宙生命探査の鍵となる「太陽系外の金星」を発見』と題して、太陽系外惑星Gliese 12 bを発見したと発表しました。
冒頭には、次のように書かれています。

地球からわずか40光年の距離に、太陽よりも低温の恒星のまわりを周回し、地球や金星と同程度の大きさを持つ太陽系外惑星「グリーゼ 12 b」が、発見されました。この惑星が恒星から受け取る日射量は、太陽系において金星が太陽から受けるものと同程度と考えられます。また、惑星の大気が宇宙空間に散逸せずに一定量留まっている可能性があります。これらのことから、過去に発見された系外惑星に比べて、「金星のような惑星の大気の特徴を、地球と比較しつつ調べるのに最も適した惑星」と言えます。今後、グリーゼ 12 bの大気の詳細な調査が、惑星が生命の居住に適した環境を持つための条件についての理解を、大きく進歩させると期待されます。

興味のある方は、以下のリンクから飛んでみましょう。
宇宙生命探査の鍵となる「太陽系外の金星」を発見

ChatGPTで論文を作れないか？

幸運なことに、私は大学の研究やサークル活動を通して宇宙開発に触れることができています。宇宙に関するニュースが毎日飛び込んでくると、そのたびにワクワクする変人であることも自覚しています...。
最近では、Cursorなどのエディタでシミュレーション用のプログラムを作っていますが、このCursorはvscodeにChatGPTがくっついているという衝撃的なエディタです。一旦自分でコードを書いてランすると、当然ながら多種多様なエラーが襲いかかってきます。
しかし、エラーでしょんぼりする前にコードの隣に開いたチャット欄にエラーをそのまま投げると、ChatGPTがものの数秒で修正案を出してくれます。これのおかげで、プログラムがある程度形になるまでの時間がかなり短縮されたと思います。
参照：次世代のエディター Cursor(カーソル) を使いこなす(2024年更新)

そんな折、先述の国立天文台のニュースが飛び込んできました。そこで、私はChatGPTにこの系外惑星をテーマとして研究論文（もどき）を作らせようではないかと考えました。日頃のレポートなどを敢えて誰の力も借りることなく己の固い頭で無理やりひねり出して作っていたからこそ、AIがどんな論文（もどき）を作るのかは純粋に興味がありました。

タイトル：The Orbital Dynamics of Gliese 12 b: A Closer Look at a Promising Exoplanet

随分と大袈裟なことを言っていますが、要するにAstropyとSkyfieldライブラリを用いて2024年1月1日から1月15日まで、0.5時間の時間ステップでGliese 12 bの軌道シミュレーションを行ったというだけの内容です。

使用コード

今回の論文もどきを作成するにあたって、以下のコードを走らせました。

gliese12b-simu.py

import numpy as np
from astropy import units as u
from astropy.time import Time, TimeDelta
from astropy.coordinates import SkyCoord
from skyfield.api import load, Star
from skyfield.data import hipparcos
import matplotlib.pyplot as plt

# JPLのエフェメリスデータをロード
eph = load('de421.bsp')

# グリーゼ12のデータを取得
with load.open(hipparcos.URL) as f:
    stars = hipparcos.load_dataframe(f)
gliese12_data = stars.loc[87937]  # グリーゼ12のデータ
gliese12 = Star.from_dataframe(gliese12_data)
print(gliese12_data)
# 時間スケールの設定
ts = load.timescale()

# グリーゼ12bの軌道要素
semi_major_axis = 0.1 * u.au               # 軌道長半径
eccentricity = 0.05                        # 離心率
inclination = 45 * u.deg                   # 軌道傾斜角
longitude_of_ascending_node = 100 * u.deg  # 昇交点経度
argument_of_periapsis = 250 * u.deg        # 近点引数
period = 12.76 * u.day                     # 軌道周期

def calculate_true_anomaly(time):
    """天体の公転軌道において、惑星や衛星が軌道上でどれだけ進んでいるかを示す角度を計算する関数"""
    mean_anomaly = 2 * np.pi * (time.jd % period.value) / period.value
    return mean_anomaly + 2 * eccentricity * np.sin(mean_anomaly)

def calculate_position(time):
    """グリーゼ12bの座標(x,y,z)をそれぞれ計算する関数"""
    true_anomaly = calculate_true_anomaly(time)
    r = semi_major_axis * (1 - eccentricity**2) / (1 + eccentricity * np.cos(true_anomaly))
    lon_asc_node_rad = longitude_of_ascending_node.to(u.rad).value
    arg_periapsis_rad = argument_of_periapsis.to(u.rad).value
    inclination_rad = inclination.to(u.rad).value
    
    x = r * (np.cos(lon_asc_node_rad) * np.cos(arg_periapsis_rad + true_anomaly) -
             np.sin(lon_asc_node_rad) * np.sin(arg_periapsis_rad + true_anomaly) * np.cos(inclination_rad))
    y = r * (np.sin(lon_asc_node_rad) * np.cos(arg_periapsis_rad + true_anomaly) +
             np.cos(lon_asc_node_rad) * np.sin(arg_periapsis_rad + true_anomaly) * np.cos(inclination_rad))
    z = r * np.sin(arg_periapsis_rad + true_anomaly) * np.sin(inclination_rad)
    return x, y, z

def plot_positions(start_time, end_time, time_step):
    """指定期間におけるグリーゼ12およびグリーゼ12bをそれぞれプロットする関数"""
    times = np.arange(start_time.jd, end_time.jd, time_step.jd)
    x_vals, y_vals, z_vals = [], [], []

    for t in times:
        current_time = Time(t, format='jd', scale='utc')
        x, y, z = calculate_position(current_time)
        x_vals.append(x.to(u.au).value)
        y_vals.append(y.to(u.au).value)
        z_vals.append(z.to(u.au).value)

    # gliese12_dataからRAとDecを取得
    ra = gliese12_data['ra_degrees'] * u.deg    # 赤経 [deg]
    dec = gliese12_data['dec_degrees'] * u.deg  # 赤緯 [deg]    
    # SkyCoordを使用して銀河座標を作成
    gliese12_coord = SkyCoord(ra, dec, frame='icrs')    
    # 銀河座標を取得
    galactic_coords = gliese12_coord.galactic
    galactic_x = galactic_coords.l  # 銀経 [deg]
    galactic_y = galactic_coords.b  # 銀緯 [deg]
    print(galactic_coords)

    # グリーゼ12を銀経銀緯座標で2Dプロットする
    plt.figure(figsize=(10, 6))
    plt.scatter(galactic_x, galactic_y, label='Gliese 12', color='blue')
    plt.xlim(-180, 180)
    plt.ylim(-90, 90)
    plt.xlabel('Galactic Longitude (degrees)')
    plt.ylabel('Galactic Latitude (degrees)')
    plt.title('M3.0 dwarf Gliese 12 (TOI-6251) in Galactic Coordinates')
    plt.legend()
    plt.grid()
    plt.show()

    # 指定期間にて、グリーゼ12の周りを公転するグリーゼ12bを3Dプロットする
    fig = plt.figure(figsize=(10, 6))
    ax = fig.add_subplot(111, projection='3d')
    # グリーゼ12bの軌道を滑らかにプロット
    ax.plot(x_vals, y_vals, z_vals, label='Gliese 12b Orbit', color='blue')
    # グリーゼ12を原点としてプロット
    ax.scatter(0,0,0, c='red', s=100, label='Gliese 12') 
    ax.set_xlabel('X (AU)')
    ax.set_ylabel('Y (AU)')
    ax.set_zlabel('Z (AU)')
    ax.set_title('Orbit of Gliese 12b around Gliese 12')
    ax.legend()
    plt.grid()
    plt.show()

# シミュレーション期間の設定
start_time = Time("2024-01-01T00:00:00", format='isot', scale='utc')
end_time = Time("2024-01-15T00:00:00", format='isot', scale='utc')
time_step = TimeDelta(0.5 * u.hour)  # 1日ごとのステップ

# 軌道プロット
plot_positions(start_time, end_time, time_step)

# 座標をプリントする関数
def print_positions(start_time, end_time):
    current_time = start_time
    print("Date\t\tX (AU)\t\tY (AU)\t\tZ (AU)")
    print("-------------------------------------------------")
    
    while current_time <= end_time:
        x, y, z = calculate_position(current_time)
        print(f"{current_time.iso}\t{x.to(u.au).value:.6f}\t{y.to(u.au).value:.6f}\t{z.to(u.au).value:.6f}")
        current_time += TimeDelta(1 * u.day)  # 1日進める

# 座標をプリント
print_positions(start_time, end_time)

得られるグラフ

こんな感じの2Dグラフと3Dグラフが得られると、きちんと動いていると思います。画像ファイルの添付方法が分からないため、Google Drive経由で示します。


(↑ クリックしてもダウンロードしないと見られないそうです。申し訳ありません。)

論文もどき

PDFの添付方法が分からないため、やっぱりGoogle Drive経由で示します。
The Orbital Dynamics of Gliese 12 b: A Closer Look at a Promising Exoplanet

結論

ここまで読んでくださりありがとうございます。
今後も暇な時間ができたら、趣味の範疇で#2をやってみたいです。