Computationally solving the two-body problem for patched conics approximation

Question

I have been trying to create an orbital simulation using the patched conics approximation like Kerbal Space Program. I initially tried solving Kepler's equation using Newton's method, however this left much to be desired as it requires many initial variables such as the semimajor axis, and eccentricity to be provided up front, where as I only have a radius vector and velocity vector. Also this approach seems to fall apart with orbits that are parabolic or hyperbolic (or rectilinear), not to mention that it is quite computationally intensive. This lead me to the Goodyear method from lamont's answer on this question.

Determining orbital position at a future point in time

However, I have been unable to get it to produce valid results. All of the radius vectors I get from it form straight lines. From what I understand, all that is required is the initial radius vector, velocity vector, gravitational parameter, and time interval. I have been giving these values in SI units, so meters, and seconds. Does it expect different units?

If all of the above fails, are there any other general solutions to the two body problem that only require the initial radius vector and velocity vector, and can return a radius vector and velocity vector at some time 't' in the future? And which of these lends themselves best to a patched conic approximation?

Answer 1

You probably already found an answer by the time you posted this question, but I still give an answer in case it helps.

Using the orbital elements is the "right way" to do it. There are 6 of them, described in this Wikipedia article :

• eccentricity $e$
• semi-major axis $a$
• inclination $i$
• longitude of the ascending node $\Omega$
• argument of periapsis $\omega$
• true anomaly $\theta$ (which is time dependant)

The first five of them completely describe the geometry of an orbit in 3D. The true anomaly is an angle refering to the position of the body/spacecraft on that orbit. You can calculate all of these parameters from initial state vectors (position and velocity), knowing the standard gravitational parameter of the attractor body. The mathematics are described in this PDF by René Schwarz on how to convert cartesian orbital state vectors to orbital elements.

Once you have the elements of an orbit, you can use them to compute the state vectors (position/radius vector and velocity vector) at any moment in time following this other PDF on how to convert orbital elements to cartesian state vectors. However, some formulas described in the document only work for elliptic orbits. For hyperbolic orbits (when $e > 1$) you'll have to use other formulas for the eccentric and true anomaly (see this Wikipedia article and the second answer on this post), as well as adapt the formula for the velocity vector from the document (see this post).

The "complicated" looking rotation matrices that set the position and velocity of the spacecraft/body in the 3D space can be developed in these steps:

• Choose the reference vectors for "right" and "up", e.g. the x direction (1, 0, 0) and y direction (0, 1, 0), I call them $\vec{u_x}$ and $\vec{u_y}$.
• Compute the ascending node vector $\vec{u_{asc}}$ : it is $\vec{u_x}$ rotated around $\vec{u_y}$ by an angle $\Omega$;
• Set the position $\vec{r} = \vec{o}$ and the velocity $\vec{v} = \dot{\vec{o}}$ as described in the document (with the modified formula in case $e > 1$);
• Rotate $\vec{r}$ and $\vec{v}$ around $\vec{u_y}$ of an angle $\Omega + \omega$;
• Rotate again $\vec{r}$ and $\vec{v}$ around $\vec{u_{asc}}$ of an angle $i$.

(Note that in René Schwarz's documents, he considers the z-axis as the "up" axis, instead of the y-axis)

The case of a parabolic orbit ($e = 1$) can be, as notovny said in a comment, politely ignored (apply the ostrich algorithm).

Note however that you may need to take account for extreme cases. i.e. when you are dealing with orbits that are really close to circular orbits or have a null inclination. In these cases the formulas for angles, such as the true anomaly in the first document, won't work (they will give NaN values because of floating point precision and roundings). Thus for these particular cases you'll have to set default vectors, e.g. the eccentric vector, and calculate the angles from it.

As for units, the formulas indeed expect SI units. (You can however use any multiple of these units as long as you account for it in the equations).