Initial velocity vector for circular orbit

Question

I'm trying to create a solar system simulation, and I'm having problems trying to figure out initial velocity vectors for random objects I've placed into the simulation.

Assume: - I'm using Gaussian grav constant, so all my units are AU/Solar Masses/Day - Using x,y,z for coordinates - One star, which is fixed at 0,0,0. Quasi-random mass is determined for it - I place a planet, at a random x,y,z coordinate, and its own quasi-random mass determined.

Before I start the nbody loop (using RK4), I would like the initial velocity of the planet to be such that it has a circular orbit around the star. Other placed planets will, of course, pull on it once the simulation starts, but I want to give it the chance to have a stable orbit...

So, in the end, I need to have an initial velocity vector (x,y,z) for the planet that means it would have a circular orbit around the star after 1 timestep.

Help? I've been beating my head against this for weeks and I don't believe I have any reasonable solution yet...

Answer 1

It is quite simple if you assume that the mass of the star M is much bigger than the total mass of all planets sum(m[i]). This simplifies the problem as it allows you to pin the star to the centre of the coordinate system. Also it is much easier to assume that the motion of all planets is coplanar, which further reduces the dimensionality of the problem to 2D.

1. First determine the magnitude of the circular orbit velocity given the magnitude of the radius vector r[i] (the radius of the orbit). It only depends on the mass of the star, because of the above mentioned assumption: v[i] = sqrt(mu / r[i]), where mu is the standard gravitational parameter of the star, mu = G * M.

2. Pick a random orbital phase parameter phi[i] by sampling uniformly from [0, 2*pi). Then the initial position of the planet in Cartesian coordinates is:
   x[i] = r[i] * cos(phi[i])
y[i] = r[i] * sin(phi[i])

3. With circular orbits the velocity vector is always perpendicular to the radial vector, i.e. its direction is phi[i] +/- pi/2 (+pi/2 for counter-clockwise (CCW) rotation and -pi/2 for clockwise rotation). Let's take CCW rotation as an example. The Cartesian coordinates of the planet's velocity are:
   vx[i] = v[i] * cos(phi[i] + pi/2) = -v[i] * sin(phi[i])
vy[i] = v[i] * sin(phi[i] + pi/2) = v[i] * cos(phi[i])

This easily extends to coplanar 3D motion by adding z[i] = 0 and vz[i] = 0, but it makes no sense, since there are no forces in the Z direction and hence z[i] and vz[i] would forever stay equal to 0 (i.e. you will be solving for a 2D subspace problem of the full 3D space).

With full 3D simulation where each planet moves in a randomly inclined initial orbit, one can work that way:

1. This step is equal to step 1 from the 2D case.

2. You need to pick an initial position on the surface of the unit sphere. See here for examples on how to do that in a uniformly random fashion. Then scale the unit sphere coordinates by the magnitude of r[i].

3. In the 3D case, instead of two possible perpendicular vectors, there is a whole tangential plane where the planet velocity lies. The tangential plane has its normal vector collinear to the radius vector and dot(r[i], v[i]) = 0 = x[i]*vx[i] + y[i]*vy[i] + z[i]*vz[i]. One could pick any vector that is perpendicular to r[i], for example e1[i] = (-y[i], x[i], 0). This results in a null vector at the poles, so there one could pick e1[i] = (0, -z[i], y[i]) instead. Then another perpendicular vector can be found by taking the cross product of r[i] and e1[i]:
   e2[i] = r[i] x e1[i] = (r[2]*e1[3]-r[3]*e1[2], r[3]*e1[1]-r[1]*e1[3], r[1]*e1[2]-r[2]*e1[1]). Now e1[i] and e2[i] can be normalised by dividing them by their norms:
   n1[i] = e1[i] / ||e1[i]||
n2[i] = e2[i] / ||e2[i]||
   where ||a|| = sqrt(dot(a, a)) = sqrt(a.x^2 + a.y^2 + a.z^2). Now that you have an orthogonal vector basis in the tangential plane, you can pick one random angle omega in [0, 2*pi) and compute the velocity vector as v[i] = cos(omega) * n1[i] + sin(omega) * n2[i], or as Cartesian components:
   vx[i] = cos(omega) * n1[i].x + sin(omega) * n2[i].x
vy[i] = cos(omega) * n1[i].y + sin(omega) * n2[i].y
vz[i] = cos(omega) * n1[i].z + sin(omega) * n2[i].z.

Note that by construction the basis in step 3 depends on the radius vector, but this does not matter since a random direction (omega) is added.

As to the choice of units, in simulation science we always tend to keep things in natural units, i.e. units where all computed quantities are dimensionless and kept in [0, 1] or at least within 1-2 orders of magnitude and so the full resolution of the limited floating-point representation could be used. If you take the star mass to be in units of Solar mass, distances to be in AUs and time to be in years, then for an Earth-like planet at 1 AU around a Sun-like star, the magnitude of the orbital velocity would be 2*pi (AU/yr) and the magnitude of the radius vector would be 1 (AU).

Answer 2

As Ed suggested, I would use the mks units, rather than some other set of units.

For the initial velocity, I would agree with part of what Ed said, but I would use the vector form of the centripetal acceleration:

m₁v²/r r(hat) = G m₁ m₂ / r2 r(hat)

Set z to 0, and convert from polar coordinates to cartesian coordinates (x,y). Then, you can assign either y or x an initial velocity, and compute what the other variable is to satisfy the circular orbit criteria. This should give you an initial (Vx,Vy) that you can start your nbody problem from. There should also be quite a bit of literature out there on numerical recipes for nbody central force problems.

Answer 3

Just let centripetal acceleration equal gravitational acceleration.

m₁v² / r = G m₁m₂ / r²

v = sqrt( G m₂ / r )

Of course the star mass m₂ must be much greater than the planet mass m₁ or you don't really have a one-body problem.

Units are a pain in the butt when setting up physics problems. I've spent days resolving errors in seconds vs timestep units. Your choice of AU/Solar Masses/Day is utterly insane. Fix that before anything else.

And, keep in mind that computers have inherently limited precision. An nbody simulation accumulates integration error, so after a million or a billion steps you will certainly not have a circle, regardless of the step duration. I don't know much about that math, but I think stable n-body systems keep themselves stable by resonances which absorb minor variations, whether introduced by nearby stars or by the FPU. So the setup might work fine for a stable, 5-body problem but still fail for a 1-body problem.