I wrote some code to animate a satellite/body orbiting around Earth using Matplotlib. However, Matplotlib either sets the orbit behind or doesn’t display the orbit hierarchy layer correctly, despite setting different zorder values for each ax figure, as you can see in the images below. So instead, I decided to use Pyvista to plot the orbit in 3D. However, I can’t find the right way to animate the satellite orbiting the cental body and the rotation of it like in Matplotlib because I can’t find documentation on how to do it. Any idea how to animate it?

Code:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import pyvista as pv
import matplotlib.animation as animation
import matplotlib.patches as patches

# Constants
M = 5.974 * 10 ** 24  # mass of the earth
G = 6.67430 * 10 ** (-20) #6.67259 * 10 ** (-20)  # gravitational constant
R = 6378.0  # radius of the moon in km 
mu = M *G ## 3.986 * 10 ** 5  #  gravitational parameter of earth km^3/sec^2 
# Orbital elements
a = 6789 # semi-major axis in km of ISS
e = 0.003280    # eccentricity of ISS
i = np.deg2rad(51.6)  # inclination
raan = np.deg2rad(72.5085)  # right ascending node
w = np.deg2rad(202.3397)  # argument of perigee
p = a * (1 - e**2)


# Initial conditions for numerical calculations
 
nu= np.deg2rad(0)  # true anomaly Is possible to make a initial value by calculatin the mean anomaly
t0 = 0           # initial time 
tf = 2 * np.pi * np.sqrt ( a **3 / mu) # final time = 2 orbit 
dt = 20           # time step 
#n = np.sqrt( mu / a ** 3)
print(tf)
# Rotation Matrix from PQW to ECI
def Q_x(raan, w, i):
    q = np.array([
        [-np.sin(raan) * np.cos(i) * np.sin(w) + np.cos(raan) * np.cos(w), -np.sin(raan) * np.cos(i) * np.cos(w) - np.cos(raan) * np.sin(w), np.sin(raan) * np.sin(i)],
        [np.cos(raan) * np.cos(i) * np.sin(w) + np.sin(raan) * np.cos(w), np.cos(raan) * np.cos(i) * np.cos(w) - np.sin(raan) * np.sin(w), -np.cos(raan) * np.sin(i)],
        [np.sin(i) * np.sin(w), np.sin(i) * np.cos(w), np.cos(i)]
    ])
    return q
# Building the position and velocity vectors functions using classic kepler elements
def r_po(p,e,nu):
    cos_v = np.cos(nu)
    sin_v = np.sin(nu)
    #r = (h**2 / (mu)) / (1 + e * cos_v)
    r = p / (1 + e * cos_v)
    perifocal_r = np.array([r * cos_v, r * sin_v, 0])
    return perifocal_r

def v_po(p,e,nu):
     x = - np.sqrt(mu/p) * np.sin(nu)
     y = np.sqrt(mu/p) * (e + np.cos(nu))
     z = 0
     perifocal_vel = np.array([x,y,z])
     return perifocal_vel
# Rotation from PQW to ECI
def p_perifocal_to_ECI(p, e, nu, w , i, raan):
    r_0 = r_po(p,e,nu)
    r_1 = Q_x(raan, w,i).dot(r_0)
    return r_1

def v_perifocal_to_ECI(p, e, nu, w , i, raan):
    v_0 = v_po(p,e,nu)
    v_1 = Q_x(raan,w,i).dot(v_0)
    return v_1
# We store the initial values
x0 = p_perifocal_to_ECI(p, e, nu, w , i, raan)[0]
y0 = p_perifocal_to_ECI(p, e, nu, w , i, raan)[1]
z0 = p_perifocal_to_ECI(p, e, nu, w , i, raan)[2]

vx0 = v_perifocal_to_ECI(p, e, nu, w , i, raan)[0]
vy0 = v_perifocal_to_ECI(p, e, nu, w , i, raan)[1]
vz0 = v_perifocal_to_ECI(p, e, nu, w , i, raan)[2]

# Define function for RK4 method
def f(t, propagation):
    x, y, z, vx, vy, vz = propagation # Assigns each value to each array 
    r = np.sqrt(x**2 + y**2 + z**2)
    ax = -G * M * x / r**3
    ay = -G * M * y / r**3
    az = -G * M * z / r**3
    return np.array([vx, vy, vz, ax, ay, az])


# Initialize arrays to store results, zeros_like fills the the array with zeros
t_values = np.arange(t0, tf, dt)
X = np.zeros_like(t_values)
Y = np.zeros_like(t_values)
Z = np.zeros_like(t_values)
VX = np.zeros_like(t_values)
VY = np.zeros_like(t_values)
VZ= np.zeros_like(t_values)
p_vector = []
v_vector = []
# Set initial values from perifocal to ECI
X[0] = x0
Y[0] = y0
Z[0] = z0
VX[0] = vx0
VY[0] = vy0
VZ[0] = vz0

# Apply Runge Kutta order 4 
for i in range(1, len(t_values)): #  Iterates over each time step index
    t = t_values[i-1]
    propagation = [X[i-1], Y[i-1], Z[i-1], VX[i-1], VY[i-1], VZ[i-1]] # we're retrieving the time value at the previous time step in the t_values array
    k1 = dt * f(t, propagation)
    k2 = dt * f(t + 0.5*dt, propagation + 0.5*k1)
    k3 = dt * f(t + 0.5*dt, propagation + 0.5*k2)
    k4 = dt * f(t + dt, propagation + k3)
    propagation += (k1 + 2*k2 + 2*k3 + k4) / 6 # propagation_+1 = propagation + 
    X[i], Y[i], Z[i], VX[i], VY[i], VZ[i] = propagation
    p_vector.append([X[i], Y[i], Z[i]])
    v_vector.append([VX[i], VY[i], VZ[i]])
# Transponse the vectors to plot them
position = np.array(p_vector).T 
v_vector = np.array(v_vector).T

# The plot with Matplotlib or Pyvista starts here
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Plot of the orbit
line, = ax.plot(position[0], position[1], position[2], zorder = 3)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
#path, = ax.plot([], [],[], color = 'black',  marker='o', markersize=8, zorder = 2)

# Plot Earth sphere
phi = np.linspace(0, 2 * np.pi, 100)
theta = np.linspace(0, np.pi, 100)
xm = R * np.outer(np.cos(phi), np.sin(theta))
ym = R * np.outer(np.sin(phi), np.sin(theta))
zm = R * np.outer(np.ones(np.size(phi)), np.cos(theta))
Earth = ax.plot_surface(xm, ym , zm , rstride=4, cstride=4, facecolor="grey", edgecolor="black", zorder = 2)
#Euler roation about Z Axis
def rotation_matrix(angle):
    c, s = np.cos(angle), np.sin(angle)
    array= np.array([[c, -s, 0],
                     [s, c, 0],
                     [0, 0, 1]])
    return array
# 3D plot of a satellite
Rm = 500
rho = np.linspace(0, 2 * np.pi, 100)
gamma = np.linspace(0, np.pi, 100)
x_o = Rm * np.outer(np.cos(rho), np.sin(gamma))
y_o = Rm * np.outer(np.sin(rho), np.sin(gamma))
z_o = Rm * np.outer(np.ones(np.size(rho)), np.cos(gamma))
sat = ax.plot_surface(x_o,y_o,z_o, rstride=4, cstride=4, facecolor="b",edgecolor="black", zorder = 2)
# Sat movement along the orbit path
def sat_move(position, i):
    global sat
    sat.remove()
    sat_x = x_o + position[0, i]
    sat_y = y_o + position[1, i]
    sat_z = z_o + position[2, i]
    sat = ax.plot_surface(sat_x, sat_y, sat_z, rstride=4, cstride=4, facecolor="b",edgecolor="black", zorder=10)
# Planet/Body rotation
def Earth_move(i):
    global Earth
    Earth.remove()
    angle = np.radians(i)
    R = rotation_matrix(angle)
    
    x_rot = R[0, 0] * xm + R[0, 1] * ym + R[0, 2] * zm
    y_rot = R[1, 0] * xm + R[1, 1] * ym + R[1, 2] * zm
    z_rot = R[2, 0] * xm + R[2, 1] * ym + R[2, 2] * zm
    Earth = ax.plot_surface(x_rot , y_rot  , z_rot  , rstride=4, cstride=4, facecolor="grey", edgecolor="black", zorder = 2)

def orbit(i):
    sat_move(position, i)
    Earth_move(i)
    return  sat, Earth #Earth
ax.set_aspect('equal') 
ani = animation.FuncAnimation(fig , orbit,  frames = len(t_values), interval=1)
plt.show()

Matplotlib:

With:

line, = ax.plot(position[0], position[1], position[2], zorder = 4)

Pyvista, but I don’t know how to animate the satellite orbiting the central body:

plot = pv.Plotter()

# Plot orbit

points = np.column_stack((X, Y,Z))
orbit_path = pv.PolyData(points)


# Plot Moon
phi = np.linspace(0, 2 * np.pi, 100)
theta = np.linspace(0, np.pi, 100)
xm = R * np.outer(np.cos(phi), np.sin(theta))
ym = R * np.outer(np.sin(phi), np.sin(theta))
zm = R * np.outer(np.ones(np.size(phi)), np.cos(theta))
moon = pv.StructuredGrid(xm, ym, zm)
plot.add_mesh(moon, color="grey", smooth_shading = True)

# Correct the aliasing
plot.enable_anti_aliasing()




plot.add_mesh(orbit_path, color="red", label = 'Orbit', style = 'points', point_size = 5)
plot.add_legend(bcolor = None, size = (0.1,0.1), face = None , loc = 'upper right')

plot.show_grid()
plot.show_axes()
plot.add_floor('-z', lighting=True, color='b', pad=1.0, opacity = 0.5) # Pad sets the size much better

# Show the plot plot.show()
plot.show()