Making the axes in RA and Dec rather than 0,0

Question

I'm trying to make a figure of mars retrograde motion using Skyfield library and I successfully done it, but the axes of the figure are zeros in the middle of the plot.

I'm thinking that the issue is from center = earth.at(t).observe(mars) projection = build_stereographic_projection(center) but I don't know how to make it focus on mars but keeping the coordinates in RA (hours) and Dec (degrees). The image below is the final result of the code.

What I want is to use RA & Dec coordinates system (where 0,0 on the Sun in Vernal Equinox). Like these pictures (see the coordinates):

and this GIF:

The code:

import numpy as np
from matplotlib import pyplot as plt
from matplotlib.collections import LineCollection
from skyfield.api import Star, load
from skyfield.data import hipparcos, mpc, stellarium
from skyfield.projections import build_stereographic_projection
from skyfield.magnitudelib import planetary_magnitude


ts = load.timescale()
start_date = ts.utc(2016, 2, 1)
end_date = ts.utc(2016, 9, 2)
t_mars = ts.linspace(start_date, end_date, 20)
t = t_mars[len(t_mars) // 2]  # middle date

year, month, day, hour, minute, second = t.utc
year = year.astype(int)
month = month.astype(int)
day = day.astype(int)

eph = load('de421.bsp')
earth = eph['earth']
mars = eph['mars']

apparent = (earth).at(t).observe(mars).apparent()
m = np.array([planetary_magnitude((earth).at(t).observe(mars).apparent()) for t in t_mars])

maxmag = np.max(m)
minmag = np.min(m)

astrometric=earth.at(t).observe(mars).apparent()

ra, dec, distance = astrometric.radec()
print("RA:", ra)
print("Dec:", dec)
x, y = ra.hours, dec.degrees


size = 40 - 30 * (m - minmag) / (maxmag - minmag)


with load.open(hipparcos.URL) as f:
    stars = hipparcos.load_dataframe(f)


url = ('https://raw.githubusercontent.com/Stellarium/stellarium/master'
       '/skycultures/modern_st/constellationship.fab')

with load.open(url) as f:
    constellations = stellarium.parse_constellations(f)

edges = [edge for name, edges in constellations for edge in edges]
edges_star1 = [star1 for star1, star2 in edges]
edges_star2 = [star2 for star1, star2 in edges]

# We will center the chart on the mars's middle position.

center = earth.at(t).observe(mars)
projection = build_stereographic_projection(center)
field_of_view_degrees = 45.0
limiting_magnitude = 5.0

star_positions = earth.at(t).observe(Star.from_dataframe(stars))
stars['x'], stars['y'] = projection(star_positions)

mars_x, mars_y = projection(earth.at(t_mars).observe(mars))


bright_stars = (stars.magnitude <= limiting_magnitude)
magnitude = stars['magnitude'][bright_stars]
marker_size = (0.5 + limiting_magnitude - magnitude) ** 2.0


xy1 = stars[['x', 'y']].loc[edges_star1].values
xy2 = stars[['x', 'y']].loc[edges_star2].values
lines_xy = np.rollaxis(np.array([xy1, xy2]), 1)


fig, ax = plt.subplots(figsize=[10, 10])
fives = (np.arange(len(t_mars)) % 0.5 == 0) & (np.arange(len(t_mars)) < 30)
ax.scatter(mars_x[fives], mars_y[fives], size[fives], 'orange',
           edgecolor='black', linewidth=0.2, zorder=3)

ax.add_collection(LineCollection(lines_xy, colors=[(0.51, 0.52, 0.54, 0.3)]))


ax.scatter(stars['x'], stars['y'][bright_stars], s=marker_size, color='k', label='Stars')


offset_x, offset_y = 0.008, -0.004

for i in np.flatnonzero(fives):
    if i == 0:
        continue  

    xi = mars_x[i]
    yi = mars_y[i]
    dx = xi - mars_x[i-1]  
    dy = yi - mars_y[i-1]  
    length = np.sqrt(dx*dx + dy*dy)
    dx /= length
    dy /= length

    ax.plot([mars_x[i-1], mars_x[i]], [mars_y[i-1], mars_y[i]], c= '#0047AB', zorder=1)

for xi, yi, tstr in zip(mars_x, mars_y, t_mars.utc_strftime('%m/%d')):
    tstr = tstr.lstrip('0')
    text = ax.text(xi + offset_x, yi - offset_y, tstr, color='#0047AB',
                   ha='right', va='top', fontsize=6, weight='bold', zorder=-1, bbox=dict(facecolor='none', edgecolor='none'))
    text.set_alpha(0.5)



angle = np.pi - field_of_view_degrees / 360.0 * np.pi
limit = np.sin(angle) / (1.0 - np.cos(angle))

ax.set(
    aspect=1.0,
    xlabel='Declination (°)',
    ylabel='Right Ascension (hours)',
    xlim=(-0.8*limit, 0.9*limit),
    ylim=(-0.4*limit, 0.5*limit),
)
ax.xaxis.set_visible(True)
ax.yaxis.set_visible(True)
ax.set_aspect(1.0)


fig.savefig('Mars-Retrograde.png')