When plotting the orbit of an object, the 3D plot doesn’t hide the lines behind the sphere. How can I fix this? I tried using zorder in matplotlib, but it either shows the entire orbit or doesn’t display the orbit in front of the sphere.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Constants
re = 6378.0  # radius of the earth in km
mu = 3.986 * 10 ** 5  # gravitational parameter of earth km^3/sec^2


# Kepler elements
a = 6792 # semi-major axis in km of ISS
e = 0.0003612  # eccentricity of ISS
i = np.deg2rad(51.6367)  # inclination
Omega = np.deg2rad(198.1689 )  # right ascending node
w = np.deg2rad(120.7165)  # argument of perigee
v_anomaly = []  # True anomaly (empty list, assuming it will be filled later)
n = np.sqrt(mu / a**3)  # mean motion in radians per second
p = a * (1 - e**2)  # semi-parameter p
m_0 = 0  # Initial value of mean anomaly
t_p = 0  # seconds

# Period
orbital_period = 2 * np.pi / n # One orbit of the moon
# Vectors
r_position = []
v_position = []

def r_po(p,e,v):
    cos_v = np.cos(v)
    sin_v = np.sin(v)
    #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,v):
     x = - np.sqrt(mu/p) * np.sin(v)
     y = np.sqrt(mu/p) * (e + np.cos(v))
     z = 0
     #velo = np.sqrt((mu/r)(2 - (1 - (e**2) / 1 + e*np.cos(v_anomaly))))
     perifocal_vel = np.array([x,y,z])
     return perifocal_vel
def Q_x(Omega, w, i):
    q = np.array([
        [-np.sin(Omega) * np.cos(i) * np.sin(w) + np.cos(Omega) * np.cos(w), -np.sin(Omega) * np.cos(i) * np.cos(w) - np.cos(Omega) * np.sin(w), np.sin(Omega) * np.sin(i)],
        [np.cos(Omega) * np.cos(i) * np.sin(w) + np.sin(Omega) * np.cos(w), np.cos(Omega) * np.cos(i) * np.cos(w) - np.sin(Omega) * np.sin(w), -np.cos(Omega) * np.sin(i)],
        [np.sin(i) * np.sin(w), np.sin(i) * np.cos(w), np.cos(i)]
    ])
    return q
# Propagation of the true anomaly

for t_op in range (int(orbital_period)):
    mean =  n*(t_op - t_p)
    e_anom = mean  # Initial guess for Eccentric Anomaly
    tolerance = 1e-8
    max_iterations = 1000

    for _ in range(max_iterations):
           f_e = e_anom - (e * np.sin(e_anom) - mean)
           f_prime_e = 1 - (e * np.cos(e_anom))

           E_new = (e_anom - f_e) / f_prime_e
           # function returns the absolute value of the given number
           if abs(E_new - e_anom) < tolerance:
               break

           e_anom = E_new
           # V_0        
    true_0 = 2.0 * np.arctan2(
                np.sqrt(1.0 + e) * np.sin(e_anom / 2.0),
                 np.sqrt(1.0 - e) * np.cos(e_anom / 2.0))
           #keep true anomaly in range
    v_0 = np.degrees(true_0) % 360
    v_anomaly.append(v_0)



full = 360 
for v in v_anomaly:
     anomaly_n= (int(v) + 180) % full - 180
     tru_ano = np.radians(anomaly_n)
    
     q_0 = Q_x(Omega, w, i)
     r_0= r_po(p, e, tru_ano)
     r_1 = q_0.dot(r_0)
     r_position.append(r_1)
r_position = np.array(r_position).T  # Transpose the r_position array for correct shape

X = r_position[0]
Y = r_position[1]
Z = r_position[2]

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Plot orbit
ax.plot(X, Y, Z, color='red', zorder = 3)

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('Orbit around Earth')


# Plot Earth
phi = np.linspace(0, 2 * np.pi, 100)
theta = np.linspace(0, np.pi, 100)
xm = re * np.outer(np.cos(phi), np.sin(theta))
ym = re * np.outer(np.sin(phi), np.sin(theta))
zm = re * np.outer(np.ones(np.size(phi)), np.cos(theta))
ax.plot_surface(xm, ym, zm, rstride=4, cstride=4, facecolor="grey", edgecolor="black", zorder = 1)
ax.set_aspect('equal')
plt.show()

For orbit: zorder = 3 and Earth = 1:

For orbit: zorder = 2.5 and Earth = 1: