PyEphemの衛星軌道を描く

またまた、横道に逸れますが、PyEphemを使って

人工衛星をシミュレーションしてみたいと思います。

動詞だからシミュレートかな？

参考URLは下記になります。

URL:http://cyclodextrin.hatenablog.com/entry/2017/09/15/013738

4Kノートなのは良いのですが、文字が米粒大になったり

地図などを描く線が凄く細くなってしまうため解像度を落として作成しました。

そういえば、Jupyter notebookは大丈夫なのですが

SpyderはSpyder notebookを使おうと開いたのですが、文字が米粒大どころかもっと小さく

とても作業できる状態ではありませんでした。

（これは解像度下げ、倍率を上げても解決しませんでした）

import numpy as np
import ephem
import time
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap
from datetime import datetime as dt
from datetime import timedelta as td

# 人工衛星の軌道要素
object_name = "Nakano"       # 衛星名
epoch_date = "23313.3333333333"   # 左記の軌道要素が観測された時刻
inclination = 50.66666           # 軌道平面の傾き（度）赤道に対する傾き
RAofAN = 138.3480               # 昇交点の赤経（度, 春分点基準）：赤道面を南側から北側へ横切る位置
eccentricity = 0.0834024        # 軌道の離心率 (必ず1未満)  　　
arg_perigee = 275.5910          # 近点引数（度）：地球の重心から見た時になす角度を、
                                        # 天体の運動方向に沿って昇交点から測った角度。記号は「ω」
mean_anomaly = 359.7789         # 平均近点角(度)： 図参照
mean_motion = 1.00258630        # 平均回転数（周回/日)
decay_rate = 0.00000228         # 軌道減衰率（周回/日^2)：図参照
orbit_num = 3000                # これまでの積算周回数：周回した回数の合計

# epoch_dateをXEphem用の形式に変換
year = int(epoch_date[:2])
if (year >= 57):
    year += 1900
else:
    year += 2000
dayno = float(epoch_date[2:])
epoch_date_sl = "1/%s/%s" % (dayno, year)

# epoch_dateから一分ごとの時刻を3600個作成
decimal_time = dt(year, 1, 1) + td(dayno - 1)
time_list = [(decimal_time + td(hours=i/60)).strftime("%Y/%m/%d %H:%M:%S") for i in range(0,3600)]

# 世界地図を描写
m = Basemap()
m.drawcoastlines()

# time_listに基づいて一分ごとの人工衛星の位置をプロット
for item in time_list:
    orbit_data = ",".join(map(str, (object_name, "E", epoch_date_sl, inclination, RAofAN, eccentricity, arg_perigee, mean_anomaly, mean_motion, decay_rate, orbit_num)))
    satellite = ephem.readdb(orbit_data)
    satellite.compute('%s' % item)
    latitude = satellite.sublat / ephem.degree
    longitude = satellite.sublong / ephem.degree
    x1,y1 = m(longitude, latitude)
    m.plot(x1, y1, "m.", markersize=1.5)

# プロット表示
plt.show()