I have a script which converts Keplerian orbital elements to cartesian state vectors. Units are Years, AU and radians.

Here is my code:

def f_mean_anomaly(time_since_perihelion, semi_major_axis):
    orbital_period = math.sqrt(semi_major_axis**3) 
    mean_anomaly = (2 * math.pi * time_since_perihelion) / orbital_period
    return mean_anomaly

def f_apoapsis(semi_major_axis, eccentricity):
    return semi_major_axis * (1 + eccentricity)

def f_periapsis(semi_major_axis, eccentricity):
    return semi_major_axis * (1 - eccentricity)

def true_anomaly(mean_anomaly, eccentricity):
    mean_anomaly -= math.floor(mean_anomaly / (2 * math.pi)) * (2 * math.pi)

    if eccentricity > 0.5:
        eccentric_anomaly = math.pi
    else:
        eccentric_anomaly = mean_anomaly

    while True:
        delta = (mean_anomaly - (eccentric_anomaly - eccentricity * math.sin(eccentric_anomaly))) / (1.0 - eccentricity * math.cos(eccentric_anomaly))
        eccentric_anomaly += delta
        if abs(delta) <= 1e-15 * eccentric_anomaly:
            break

    true_anomaly = 2.0 * math.atan(math.sqrt((1.0 - eccentricity) / (1.0 + eccentricity)) * math.tan(eccentric_anomaly / 2))
    if true_anomaly < 0.0:
        true_anomaly += (2 * math.pi)

    return true_anomaly

def eccentric_anomaly(mean_anomaly, eccentricity):
    mean_anomaly -= math.floor(mean_anomaly / (2 * math.pi)) * (2 * math.pi)

    if eccentricity > 0.5:
        eccentric_anomaly = math.pi
    else:
        eccentric_anomaly = mean_anomaly

    while True:
        delta = (mean_anomaly - (eccentric_anomaly - eccentricity * math.sin(eccentric_anomaly))) / (1.0 - eccentricity * math.cos(eccentric_anomaly))
        eccentric_anomaly += delta
        if abs(delta) <= 1e-15 * eccentric_anomaly:
            break

    nu = 2.0 * math.atan(math.sqrt((1.0 + eccentricity) / (1.0 - eccentricity)) * math.tan(0.5 * eccentric_anomaly))
    if nu < 0.0:
        nu += (2 * math.pi)

    return eccentric_anomaly

def kepler_to_cartesian(true_anomaly, semi_major_axis, eccentricity, longitude_of_ascending_node, argument_periapsis, inclination):
    phi = true_anomaly
    a = semi_major_axis
    e = eccentricity
    omega = longitude_of_ascending_node
    omega_prime = argument_periapsis
    i = inclination

    radius = a * ((1 - e**2) / (1 + e * math.cos(phi)))
    r = radius

    x_value = r * (math.cos(omega) * math.cos(omega_prime + phi) - math.sin(omega) * math.sin(omega_prime + phi) * math.cos(i))
    y_value = r * (math.sin(omega) * math.cos(omega_prime + phi) + math.cos(omega) * math.sin(omega_prime + phi) * math.cos(i))
    z_value = r * (math.sin(omega_prime + phi) * math.sin(i))

    return x_value, y_value, z_value

and this is what I am trying to reproduce using the table below:

When I plot the RA using RA = x_value / (0.04 * parsec) # conversion from RA to x in the galactic center. This is what I get for

data = { 2023: (717.1470698197126, 3877.424985077577, 3755.4666153429635), 2021: (-4621.544894435808, -1177.926886174313, 3665.413967770432), 2020: (-3740.6813990975033, -3035.6632272995867, 529.702793927504), 2018: (62.25366882760013, -2147.079670435817, -2580.91000905229), 2015: (2876.854950728359, 2638.4632388430355, -51.787629080512914), 2012: (-4010.5553090894045, 480.7554395841735, 4939.903684859025), }

Is there something wrong with my conversion because this is nowhere close to what I am going for?

References used: https://downloads.rene-schwarz.com/download/M001-Keplerian_Orbit_Elements_to_Cartesian_State_Vectors.pdf, space.stackexchange.com/questions/60006/, space.stackexchange.com/a/19335/48051, space.stackexchange.com/a/8915/48051, space.stackexchange.com/a/59600/48051