How to calculate when Mars was in retrograde around 129 BCE?

Question

I am trying to make an accurate illustration of the sky from when Ptolemy was observing the planets in Greece.

I am trying to find out when Mars was in retrograde during those years. I have read both SE questions and answer about how to calculate it mathematically but both of those are a little over my head.

Some of the silly astrology sites have the retrograde dates listed back the 1930s, I want to know what were the dates Mars was in retrograde around 129 BC.

https://physics.stackexchange.com/questions/249493/mathematically-calculate-if-a-planet-is-in-retrograde

https://physics.stackexchange.com/questions/288023/mathematically-calculate-if-a-planet-is-in-apparent-retrograde-motion

Answer 1

Mars' retrograde motion happens in the summer of 129 BC.

The JPL Ephemerides are an excellent source of solar system predictions, and are widely used. DE422 is about 600 MB and covers your period of interest.

Since your profile shows that you can program a computer, I've used the Python package Skyfield to download it, do the interpolation and calculate the apparent positions of Mars and the Sun from Alexandria in 129 BC. If I understand correctly, this would be year -128 in UTC, since there is no year zero.

You can also get similar data online using the JPL Horizons web interface.

I'll leave the problem of checking when the Sun was below the horizon but Mars was above the horizon in Alexandria on each day to you, but I've included the altitude methods in the script.

Mars is red, Sun is yellow:

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import Loader, Topos

load = Loader('~/Documents/SkyData')  # avoids multiple copies of large files
ts   = load.timescale()

data    = load('de422.bsp')
earth   = data['earth']
mars    = data['mars']
sun     = data['sun']

ts      = load.timescale()

Alexandria  = earth + Topos(latitude_degrees  = 31.2,
                            longitude_degrees = 29.9,
                            elevation_m       = 5.0)

# NOTE! The year 129 BC is year -128 in UTC I believe
# SEE   https://en.wikipedia.org/wiki/Year_zero

# Check when Mars is "up" and Sun is "down" on a given day

hours    = np.arange(0, 24, 0.1)
one_day  = ts.utc(-128, 7, 1, hours)  # July 1st

mars_obs = Alexandria.at(one_day).observe(mars)
sun_obs  = Alexandria.at(one_day).observe(sun)

mars_alt, mars_az, mars_dist = mars_obs.apparent().altaz()
sun_alt,  sun_az,  sun_dist  = sun_obs.apparent().altaz()

# Check RA and DEC in 129 BC

days     = np.arange(1, 366)
one_year = ts.utc(-128, 1, days)

mars_obs = Alexandria.at(one_year).observe(mars)
sun_obs  = Alexandria.at(one_year).observe(sun)

mars_RA, mars_Dec, mars_dist = mars_obs.apparent().radec()
sun_RA,  sun_Dec,  sun_dist  = sun_obs.apparent().radec()

if True:
    plt.figure()

    plt.subplot(3, 1, 1)

    marsalt, sunalt = mars_alt.degrees, sun_alt.degrees

    plt.plot(hours, marsalt, '-r')
    plt.plot(hours, sunalt,  '-y')
    plt.plot(hours, np.zeros_like(hours), '-k')
    plt.title('altitude of Mars and Sun from Alexandria, July 1, 129 BC')
    plt.xlabel('hours (UTC)')
    plt.ylabel('altitude (degs)')
    plt.ylim(-10, None)
    plt.xlim(0, 24)

    marsra, marsdec = mars_RA._degrees, mars_Dec._degrees
    sunra,  sundec  = sun_RA._degrees,  sun_Dec._degrees

    # get rid of ugly lines
    dsunra, dmarsra = sunra[1:] - sunra[:-1], marsra[1:] - marsra[:-1]
    sunra, marsra = sunra[:-1], marsra[:-1]
    sunra[dsunra<-100]   = np.nan
    marsra[dmarsra<-100] = np.nan

    plt.subplot(3, 1, 2)
    plt.plot(days, mars_RA._degrees, '-r')
    plt.plot(days, sun_RA._degrees,  '-y')
    plt.title('RA of Mars and Sun from Alexandria, year of 129 BC')
    plt.xlabel('days')
    plt.ylabel('RA (degs)')
    plt.xlim(0, 365)

    plt.subplot(3, 1, 3)
    plt.plot(mars_RA._degrees, mars_Dec.degrees, '-r')
    plt.plot(sun_RA._degrees,  sun_Dec.degrees,  '-y')
    plt.plot([0, 360],  [0, 0], '-k')
    plt.title('RA vs Dec of Mars and Sun from Alexandria, year of 129 BC')
    plt.xlabel('RA (degs)')
    plt.ylabel('Dec (degs)')
    plt.xlim(0, 360)

    plt.show()