PyEphemによる人工衛星の軌道計算と可視化 — ISSの追跡からパス予測まで

ISS (ZARYA)
1 25544U 98067A   24045.51782528  .00011186  00000-0  20208-3 0  9993
2 25544  51.6412 208.6063 0005140  43.4286  70.0747 15.49562695441075

pip install ephem

import ephem

import ephem
from datetime import datetime

# 文字列から生成（UTC）
d1 = ephem.Date("2026/4/25 12:00:00")

# datetimeから生成（UTC）
d2 = ephem.Date(datetime(2026, 4, 25, 12, 0, 0))

# 現在時刻を取得
d3 = ephem.now()

print(f"d1 = {d1}")
print(f"d2 = {d2}")
print(f"d3 = {d3}")

# ephem.Date -> datetime に変換
dt = ephem.Date(d1).datetime()
print(f"datetime型: {dt}")

import ephem

# ISSのTLEデータ（サンプル）
tle_name = "ISS (ZARYA)"
tle_line1 = "1 25544U 98067A   24045.51782528  .00011186  00000-0  20208-3 0  9993"
tle_line2 = "2 25544  51.6412 208.6063 0005140  43.4286  70.0747 15.49562695441075"

# 衛星オブジェクトを作成
iss = ephem.readtle(tle_name, tle_line1, tle_line2)

# 特定の時刻で位置を計算
iss.compute("2024/2/15 12:00:00")

# 結果を表示
print(f"衛星名: {iss.name}")
print(f"緯度:   {iss.sublat}")          # 直下点の緯度
print(f"経度:   {iss.sublong}")         # 直下点の経度
print(f"高度:   {iss.elevation / 1000:.1f} km")  # 高度 (mからkmに変換)
print(f"赤経:   {iss.ra}")              # 赤経（天球座標）
print(f"赤緯:   {iss.dec}")             # 赤緯（天球座標）

import ephem
import math

iss = ephem.readtle(
    "ISS (ZARYA)",
    "1 25544U 98067A   24045.51782528  .00011186  00000-0  20208-3 0  9993",
    "2 25544  51.6412 208.6063 0005140  43.4286  70.0747 15.49562695441075"
)
iss.compute("2024/2/15 12:00:00")

# ラジアン値
lat_rad = float(iss.sublat)
lon_rad = float(iss.sublong)

# 度に変換
lat_deg = math.degrees(lat_rad)
lon_deg = math.degrees(lon_rad)

print(f"緯度: {iss.sublat} = {lat_deg:.4f}°")
print(f"経度: {iss.sublong} = {lon_deg:.4f}°")

import ephem
import math
from datetime import datetime, timedelta

# ISSのTLEデータ
iss = ephem.readtle(
    "ISS (ZARYA)",
    "1 25544U 98067A   24045.51782528  .00011186  00000-0  20208-3 0  9993",
    "2 25544  51.6412 208.6063 0005140  43.4286  70.0747 15.49562695441075"
)

# 計算期間: エポック付近の1周回（93分間）
start = datetime(2024, 2, 15, 0, 0, 0)
dt_minutes = 1  # 1分間隔
n_steps = 93    # 93ステップ（約1周回）

lats = []
lons = []
alts = []
times = []

for i in range(n_steps):
    t = start + timedelta(minutes=i * dt_minutes)
    iss.compute(t)
    lats.append(math.degrees(float(iss.sublat)))
    lons.append(math.degrees(float(iss.sublong)))
    alts.append(iss.elevation / 1000)  # km
    times.append(t)

print(f"計算点数: {len(lats)}")
print(f"緯度範囲: {min(lats):.2f}° ～ {max(lats):.2f}°")
print(f"経度範囲: {min(lons):.2f}° ～ {max(lons):.2f}°")
print(f"高度範囲: {min(alts):.1f} km ～ {max(alts):.1f} km")

import ephem
import math
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime, timedelta

# ISSのTLE
iss = ephem.readtle(
    "ISS (ZARYA)",
    "1 25544U 98067A   24045.51782528  .00011186  00000-0  20208-3 0  9993",
    "2 25544  51.6412 208.6063 0005140  43.4286  70.0747 15.49562695441075"
)

start = datetime(2024, 2, 15, 0, 0, 0)
minutes = np.arange(0, 93, 0.5)  # 0.5分刻みで1周回
alts = []

for m in minutes:
    t = start + timedelta(minutes=float(m))
    iss.compute(t)
    alts.append(iss.elevation / 1000)

plt.figure(figsize=(10, 4))
plt.plot(minutes, alts, color="cyan", linewidth=1.5)
plt.xlabel("Time [min]")
plt.ylabel("Altitude [km]")
plt.title("ISS Altitude over One Orbit")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

import ephem
import math
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime, timedelta

# ISSのTLE
iss = ephem.readtle(
    "ISS (ZARYA)",
    "1 25544U 98067A   24045.51782528  .00011186  00000-0  20208-3 0  9993",
    "2 25544  51.6412 208.6063 0005140  43.4286  70.0747 15.49562695441075"
)

# 3周回分（約280分）を計算
start = datetime(2024, 2, 15, 0, 0, 0)
minutes = np.arange(0, 280, 0.5)

lats = []
lons = []

for m in minutes:
    t = start + timedelta(minutes=float(m))
    iss.compute(t)
    lats.append(math.degrees(float(iss.sublat)))
    lons.append(math.degrees(float(iss.sublong)))

lats = np.array(lats)
lons = np.array(lons)

# 経度の不連続点を処理（-180/+180の境界）
# 大きなジャンプがある箇所にNaNを挿入
lon_diff = np.abs(np.diff(lons))
jump_idx = np.where(lon_diff > 300)[0]
for idx in sorted(jump_idx, reverse=True):
    lats = np.insert(lats, idx + 1, np.nan)
    lons = np.insert(lons, idx + 1, np.nan)

# 描画
fig, ax = plt.subplots(figsize=(14, 7))
ax.plot(lons, lats, color="orange", linewidth=1.0, label="ISS Ground Track")

# 開始点をマーク
ax.plot(lons[0], lats[0], "o", color="lime", markersize=8, label="Start")

# 地図のグリッドと範囲
ax.set_xlim(-180, 180)
ax.set_ylim(-90, 90)
ax.set_xlabel("Longitude [deg]")
ax.set_ylabel("Latitude [deg]")
ax.set_title("ISS Ground Track (3 orbits)")
ax.set_aspect("equal")
ax.grid(True, alpha=0.3)

# 赤道線
ax.axhline(y=0, color="white", linewidth=0.5, alpha=0.5)

# 傾斜角の限界線
ax.axhline(y=51.64, color="cyan", linewidth=0.5, linestyle="--", alpha=0.5, label="Inclination limit")
ax.axhline(y=-51.64, color="cyan", linewidth=0.5, linestyle="--", alpha=0.5)

ax.legend(loc="lower left")
plt.tight_layout()
plt.show()

import ephem
import math
import numpy as np
import matplotlib.pyplot as plt

try:
    import cartopy.crs as ccrs
    import cartopy.feature as cfeature
    HAS_CARTOPY = True
except ImportError:
    HAS_CARTOPY = False
    print("Cartopyがインストールされていません。pip install cartopy でインストールしてください。")

from datetime import datetime, timedelta

if HAS_CARTOPY:
    # ISSのTLE
    iss = ephem.readtle(
        "ISS (ZARYA)",
        "1 25544U 98067A   24045.51782528  .00011186  00000-0  20208-3 0  9993",
        "2 25544  51.6412 208.6063 0005140  43.4286  70.0747 15.49562695441075"
    )

    start = datetime(2024, 2, 15, 0, 0, 0)
    minutes = np.arange(0, 280, 0.5)
    lats, lons = [], []

    for m in minutes:
        t = start + timedelta(minutes=float(m))
        iss.compute(t)
        lats.append(math.degrees(float(iss.sublat)))
        lons.append(math.degrees(float(iss.sublong)))

    lats = np.array(lats)
    lons = np.array(lons)

    fig = plt.figure(figsize=(14, 7))
    ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
    ax.set_global()

    # 地図の装飾
    ax.add_feature(cfeature.LAND, facecolor="darkgreen", alpha=0.3)
    ax.add_feature(cfeature.OCEAN, facecolor="darkblue", alpha=0.3)
    ax.add_feature(cfeature.COASTLINE, linewidth=0.5, edgecolor="gray")
    ax.gridlines(draw_labels=True, alpha=0.3)

    # 地上軌跡をプロット（セグメント分割で不連続対応）
    segments = []
    seg_lats, seg_lons = [lats[0]], [lons[0]]
    for i in range(1, len(lons)):
        if abs(lons[i] - lons[i-1]) > 300:
            segments.append((seg_lons, seg_lats))
            seg_lats, seg_lons = [], []
        seg_lats.append(lats[i])
        seg_lons.append(lons[i])
    segments.append((seg_lons, seg_lats))

    for seg_lon, seg_lat in segments:
        ax.plot(seg_lon, seg_lat, color="orange", linewidth=1.5,
                transform=ccrs.Geodetic())

    ax.plot(lons[0], lats[0], "o", color="lime", markersize=10,
            transform=ccrs.PlateCarree(), label="Start")

    ax.set_title("ISS Ground Track on World Map")
    ax.legend(loc="lower left")
    plt.tight_layout()
    plt.show()

import ephem
from datetime import datetime

# ISSのTLE
iss = ephem.readtle(
    "ISS (ZARYA)",
    "1 25544U 98067A   24045.51782528  .00011186  00000-0  20208-3 0  9993",
    "2 25544  51.6412 208.6063 0005140  43.4286  70.0747 15.49562695441075"
)

# 観測地点を設定（東京）
observer = ephem.Observer()
observer.lat = "35.6762"    # 北緯 35.6762°
observer.lon = "139.6503"   # 東経 139.6503°
observer.elevation = 40     # 標高 40 m
observer.date = "2024/2/15 00:00:00"  # UTC

# 次のパスを予測
pass_info = observer.next_pass(iss)

# 結果を展開
rise_time, rise_az, max_alt_time, max_alt, set_time, set_az = pass_info

print("=== 次回のISSパス予測（東京）===")
print(f"AOS (出現):   {ephem.Date(rise_time).datetime()} UTC")
print(f"  方位角:     {rise_az}  ({float(rise_az) * 180 / 3.14159:.1f}°)")
print(f"TCA (最接近): {ephem.Date(max_alt_time).datetime()} UTC")
print(f"  最大仰角:   {max_alt}  ({float(max_alt) * 180 / 3.14159:.1f}°)")
print(f"LOS (消失):   {ephem.Date(set_time).datetime()} UTC")
print(f"  方位角:     {set_az}  ({float(set_az) * 180 / 3.14159:.1f}°)")
print(f"パス時間:     {(set_time - rise_time) * 24 * 60:.1f} 分")

import ephem
import math

# ISSのTLE
iss = ephem.readtle(
    "ISS (ZARYA)",
    "1 25544U 98067A   24045.51782528  .00011186  00000-0  20208-3 0  9993",
    "2 25544  51.6412 208.6063 0005140  43.4286  70.0747 15.49562695441075"
)

# 観測地点（東京）
observer = ephem.Observer()
observer.lat = "35.6762"
observer.lon = "139.6503"
observer.elevation = 40

# 24時間以内のパスをすべて列挙
observer.date = "2024/2/15 00:00:00"
end_date = ephem.Date("2024/2/16 00:00:00")

passes = []
while observer.date < end_date:
    try:
        pass_info = observer.next_pass(iss)
        rise_time, rise_az, max_alt_time, max_alt, set_time, set_az = pass_info

        if rise_time is None:
            break

        passes.append({
            "rise_time": ephem.Date(rise_time).datetime(),
            "max_alt_deg": math.degrees(float(max_alt)),
            "set_time": ephem.Date(set_time).datetime(),
            "duration_min": (set_time - rise_time) * 24 * 60,
            "rise_az_deg": math.degrees(float(rise_az)),
            "set_az_deg": math.degrees(float(set_az)),
        })

        # 次のパスを探すために時刻を進める
        observer.date = set_time + ephem.minute
    except Exception:
        break

print(f"24時間以内のパス数: {len(passes)}\n")
for i, p in enumerate(passes, 1):
    print(f"パス {i}:")
    print(f"  出現: {p['rise_time'].strftime('%H:%M:%S')} UTC  "
          f"方位角: {p['rise_az_deg']:.0f}°")
    print(f"  最大仰角: {p['max_alt_deg']:.1f}°")
    print(f"  消失: {p['set_time'].strftime('%H:%M:%S')} UTC  "
          f"方位角: {p['set_az_deg']:.0f}°")
    print(f"  継続時間: {p['duration_min']:.1f} 分")
    print()

import ephem
import math
import numpy as np
import matplotlib.pyplot as plt
from datetime import timedelta

# ISSのTLE
iss = ephem.readtle(
    "ISS (ZARYA)",
    "1 25544U 98067A   24045.51782528  .00011186  00000-0  20208-3 0  9993",
    "2 25544  51.6412 208.6063 0005140  43.4286  70.0747 15.49562695441075"
)

# 観測地点（東京）
observer = ephem.Observer()
observer.lat = "35.6762"
observer.lon = "139.6503"
observer.elevation = 40
observer.date = "2024/2/15 00:00:00"

# 次のパスを取得
pass_info = observer.next_pass(iss)
rise_time, rise_az, max_alt_time, max_alt, set_time, set_az = pass_info

# パス期間中の仰角・方位角を10秒間隔で計算
dt_sec = 10
duration_sec = (set_time - rise_time) * 86400  # 秒に変換
n_points = int(duration_sec / dt_sec) + 1

az_list = []
el_list = []
time_list = []

for i in range(n_points):
    t = rise_time + i * dt_sec / 86400.0  # ephem.Dateの単位は日
    observer.date = t
    iss.compute(observer)
    az_list.append(math.degrees(float(iss.az)))
    el_list.append(math.degrees(float(iss.alt)))
    time_list.append(ephem.Date(t).datetime())

az_arr = np.array(az_list)
el_arr = np.array(el_list)

# ポーラープロット
fig, ax = plt.subplots(figsize=(8, 8), subplot_kw={"projection": "polar"})

# 方位角をラジアンに変換（北を上にするために調整）
az_rad = np.radians(az_arr)

# 仰角を「天頂を中心」にするため、90 - elevation をプロット
r = 90 - el_arr

# カラーマップで時間経過を表現
colors = np.linspace(0, 1, len(az_rad))
scatter = ax.scatter(az_rad, r, c=colors, cmap="cool", s=15, zorder=5)

# 開始と終了をマーク
ax.plot(az_rad[0], r[0], "^", color="lime", markersize=12, label="AOS", zorder=10)
ax.plot(az_rad[-1], r[-1], "v", color="red", markersize=12, label="LOS", zorder=10)

# 天頂を中心、地平線を外周に設定
ax.set_theta_zero_location("N")  # 北を上に
ax.set_theta_direction(-1)       # 時計回り
ax.set_ylim(0, 90)
ax.set_yticks([0, 15, 30, 45, 60, 75, 90])
ax.set_yticklabels(["90°", "75°", "60°", "45°", "30°", "15°", "0°"])
ax.set_title(f"ISS Pass - Max El: {max(el_arr):.1f}°\n"
             f"AOS: {time_list[0].strftime('%H:%M:%S')} UTC  "
             f"LOS: {time_list[-1].strftime('%H:%M:%S')} UTC",
             pad=20)
ax.legend(loc="lower right")

plt.colorbar(scatter, ax=ax, label="Time progression", shrink=0.8)
plt.tight_layout()
plt.show()

import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import numpy as np

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))

# 仰角の時間変化
ax1.plot(time_list, el_arr, color="cyan", linewidth=2)
ax1.axhline(y=10, color="yellow", linestyle="--", alpha=0.5, label="Min usable el (10°)")
ax1.fill_between(time_list, el_arr, alpha=0.2, color="cyan")
ax1.set_ylabel("Elevation [deg]")
ax1.set_title("ISS Pass - Elevation over Time")
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M:%S"))

# 方位角の時間変化
ax2.plot(time_list, az_arr, color="orange", linewidth=2)
ax2.set_ylabel("Azimuth [deg]")
ax2.set_xlabel("Time [UTC]")
ax2.set_title("ISS Pass - Azimuth over Time")
ax2.grid(True, alpha=0.3)
ax2.xaxis.set_major_formatter(mdates.DateFormatter("%H:%M:%S"))

plt.tight_layout()
plt.show()

import ephem
import math
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime, timedelta

# 複数衛星のTLEデータ（サンプル）
satellites = {
    "ISS": (
        "ISS (ZARYA)",
        "1 25544U 98067A   24045.51782528  .00011186  00000-0  20208-3 0  9993",
        "2 25544  51.6412 208.6063 0005140  43.4286  70.0747 15.49562695441075"
    ),
    "NOAA 19": (
        "NOAA 19",
        "1 33591U 09005A   24045.48691712  .00000182  00000-0  11544-3 0  9991",
        "2 33591  99.1625 39.0956 0013367 205.5765 154.4621 14.12591605781721"
    ),
    "Hubble": (
        "HST",
        "1 20580U 90037B   24045.17305220  .00001408  00000-0  73018-4 0  9991",
        "2 20580  28.4711 256.5033 0002576 110.3019 249.7904 15.09425786431641"
    ),
}

# 各衛星の地上軌跡を計算（3周回分）
start = datetime(2024, 2, 15, 0, 0, 0)
minutes = np.arange(0, 280, 0.5)

colors = {"ISS": "orange", "NOAA 19": "cyan", "Hubble": "lime"}

fig, ax = plt.subplots(figsize=(14, 7))

for name, (tle_name, line1, line2) in satellites.items():
    sat = ephem.readtle(tle_name, line1, line2)
    lats, lons = [], []

    for m in minutes:
        t = start + timedelta(minutes=float(m))
        sat.compute(t)
        lats.append(math.degrees(float(sat.sublat)))
        lons.append(math.degrees(float(sat.sublong)))

    lats = np.array(lats)
    lons = np.array(lons)

    # 経度の不連続点をNaNで処理
    lon_diff = np.abs(np.diff(lons))
    jump_idx = np.where(lon_diff > 300)[0]
    for idx in sorted(jump_idx, reverse=True):
        lats = np.insert(lats, idx + 1, np.nan)
        lons = np.insert(lons, idx + 1, np.nan)

    ax.plot(lons, lats, color=colors[name], linewidth=1.0,
            label=f"{name} (i={line2.split()[2]}°)", alpha=0.8)

ax.set_xlim(-180, 180)
ax.set_ylim(-90, 90)
ax.set_xlabel("Longitude [deg]")
ax.set_ylabel("Latitude [deg]")
ax.set_title("Multi-Satellite Ground Tracks (3 orbits)")
ax.set_aspect("equal")
ax.grid(True, alpha=0.3)
ax.legend(loc="lower left")
plt.tight_layout()
plt.show()

import ephem
import math
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime, timedelta

# NOAA 19のTLE
noaa19 = ephem.readtle(
    "NOAA 19",
    "1 33591U 09005A   24045.48691712  .00000182  00000-0  11544-3 0  9991",
    "2 33591  99.1625 39.0956 0013367 205.5765 154.4621 14.12591605781721"
)

# 1日分の地上軌跡を計算
start = datetime(2024, 2, 15, 0, 0, 0)
minutes = np.arange(0, 1440, 0.5)  # 24時間
lats, lons = [], []

for m in minutes:
    t = start + timedelta(minutes=float(m))
    noaa19.compute(t)
    lats.append(math.degrees(float(noaa19.sublat)))
    lons.append(math.degrees(float(noaa19.sublong)))

lats = np.array(lats)
lons = np.array(lons)

# 昇交点の経度を抽出（赤道を南→北に横切る点）
ascending_lons = []
for i in range(1, len(lats)):
    if lats[i-1] < 0 and lats[i] >= 0:
        ascending_lons.append(lons[i])

# 経度の不連続点をNaNで処理
lon_diff = np.abs(np.diff(lons))
jump_idx = np.where(lon_diff > 300)[0]
lats_plot = lats.copy()
lons_plot = lons.copy()
for idx in sorted(jump_idx, reverse=True):
    lats_plot = np.insert(lats_plot, idx + 1, np.nan)
    lons_plot = np.insert(lons_plot, idx + 1, np.nan)

fig, ax = plt.subplots(figsize=(14, 7))
ax.plot(lons_plot, lats_plot, color="cyan", linewidth=0.8, alpha=0.7)

# 昇交点をマーク
for asc_lon in ascending_lons:
    ax.plot(asc_lon, 0, "o", color="yellow", markersize=6)

ax.set_xlim(-180, 180)
ax.set_ylim(-90, 90)
ax.set_xlabel("Longitude [deg]")
ax.set_ylabel("Latitude [deg]")
ax.set_title("NOAA 19 Ground Track (24 hours) - Sun-Synchronous Orbit")
ax.set_aspect("equal")
ax.grid(True, alpha=0.3)
ax.axhline(y=0, color="white", linewidth=0.5, alpha=0.5)
plt.tight_layout()
plt.show()

print(f"昇交点の数: {len(ascending_lons)}")
if len(ascending_lons) >= 2:
    lon_spacings = [ascending_lons[i+1] - ascending_lons[i]
                    for i in range(len(ascending_lons)-1)]
    print(f"昇交点間の経度間隔: {np.mean(lon_spacings):.1f}° (平均)")

from sgp4.api import Satrec, WGS72
from skyfield.api import load, EarthSatellite, wgs84

# Timescaleオブジェクトを作成
ts = load.timescale()

# TLEデータ
tle_line1 = "1 25544U 98067A   24045.51782528  .00011186  00000-0  20208-3 0  9993"
tle_line2 = "2 25544  51.6412 208.6063 0005140  43.4286  70.0747 15.49562695441075"

# 衛星オブジェクトを作成
satellite = EarthSatellite(tle_line1, tle_line2, "ISS (ZARYA)", ts)

# 時刻を指定して位置を計算
t = ts.utc(2024, 2, 15, 12, 0, 0)
geocentric = satellite.at(t)

# 直下点の緯度・経度・高度
subpoint = wgs84.subpoint(geocentric)
print(f"緯度: {subpoint.latitude.degrees:.4f}°")
print(f"経度: {subpoint.longitude.degrees:.4f}°")
print(f"高度: {subpoint.elevation.km:.1f} km")

# 観測地点からの仰角・方位角
tokyo = wgs84.latlon(35.6762, 139.6503, elevation_m=40)
difference = satellite - tokyo
topocentric = difference.at(t)
alt, az, distance = topocentric.altaz()

print(f"\n東京からの観測:")
print(f"仰角:   {alt.degrees:.2f}°")
print(f"方位角: {az.degrees:.2f}°")
print(f"距離:   {distance.km:.1f} km")

from skyfield.api import load, EarthSatellite, wgs84
import numpy as np

ts = load.timescale()

tle_line1 = "1 25544U 98067A   24045.51782528  .00011186  00000-0  20208-3 0  9993"
tle_line2 = "2 25544  51.6412 208.6063 0005140  43.4286  70.0747 15.49562695441075"

satellite = EarthSatellite(tle_line1, tle_line2, "ISS (ZARYA)", ts)

# 1分間隔で1日分の時刻配列を一括で生成
minutes = np.arange(0, 1440)
t_array = ts.utc(2024, 2, 15, 0, minutes)

# 全時刻の位置を一括計算（ループ不要）
geocentric = satellite.at(t_array)
subpoints = wgs84.subpoint(geocentric)

lats = subpoints.latitude.degrees   # NumPy配列
lons = subpoints.longitude.degrees  # NumPy配列
alts = subpoints.elevation.km       # NumPy配列

print(f"計算点数: {len(lats)}")
print(f"緯度範囲: {lats.min():.2f}° ～ {lats.max():.2f}°")

import ephem
import math
from datetime import datetime

# ISSのTLE
iss = ephem.readtle(
    "ISS (ZARYA)",
    "1 25544U 98067A   24045.51782528  .00011186  00000-0  20208-3 0  9993",
    "2 25544  51.6412 208.6063 0005140  43.4286  70.0747 15.49562695441075"
)

# 観測地点（東京）
observer = ephem.Observer()
observer.lat = "35.6762"
observer.lon = "139.6503"
observer.elevation = 40

# 太陽オブジェクト
sun = ephem.Sun()

# 3日間の可視パスを探索
observer.date = "2024/2/15 00:00:00"
end_date = ephem.Date("2024/2/18 00:00:00")

visible_passes = []

while observer.date < end_date:
    try:
        pass_info = observer.next_pass(iss)
        rise_time, rise_az, max_alt_time, max_alt, set_time, set_az = pass_info

        if rise_time is None:
            break

        # 最大仰角が10度未満のパスはスキップ
        if math.degrees(float(max_alt)) < 10:
            observer.date = set_time + ephem.minute
            continue

        # パスの中間時刻で太陽の位置を確認
        mid_time = (rise_time + set_time) / 2
        observer.date = mid_time
        sun.compute(observer)
        sun_alt_deg = math.degrees(float(sun.alt))

        # ISSが影にいるか確認
        iss.compute(observer)
        is_eclipsed = iss.eclipsed

        # 可視条件: 太陽が地平線下6～18度（市民薄明～天文薄明）
        # かつ ISSが太陽に照らされている
        is_twilight = -18 < sun_alt_deg < -6

        if is_twilight and not is_eclipsed:
            visible_passes.append({
                "rise_time": ephem.Date(rise_time).datetime(),
                "set_time": ephem.Date(set_time).datetime(),
                "max_alt_deg": math.degrees(float(max_alt)),
                "sun_alt_deg": sun_alt_deg,
                "duration_min": (set_time - rise_time) * 24 * 60,
            })

        observer.date = set_time + ephem.minute
    except Exception:
        observer.date += ephem.minute * 10

print(f"3日間の可視パス数: {len(visible_passes)}\n")
for i, p in enumerate(visible_passes, 1):
    print(f"パス {i}:")
    print(f"  出現: {p['rise_time'].strftime('%Y-%m-%d %H:%M:%S')} UTC")
    print(f"  消失: {p['set_time'].strftime('%Y-%m-%d %H:%M:%S')} UTC")
    print(f"  最大仰角: {p['max_alt_deg']:.1f}°")
    print(f"  太陽高度: {p['sun_alt_deg']:.1f}°（パス中間時）")
    print(f"  継続時間: {p['duration_min']:.1f} 分")
    print()

import ephem
import urllib.request

# CelesTrakからISSの最新TLEを取得
url = "https://celestrak.org/NORAD/elements/gp.php?CATNR=25544&FORMAT=TLE"

try:
    response = urllib.request.urlopen(url)
    tle_data = response.read().decode("utf-8").strip().split("\n")

    if len(tle_data) >= 3:
        tle_name = tle_data[0].strip()
        tle_line1 = tle_data[1].strip()
        tle_line2 = tle_data[2].strip()

        iss = ephem.readtle(tle_name, tle_line1, tle_line2)
        iss.compute(ephem.now())

        print(f"衛星名: {iss.name}")
        print(f"TLEエポック: {tle_line1[18:32]}")
        print(f"現在の緯度: {iss.sublat}")
        print(f"現在の経度: {iss.sublong}")
        print(f"現在の高度: {iss.elevation/1000:.1f} km")
    else:
        print("TLEデータの取得に失敗しました")
except Exception as e:
    print(f"ネットワークエラー: {e}")

import ephem
import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from datetime import datetime, timedelta

# ISSのTLE
iss = ephem.readtle(
    "ISS (ZARYA)",
    "1 25544U 98067A   24045.51782528  .00011186  00000-0  20208-3 0  9993",
    "2 25544  51.6412 208.6063 0005140  43.4286  70.0747 15.49562695441075"
)

# 1周回分の軌道を計算
start = datetime(2024, 2, 15, 0, 0, 0)
R_earth = 6371  # km

minutes = np.arange(0, 93, 0.3)
xs, ys, zs = [], [], []

for m in minutes:
    t = start + timedelta(minutes=float(m))
    iss.compute(t)
    lat = float(iss.sublat)
    lon = float(iss.sublong)
    alt_km = iss.elevation / 1000
    r = R_earth + alt_km

    xs.append(r * math.cos(lat) * math.cos(lon))
    ys.append(r * math.cos(lat) * math.sin(lon))
    zs.append(r * math.sin(lat))

# 地球の球面を描画
u = np.linspace(0, 2 * np.pi, 50)
v = np.linspace(0, np.pi, 25)
xe = R_earth * np.outer(np.cos(u), np.sin(v))
ye = R_earth * np.outer(np.sin(u), np.sin(v))
ze = R_earth * np.outer(np.ones(np.size(u)), np.cos(v))

fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection="3d")

# 地球
ax.plot_surface(xe, ye, ze, alpha=0.15, color="blue")

# 軌道
ax.plot(xs, ys, zs, color="orange", linewidth=2, label="ISS Orbit")
ax.plot([xs[0]], [ys[0]], [zs[0]], "o", color="lime", markersize=8, label="Start")

# 赤道面のリング
theta = np.linspace(0, 2 * np.pi, 100)
ax.plot(R_earth * np.cos(theta), R_earth * np.sin(theta),
        np.zeros_like(theta), color="white", linewidth=0.5, alpha=0.5)

ax.set_xlabel("X [km]")
ax.set_ylabel("Y [km]")
ax.set_zlabel("Z [km]")
ax.set_title("ISS Orbit in 3D (1 orbit)")
ax.legend()

# アスペクト比を均等に
max_range = 7500
ax.set_xlim(-max_range, max_range)
ax.set_ylim(-max_range, max_range)
ax.set_zlim(-max_range, max_range)

plt.tight_layout()
plt.show()