仅用三角函数在Python中绘制卫星地面轨迹的连续性问题求助

大家好，我最近在尝试不依赖Astropy这类专用天文库，纯靠基础Python工具从TLE文件计算并绘制卫星的地面轨迹，但遇到了一个棘手的问题——轨迹的连续性出了问题：绘制出来的路径会突然被截断，然后直接跳到地球的另一侧继续，形成了类似十字的奇怪图案。我怀疑问题出在经度的计算逻辑上，但反复检查也没找到具体的错误点。

下面是我的完整代码：

import matplotlib.pyplot as plt
import numpy as np
from astropy.time import Time

#TLE data for the satellite "Beesat9"
#BEESAT 9
tle_line1 = "1 44412U 19038AC  23251.87820378  .00033749  00000+0  10095-2 0  9994"
tle_line2 = "2 44412  97.6799 235.5575 0014498  22.8588 337.3296 15.34383582230971"

# To get the individual values
tle_elements1 = tle_line1.split()
tle_elements2 = tle_line2.split()

# Extract the orbital elements values from tle_elements
m =float(tle_elements2[6]) #°
i = np.radians(float(tle_elements2[2])) # rad
w = np.radians(float(tle_elements2[5])) # rad
Omega = np.radians(float(tle_elements2[3])) # rad
e = float("0." + tle_elements2[4])
M_mo=round(float(tle_elements2[7]), 8) # rev/day
Epoch = str(float(tle_elements1[3]))[2:] # fraction day format xxx.xxxx
Epoch = float(Epoch)

# Get Julian Day
today = Time.now()
jd = today.jd

#Julian Centuries converter
def JC_converter(jd):
    t = (jd - 2451545.0) / 36525
    return t

jc = JC_converter(jd)

# Siderial time for Greenwich
def greenW(jc):
    Theta= 1.753368559 + 628.3319707 * jc + 0.0000067707 * jc**2
    n = Theta / (2*np.pi)
    n = np.floor(n)
    Theta = Theta - n*2*np.pi
    return Theta

Theta_k = greenW(jc)

# Next pass
w_sid= (2*np.pi) / 86164.0916 # angular siderial velocity
t_eqnx = (2*np.pi - Theta_k) / (w_sid) #in w_siderial seconds

# Transforms epoch to seconds rounding to eigth decimal since midnight
def utc(epoch):
    n_epoch = round(epoch - int(epoch),8)
    #days to sec
    seconds = n_epoch*86400
    return seconds

epoch_sec = utc(Epoch) # 75876.806592 sec

# Calculate the satellite's position at each time point
latitude = []
longitude = []

# M_0 from M = M_0 + n(t - tp), M_0 = mean_prime
mean_prime = np.radians(m)
motion_anomaly = M_mo / 86400 # to rev/sec

# Time of Periapsis passage
tp = epoch_sec - (mean_prime / motion_anomaly)

# Orbital Period T of the cubesat
orbital_period = (2*np.pi / M_mo )*86400 # 1 orbital period 86400 seconds in a day
time_elapsed = 0

for time in range(0, int(orbital_period),10):
    # Calculate mean anomaly propagation
    mean_anomaly = mean_prime + motion_anomaly * (time - tp)

    # Calculate eccentric anomaly using Newton's method
    E = mean_anomaly  # Initial guess for Eccentric Anomaly
    tolerance = 1e-8
    max_iterations = 1000
    for _ in range(max_iterations):
        f_e = E - e * np.sin(E) - mean_anomaly
        f_prime_e = 1.0 - e * np.cos(E)
        E_new = E - f_e / f_prime_e
        if abs(E_new - E) < tolerance:
            break
        E = E_new

    # Calculation new Epoch and time of the position of the cubesat
    t_new =  tp + ((E - (e*np.sin(E))) /motion_anomaly)

    #True Anomaly
    true_anomaly = 2.0 * np.arctan2(
        np.sqrt(1.0 + e) * np.sin(E / 2.0),
        np.sqrt(1.0 - e) * np.cos(E / 2.0))

    # Calculate satellite's latitude
    satellite_latitude = np.arcsin(np.sin(i) * (np.sin(w + true_anomaly)))
    satellite_latitude = np.degrees(satellite_latitude) # Converting to degrees

    # Calculate satellite's longitude
    alpha = np.arctan(np.cos(i) * np.tan(w + true_anomaly)) + Omega
    satellite_longitude = alpha - w_sid*(t_new - t_eqnx)
    satellite_longitude = np.degrees(satellite_longitude) # Converting to degrees

    # Append latitude and longitude to the lists
    latitude.append(satellite_latitude)
    longitude.append(satellite_longitude)

    # Update time elapsed
    time_elapsed += 10

    # one complete orbit has been completed
    if time_elapsed >= orbital_period:
        break

# Load background image
background_image_path = r'D:\Satellite Code\earth.jpg'
background_img = plt.imread(background_image_path)

# Create the plot
plt.figure(figsize=(15.2, 8.2))
plt.imshow(background_img, extent=[-180, 180, -90, 90])

# Plot the ground track
plt.scatter(longitude, latitude, label="BEESAT 9 Ground Track", color='red', marker='o', s=4)
plt.xlabel("Longitude (degrees)")
plt.ylabel("Latitude (degrees)")
plt.title("BEESAT 9 Ground Track ")

# Show the plot
plt.legend()
plt.grid(True, color='w', linestyle=":", alpha=0.4)
plt.show()

目前我猜测问题可能出在经度计算的象限处理或者经度值的范围归一化上，比如np.arctan的返回范围有限，没法正确覆盖所有象限，或者计算后的经度没有处理成-180到180的范围，导致出现跳变。有没有大佬能帮我看看哪里出问题了？

备注：内容来源于stack exchange，提问作者ed190