Using 2D position, velocity, and mass to determine the parametric position equations for an orbiting body

Question

I have a gravity-related question. I am programming an orbit simulator. I have everything up and running, but I would like to render the smaller body's orbital path (the larger body is fixed). To do this, I need the parametric equations for the body's position for any given time (i.e. I need an $x(t)$ and $y(t)$ function). I have access to each body's position, velocity, and mass. I have no other accessible variables.

P.S. Here is a link to what I have so far. The equations currently used to plot the orbit are $x(t)=x+v_xt$ and $y(t)=y+v_yt$.

Answer 1

This answer is an complement to Chris White's answer.

Fist of there is no explicit equations for the position of an object following a Kepler orbit as a function of time. However, when the initial conditions are known, the path the object will follow can be found, as well as the velocity, acceleration, ect. at any given position. This path can be described by the following equation: $$ r=\frac{a(1-e^2)}{1+e\cos{\theta}} $$ where $r$ is the distance between the (small) object and the focus of the orbit (which is equal to the position of the large fixed object or the center of mass of the two if the second object it not fixed), $a$ is the semi-major axis, $e$ is the eccentricity and $\theta$ is the true anomaly (the angle between the objects position, focus and the position of the object at closest approach/periapsis).

The semi-major axis and the eccentricity can be deduced from the initial conditions: $$ a=\frac{\mu r}{2\mu-rv^2} $$ $$ e=\sqrt{1+\frac{rv_{\theta}^2}{\mu}\left(\frac{rv^2}{\mu}-2\right)} $$ where $\mu$ is the gravitational parameter of the large fixed object ($\mu=GM$), $v$ is the velocity of the object and $v_\theta$ is the tangential component of the velocity.

This tangential component of the velocity can be found by using some linear algebra (the resulting equivalent of the cross product in 2D): $$ v_{\theta}=\frac{xv_y-yv_x}{\sqrt{x^2+y^2}} $$

The position of the periapsis may be have any angle relative to the focus and the axis in which you represent positions. So you will have to offset/rotate the path you have found by a specific angle to match the initial conditions.

This offset angle of $\theta$, lets call it $\Delta\theta$, can be calculated as followed: $$ \Delta\theta={\rm sign}\left(v_{\theta}v_r\right)\cos^{-1}\left(\frac{a(1-e^2)-\sqrt{x^2+y^2}}{e\sqrt{x^2+y^2}}\right)-{\rm atan2}\left(y,x\right) $$ where the function ${\rm sign}(x)$ refers to whether $x$ is positive or negative $\left(\frac{x}{|x|}\right)$ and $v_r$ is the radial velocity, which can be found by: $$ v_r=\frac{xv_x+yv_y}{\sqrt{x^2+y^2}} $$ This offset angle is chosen such that the resulting path is defined by: $$ x_{path}=r(\theta)\cos\left(\theta+\Delta\theta\right) $$ $$ y_{path}=-r(\theta)\sin\left(\theta+\Delta\theta\right) $$

I have one final note about choosing the range of $\theta$, since when $e\ge1$ then the path will no longer be an ellipse, but a parabola or a hyperbola, which extends to infinity. To avoid trying to draw this you can limit the range of $\theta$. For example by choosing a maximum radius, $r_{max}$. The range of $\theta$ will then be: $$ \theta_{max}={\rm real}\left(\cos^{-1}\left(\frac{a(1-e^2)-r_{max}}{er_{max}}\right)\right) $$ $$ \theta\in[-\theta_{max},\theta_{max}] $$

Answer 2

It seems you've done the hard part already, which is to evolve the object's position as a function of time. And moreover, the simulation seems stable over a number of orbits. (But eventually things start to go wrong; you may want to look at an answer I wrote to What is the correct way of integrating in astronomy simulations?)

So my understanding is all you really need is the full orbit plotted. In that case, there's no need to make it a function of actual time (that's hard). Instead, we can fall back to Kepler's First Law, which says that the separation $r$ between the bodies obeys $$ r = \frac{r_\mathrm{max}(1-e)}{1+e\cos(\theta)}. $$ $r_\mathrm{max}$ is the maximum separation, which I believe is the initial separation in your particular simulation. The eccentricity $e$ (formula here) is given by $$ e^2 = 1 + \frac{2Er^4\dot{\theta}^2}{G^2M^2m}, $$ where $E$ is the orbital energy, $G$ is the gravitational constant, $M$ is the mass that is so heavy it essentially doesn't move, and $m$ is the lighter mass. This can more conveniently be rewritten $$ e^2 = 1 + 2 \left(\frac{rv}{GM}\right)^2 \left(\frac{v^2}{2} - \frac{GM}{r}\right), $$ where $v$ is the velocity of the moving particle.

Note that the convention is for $\theta$ to measure the angle from closest approach (i.e., $\theta = 0$ is the negative $y$-axis in your simulation). If you want $\theta$ to increase with time, then the (nonstandard) transformation to Cartesian coordinates is \begin{align} x & = -r \sin(\theta) \\ y & = -r \cos(\theta), \end{align} assuming $r = 0$ corresponds to the massive object. In any event, to plot the orbit all you need to do is calculate the constant $e$ at one point in the orbit, sample $\theta$ with as many points as you like, calculate the separations $r(\theta)$, and convert the $(r, \theta)$ pairs to $(x, y)$.

Answer 3

Orbital simulations can be handled by using the following relations: \begin{eqnarray} \mathbf F&=&m\mathbf a=m\frac{d^2\mathbf x}{dt^2}\tag{a} \\ \mathbf v&=&\frac{d\mathbf x}{dt}\tag{b}\\ \mathbf a&=&\frac{d\mathbf v}{dt}\tag{c} \end{eqnarray} The force acting on any two bodies, mass $M$ and $m$ is given by Newton's gravitational law $$ \mathbf F=G\frac{mM}{r^2}\hat{\mathbf r}=G\frac{mM}{r^3}\left(x\hat{\mathbf x}+y\hat{\mathbf y}\right)\tag{d} $$ with $r=\sqrt{x^2+y^2}$. Then the algorithm is thus:

1. Compute $\mathbf F$ from $\mathbf x$
2. Determine the new velocity $\mathbf v$ via $\mathbf v^{n+1}=\mathbf v^n+\mathbf a\,dt$
3. Determine the new position $\mathbf x$ via $\mathbf x^{n+1}=\mathbf x^n+\mathbf v\,dt$
4. Go back to Step 1.

This can be done using the exact equations given above, but will not be numerically stable for very long. Higher order methods like Runge-Kutta 3/4/5 will keep your stability for longer periods of time, but will not be perfect.

If this is not the information needed, let me know and I'll try helping you out with what I can, but it seems to me that you just needed to know Equation (d) in relation to the others.