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

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$.

• This question appears to be off-topic because it is about programming simulations and not a physics concept. – Brandon Enright Feb 13 '14 at 22:57
• It's a bit of a grey area, but the question is intrinsically physics related. It's not about programming, it's about finding parametric equations for a body's orbit using known variables. – Conner Ruhl Feb 13 '14 at 23:07
• Hi Conner, you could simulate having a central body fixed but your simulation won't be accurate or the equations will be quite complicated. See physics.stackexchange.com/questions/98011/… – Brandon Enright Feb 13 '14 at 23:08
• There isn't any solid information present in that thread. Are you saying there is no method of taking the information I have available and developing a set of position equations (integration perhaps)? – Conner Ruhl Feb 13 '14 at 23:12
• @BrandonEnright I agree with Conner here, that this is physics in a programming context, rather than an actual programming question, and that it's on topic here. – David Z Feb 14 '14 at 0:10

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)$.

• Your third equation is only correct when at periapsis, apoapsis or in a circular orbit. Namely the first velocity therm should be the radial velocity. And you also have to keep in mind that the ellipse might be rotated. Will post an answer to address these things as well. – fibonatic Feb 14 '14 at 11:15
• I meant tangential instead of radial. – fibonatic Feb 14 '14 at 12:00

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_r$ 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}]$$

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.

• The question is not about the simulation but about the path along which the object will travel. – fibonatic Feb 14 '14 at 13:57