satellite-simulator/visualizer.py at main · JamesDCowley/satellite-simulator · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# 1. Load your satellite position data (Space-separated)
# The file is expected to have 3 columns: X Y Z (in meters)
data = np.loadtxt('positions.csv')
x_sat, y_sat, z_sat = data[:, 0], data[:, 1], data[:, 2]

# 1b. Load rotation data (three angles per time step in radians)
rot = np.loadtxt('rotations.csv')
rot_x, rot_y, rot_z = rot[:, 0], rot[:, 1], rot[:, 2]

# 2. Setup Figure and 3D Axes
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')

# 3. Create the Earth (a sphere at the origin)
earth_radius = 6371000
u, v = np.mgrid[0:2*np.pi:30j, 0:np.pi:20j]
x_earth = earth_radius * np.cos(u) * np.sin(v)
y_earth = earth_radius * np.sin(u) * np.sin(v)
z_earth = earth_radius * np.cos(v)
ax.plot_surface(x_earth, y_earth, z_earth, color='royalblue', alpha=0.6, edgecolor='w', linewidth=0.2)

# 4. Create satellite objects (point, orbit trail, and orientation arrow)
sat_point, = ax.plot([], [], [], 'ro', markersize=5, label='Satellite')
orbit_trail, = ax.plot([], [], [], 'y-', alpha=0.7, linewidth=1, label='Orbit Path')
orient_arrow, = ax.plot([], [], [], 'g-', linewidth=2, label='Top-face direction')
angle_text = ax.text2D(0.03, 0.95, '', transform=ax.transAxes)
time_text = ax.text2D(0.03, 0.90, '', transform=ax.transAxes)

# 5. Final Plot Styling
limit = max(np.max(np.abs(data)), earth_radius * 1.5)
ax.set_xlim(-limit, limit)
ax.set_ylim(-limit, limit)
ax.set_zlim(-limit, limit)
ax.set_title("Satellite Orbit Animation")
ax.legend()

# 6. Animation Functions
def init():
    sat_point.set_data([], [])
    sat_point.set_3d_properties([])
    orbit_trail.set_data([], [])
    orbit_trail.set_3d_properties([])
    orient_arrow.set_data([], [])
    orient_arrow.set_3d_properties([])
    angle_text.set_text('')
    time_text.set_text('')
    return sat_point, orbit_trail, orient_arrow, angle_text, time_text

def euler_to_top_direction(rx, ry, rz):
    # rotation vector is (x, y, z) in radians; assume order roll(X), pitch(Y), yaw(Z)
    cx, sx = np.cos(rx), np.sin(rx)
    cy, sy = np.cos(ry), np.sin(ry)
    cz, sz = np.cos(rz), np.sin(rz)

    # rotation matrices (intrinsic XYZ)
    Rx = np.array([[1, 0, 0], [0, cx, -sx], [0, sx, cx]])
    Ry = np.array([[cy, 0, sy], [0, 1, 0], [-sy, 0, cy]])
    Rz = np.array([[cz, -sz, 0], [sz, cz, 0], [0, 0, 1]])

    R = Rz @ Ry @ Rx
    top_local = np.array([0.0, 0.0, 1.0])
    return R @ top_local


def update(frame):
    pos = np.array([x_sat[frame], y_sat[frame], z_sat[frame]])

    # current position
    sat_point.set_data([pos[0]], [pos[1]])
    sat_point.set_3d_properties([pos[2]])

    # trailing line (shows path traveled so far)
    orbit_trail.set_data(x_sat[:frame], y_sat[:frame])
    orbit_trail.set_3d_properties(z_sat[:frame])

    # orientation arrow on top face
    dir_top = euler_to_top_direction(rot_x[frame], rot_y[frame], rot_z[frame])
    arrow_length = earth_radius * 0.1
    arrow_end = pos + dir_top * arrow_length

    orient_arrow.set_data([pos[0], arrow_end[0]], [pos[1], arrow_end[1]])
    orient_arrow.set_3d_properties([pos[2], arrow_end[2]])

    angle_text.set_text(f"Orientation (rad): ({rot_x[frame]:.2f}, {rot_y[frame]:.2f}, {rot_z[frame]:.2f})")

    # elapsed time
    total_seconds = int(frame * 0.1)  # set to dt (seconds per frame) if different
    hours = total_seconds // 3600
    minutes = (total_seconds % 3600) // 60
    seconds = total_seconds % 60
    time_text.set_text(f"T+{hours:02d}:{minutes:02d}:{seconds:02d}")

    return sat_point, orbit_trail, orient_arrow, angle_text, time_text

# Increase interval for slower motion; decrease for faster
ani = FuncAnimation(fig, update, frames=len(x_sat), init_func=init,
                    interval=100, blit=True)

ax.set_box_aspect([1,1,1])
plt.show()