Heliocentric position of object - kepler orbit

orbit.js

  const sin = Math.sin;
  const cos = Math.cos;  

  /**
   * Get heliocentric position of object at a given JD.
   * @param {Number} jd Date value in JD.
   * @return {Array.<Number>} [X, Y, Z] coordinates
   */
  getPositionAtTime(jd) {
    const eph = this._ephem;

    // Pull in ephemeris
    const p = eph.get('w_bar', 'rad'); // Longitude of Perihelion
    const o = eph.get('om', 'rad');    // Longitude of Ascending Node
    const ma = eph.get('ma', 'rad');   // Mean Anomaly
    const n = eph.get('n', 'rad');     // Mean Motion
    const e = eph.get('e');            // Eccentricity
    const a = eph.get('a');            // Semimajor Axis
    const i = eph.get('i');            // Inclination

    // Calculate mean anomaly at JD
    const epoch = eph.get('epoch');
    const M = ma + n * (jd - epoch);

    // Estimate eccentric and true anom using iterative approximation
    let E0 = M;
    let lastdiff;
    do {
      const E1 = M + e * sin(E0);
      lastdiff = Math.abs(E1 - E0);
      E0 = E1;
    } while (lastdiff > 1e-6);
    const E = E0;
    const v = 2 * Math.atan(Math.sqrt((1 + e) / (1 - e)) * Math.tan(E / 2));

    // Radius vector in AU
    const r = a * (1 - e * e) / (1 + e * cos(v));

    // Convert to heliocentric coords
    const X = r * (cos(o) * cos(v + p - o) - sin(o) * sin(v + p - o) * cos(i));
    const Y = r * (sin(o) * cos(v + p - o) + cos(o) * sin(v + p - o) * cos(i));
    const Z = r * (sin(v + p - o) * sin(i));
    return [X, Y, Z];
  }