tl;dr: use a parametric equation.
If the Earth were not rotating, then we would have something like
\begin{align}
x & = \cos \omega (t-t_0)\\
y & = \sin \omega (t-t_0) \ \cos i\\
z & = \sin \omega (t-t_0) \ \sin i\\
\end{align}
where the radius of the orbit is 1, $\omega$ is $2 \pi/T$ and $T$ is the orbital period, and $i$ is the inclination of the orbit.
Then we would have
\begin{align}
lon & = \arctan2(y, x) + const\\
lat & = \arcsin(z)\\
\end{align}
If the Earth is rotating then
$$lon = \arctan2(y, x) - \omega_E (t-t_0) + const$$
where $\omega_E$ is $2 \pi/T_D$ and $T_D$ is a sidereal day (23h, 56m, 4s roughly).
Solving this for longitude as a function of latitude looks like some serious work and I am not sure there is an analytical solution.
Instead you can try the parametric equation approach where you first make a hidden table of times, and then solve for $lon(t)$ and $lat(t)$ and plot $lat$ vs $lon$
Here is a plot, I haven't adjusted $t_0$ or $const$ and just used rough values for $\omega$, $\omega_E$ and $i$ but it should be enough to get you stared.
$t_0$ and $const$ represent the known starting conditions of the spacecraft that you are trying to plot; $t_0$ is the time at which it crosses the equator going north, and $const$ is the longitude on the Earth below the spacecraft at that time.
Here is some further reading:
Python script:
import numpy as np
import matplotlib.pyplot as plt
twopi = 2*np.pi
to_degs, to_rads = 180/np.pi, np.pi/180.
omega = twopi/(9260)
omega_E = twopi/(233600 + 56*60 + 4)
time = 60 * np.arange(101.) # 100 minutes
t0 = 1000. # arbitrary, you can fit this later
inc = 51.
const = 1.0 # arbitrary, you can fit this later
x = np.cos(omega * (time-t0))
y = np.sin(omega * (time-t0)) * np.cos(to_radsinc)
z = np.sin(omega (time-t0)) * np.sin(to_rads*inc)
lon = np.arctan2(y, x) - omega_E * (time-t0) + const
lat = np.arcsin(z)
if True:
plt.figure()
plt.plot(to_degslon, to_degslat, '.k')
plt.xlim(-180, 180)
plt.ylim(-60, 60)
#plt.gca().set_aspect('equal')
plt.show()
