Mathematically calculate if a Planet is in Retrograde I am trying to find out if a planet is in an apparent retrograde motion with respect to the earth at any given point in time. Given the time in Julian days.

I already have the geocentric and heliocentric coordinates and velocities of the planet using the JPL ephemeris, but not sure of what mathematical formula to use to identify the retrograde motion.

I read this article "https://physics.stackexchange.com/questions/249493/mathematically-calculate-if-a-planet-is-in-retrograde" which had a similar question but I am not sure of:


Is the answer correct as it has not been marked as correct?
How to represent earth in the XY plane only when it has XYZ planes and why do we even need to do that? Can't we just use the heliocentric equatorial plane which I think is the same as the one used by NASA JPL ephemeris.
Here is a quick answer. Apparent retrograde motion of a planet as seen from Earth means that its ecliptic longitude is decreasing. (Thanks to @barrycarter for correcting my previous statement.) You have to calculate the apparent angular position of the planet from the Earth in ecliptic coordinates (as explained here and here) and find places where the ecliptic longitude is decreasing. 

You can subtract the positions of Earth from the position of the planet in either coordinates. Then you express that vector in spherical coordinates. It doesn't matter if it's heliocentric or geocentric. What matters is the direction that the axes of the coordinate system is pointing, not where the center is. 

If you want RA and dec, use ICRF coordinates, and these can be transformed to ecliptic coordinates as well. An easy way to do this is to use the python package Skyfield.

Here's a sample script and a plot. This is Mars from 2010.0 to 2020.0 and the thick red line represents retrograde motion.



import numpy as np
import matplotlib.pyplot as plt

from skyfield.api import load

planets = load('de421.bsp')

earth   = planets['earth']
mars    = planets['mars']

year_zero = 2010
days = np.linspace(1, 3650, 10000)
years = year_zero + days / 365.2564
ts = load.timescale()
t = ts.utc(year_zero, 1, days)

# thanks to @barrycarter's comment, do it the right way!
eclat, eclon, ecd = earth.at(t).observe(mars).ecliptic_latlon()

eclondgs      = (180./np.pi) * eclon.radians
eclondel      = eclondgs[1:] - eclondgs[:-1]

eclondel[eclondel &lt; -300] += 360. # this is a fudge for now
eclondel[eclondel &gt; +300] -= 360. # this is a fudge for now

prograde   = eclondel &gt; 0.

eclon_prograde   = eclondgs.copy()[:-1]
eclon_retrograde = eclondgs.copy()[:-1]

eclon_prograde[-prograde]  = np.nan
eclon_retrograde[prograde] = np.nan

plt.figure()
plt.plot(years[:-1], eclon_prograde,   '-g', linewidth=1)
plt.plot(years[:-1], eclon_retrograde, '-r', linewidth=3)
plt.ylim(0, 360)
plt.title('Mars Ecliptic longitude (degrees) AD 2010.0 to 2020.0',
          fontsize=16)
plt.show()