Plotting ground Tracks with python only using trigonometry

I'm trying to plot a ground track in Python from a TLE file without using explicit astropy or other packages. I've encountered a problem with the code that I can't figure out. The issue pertains to the continuity of the ground track path. The code appears to cut the path and then continues on the other side, creating a cross-like pattern. I think the problem has to do with the calculation of the longitude but I'm not sure. I leave my code and the plot below

'''

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()

'''

Not sure why you wouldn't want to use AstroPy or https://rhodesmill.org/skyfield/, but OK, that's your ground rules. :) I'm assuming you aren't too worried about accuracy since you are doing basic 2-body propagation of mean orbital elements.

These lines do look a bit suspect:

# 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

alpha is using arctan, not arctan2, so it will only return an angle between -pi/2 and pi/2. I would see if you can rewrite your equation using arctan2 to preserve the quadrant information.