Calculate xyz velocity vectors for circular orbit

Question

I'm trying to calculate the xyz velocity vectors for a circular orbit given a set of position vectors, but I can't seem to get it right.

The formula for obtaining the velocity magnitude is the square root of gravitational constant multiplied by the mass of the primary divided by the distance from the primary to the object orbiting it. I'm using solar masses for mass, astronomical units for distance and years for time. To simplify things, the primary is our sun, which means that G * M is equal to 39.5. The position of the Sun is 0, 0, 0 and the vectors you see for Earth are in au and were obtained from the JPL Horizons tool.

Here's the code:

    function getDistanceParams(p1, p2) {
      const dx = p2.x - p1.x;
      const dy = p2.y - p1.y;
      const dz = p2.z - p1.z;

      return { dx, dy, dz, dSquared: dx * dx + dy * dy + dz * dz };
    }

    function getCircularOrbit() {
      const dParams = getDistanceParams(
        { x: 0, y: 0, z: 0 },
        {
          x: 0.9197324105567349,
          y: -0.4147318273536994,
          z: -1.750037390352759e-5
        }
      );

      const vMag = Math.sqrt(39.5 / Math.sqrt(dParams.dSquared));

      return {vx: dParams.dx * vMag, vy: dParams.dy * vMag, vz: dParams.dz * vMag};
    }

    getCircularOrbit();

/*

The above outputs this:

{ vx: 5.754832301294686,
  vy: -2.59500707927136,
  vz: -0.00010950110691846996 }

The actual velocity vectors:

  vx: 2.4645428337894026,
  vy: 5.7097644117945805,
  vz: -3.3177815766459033e-4,

*/

Where am I going wrong???

Update: Code that works

In case anyone is curious, here's the code that ended up doing the job:

export function getDistanceParams(p1, p2) {
  const dx = p2.x - p1.x;
  const dy = p2.y - p1.y;
  const dz = p2.z - p1.z;

  return { dx, dy, dz, dSquared: dx * dx + dy * dy + dz * dz };
}

export function getIdealCircularOrbit(primary, secondary, g) {
  const dParams = getDistanceParams(primary, secondary);

  const d = Math.sqrt(dParams.dSquared);

  const vMag = Math.sqrt(g * primary.m / d);

  return {
    x: secondary.x,
    y: secondary.y,
    z: secondary.z,
    vx: primary.vx + -dParams.dy * vMag / d,
    vy: primary.vy + dParams.dx * vMag / d,
    vz: primary.vz + dParams.dz * vMag / d
  };
}

Answer 1

In your code, you seem to be calculating the magnitude $|\vec{v}|$ of the velocity (the speed) from your gravitational formula and multiplying it by the components of the position vector $\vec{r}$. There are two errors here.

1. You need to multiply the speed by a unit vector in the direction of the velocity, in order to get the velocity vector $\vec{v}$. In general the position vector is not a unit vector (actually in this case, it is pretty close, because you chose astronomical units, but that is just luck).
2. The position vector $\vec{r}$ is pointing in a radial direction, whereas the velocity vector $\vec{v}$ (for circular motion) lies in a tangential direction.

You are not actually inputting to your program enough information to work out the direction of the velocity vector. You need to know the plane in which the rotation is taking place, and the sense of the rotation (clockwise or anticlockwise). The key equation is $$ \vec{v}=\vec{\omega}\times\vec{r}=\omega\hat{n}\times\vec{r} $$ where $\vec{\omega}=\omega\hat{n}$ is the angular velocity vector and $\times$ is the vector cross product. Assuming that you know the unit vector $\hat{n}$ pointing in the direction of $\vec{\omega}$, you can deduce the direction of $\vec{v}$ from the above formula. Alternatively, you can calculate the magnitude $\omega$ from the speed and the radius of the orbit, and use it together with $\hat{n}$ directly in the above formula.

Of course, if you are prepared to make the approximation that the rotation is entirely in the $xy$ plane (which is not exactly correct, from your numbers, but quite close) then you know $\hat{n}$ and can easily work out the direction of $\vec{v}$ from $\vec{r}$.