Code: Select all
import numpy as np
from scipy.interpolate import CubicSpline
from datetime import datetime, timedelta
import pandas as pd
# ─── WGS84 constants ───────────────────────────────────────────────
a = 6378.137 # semi-major axis (km)
b = 6356.7523142 # semi-minor axis (km)
def wgs84_radius(z, r):
"""Local Earth radius at geocentric latitude derived from Z/r."""
sin_lat = z / r
cos_lat = np.sqrt(1 - sin_lat**2)
# WGS84 ellipsoid radius at geocentric latitude
num = (a**2 * cos_lat)**2 + (b**2 * sin_lat)**2
den = (a * cos_lat)**2 + (b * sin_lat)**2
return np.sqrt(num / den)
def altitude_above_surface(x, y, z):
"""Convert EME2000 X,Y,Z (km) to altitude above WGS84 surface (km)."""
r = np.sqrt(x**2 + y**2 + z**2)
R_earth = wgs84_radius(z, r)
return r - R_earth
# ─── Parse the OEM file ────────────────────────────────────────────
def parse_oem(filepath):
times, positions, velocities = [], [], []
with open(filepath, 'r') as f:
for line in f:
line = line.strip()
# Data lines: ISO timestamp followed by 6 floats
parts = line.split()
if len(parts) == 7:
try:
t = datetime.strptime(parts[0], '%Y-%m-%dT%H:%M:%S.%f')
x, y, z = float(parts[1]), float(parts[2]), float(parts[3])
vx, vy, vz = float(parts[4]), float(parts[5]), float(parts[6])
times.append(t)
positions.append([x, y, z])
velocities.append([vx, vy, vz])
except ValueError:
continue
return times, np.array(positions), np.array(velocities)
# ─── Main processing ───────────────────────────────────────────────
times, pos, vel = parse_oem('ISS.OEM_J2K_EPH.txt') # replace with your filename
# Convert times to seconds since epoch for interpolation
t0 = times[0]
t_sec = np.array([(t - t0).total_seconds() for t in times])
# Build cubic spline interpolators for X, Y, Z
cs_x = CubicSpline(t_sec, pos[:, 0])
cs_y = CubicSpline(t_sec, pos[:, 1])
cs_z = CubicSpline(t_sec, pos[:, 2])
# Using first three orbits starting from beginning of file
period_starts = [
datetime(2026, 7, 27, 12, 0, 0), # Period 1 start
datetime(2026, 7, 27, 13, 32, 0), # Period 2 start (~1 orbit later)
datetime(2026, 7, 27, 15, 4, 0), # Period 3 start (~2 orbits later)
]
results = []
for i, start in enumerate(period_starts):
print(f"\n{'='*60}")
print(f" PERIOD {i+1}: Starting {start.isoformat()}")
print(f"{'='*60}")
print(f"{'UTC Time':<30} {'Alt (km)':>10} {'r (km)':>10}")
print(f"{'-'*52}")
current = start
while current <= start + timedelta(minutes=90):
t_s = (current - t0).total_seconds()
if 0 <= t_s <= t_sec[-1]:
x = cs_x(t_s)
y = cs_y(t_s)
z = cs_z(t_s)
r = np.sqrt(x**2 + y**2 + z**2)
alt = altitude_above_surface(x, y, z)
print(f"{current.isoformat():<30} {alt:>10.3f} {r:>10.3f}")
results.append({
'Period': i + 1,
'UTC_Time': current.isoformat(),
'X_km': round(float(x), 3),
'Y_km': round(float(y), 3),
'Z_km': round(float(z), 3),
'Geocentric_Dist_km': round(float(r), 3),
'Alt_Above_Surface_km': round(float(alt), 3)
})
current += timedelta(minutes=3) # 3-minute intervals
df = pd.DataFrame(results)
df.to_csv('ISS_altitude_3min_three_periods.csv', index=False)
print(f"\n✅ Saved {len(df)} records to ISS_altitude_3min_three_periods.csv")