Compute heliocentric ecliptic coordinates from Kepler's laws

I'm writing a program to compute the heliocentric ecliptic position of planets and asteroids. I want to output rectangular coordinates oriented according to the ecliptic: positive X in the direction of vernal equinox, positive Z toward the north ecliptic pole.

My solution roughly matches other publicly available Kepler simulations (just eyeballing), but I think it is slightly off. What is wrong and how do I verify it?

Here is my process for calculating position at a given JD:

1. Pull in orbital elements from JPL SSD website. Specifically, the EM-Barycenter elements from this table.

2. Compute period. For example, in the case of Earth:

METERS_IN_AU = 149597870700
SECONDS_IN_DAY = 86400

GM_SUN = 1.3271244004193938E+20              // m^3/s^2
SEMIMAJOR_AXIS = 1.0000026 * METERS_IN_AU    // meters

period = 2 * PI * sqrt(SEMIMAJOR_AXIS^3 / GM_SUN)
periodInDays = period / SECONDS_IN_DAY

GM is from NASA NAIF. The above formula returns periodInDays=365.2583228291575, as opposed to the commonly cited 365.2422 days in a year.

1. Compute meanMotion = 2 * PI / periodInDays (radians per day)

2. Compute mean anomaly at JD=2458563.415278 (spring equinox):

epoch = time of orbital elements (jd)
meanLon = mean longitude (radians)
lonPeri = longitude of perihelion (radians)

meanAnomalyAtEpoch = meanLon - lonPeri;
M = meanAnomalyAtEpoch + meanMotion * (2458563.415278 - epoch);

This can produce an unusually large angle. In this case, M=120.6878 radians!

1. Estimate eccentric anomaly E using iterative approximation:
let E = M;
let lastdiff;
do {
const E1 = M + e * sin(E);
lastdiff = abs(E1 - E);
E = E1;
} while (lastdiff > 1e-6);
1. Compute true anomaly from eccentric anomaly E and eccentricity e:
v = 2 * atan(sqrt((1 + e) / (1 - e)) * tan(E / 2))
1. Compute radius in AU:
r = a * (1 - e * e) / (1 + e * cos(v));
1. Compute heliocentric coordinates:
o = longitude of ascending node (radians)
p = longitude of perihelion (radians)
i = inclination (radians)

X = r * (cos(o) * cos(v + p - o) - sin(o) * sin(v + p - o) * cos(i))
Y = r * (sin(o) * cos(v + p - o) + cos(o) * sin(v + p - o) * cos(i))
Z = r * (sin(v + p - o) * sin(i))

In the case of Earth with these elements, I get the following outputs for JD=2443568.0 (code here):

Period (days): 365.25832830801806
Mean motion (rad/day): 0.01720203160400233
Mean anomaly at epoch (rad): -0.043163916976385774
Mean anomaly at JD 2458563.415278 (rad): 120.68783750518942
Eccentric anomaly (rad): 120.70404042143021
True anomaly (rad): 1.3397564871245002