# Computationally solving the two-body problem for patched conics approximation

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?

• Just a guess: you mention gravitational parameter, and then kilograms. Gravitational parameter ($GM$) doesn't contain mass in its dimension, its dimension is $L^3 T^{-2}$. If the code expected the gravitational parameter in m3/s2, and you gave the mass in kg, that would be off by many orders of magnitude. Dec 5, 2021 at 13:20
• Thats my mistake. I edited the question to remove kilograms. I never used kilograms in the calculation. I used the gravitational parameter of Earth in m^3/s^2, which was roughly 3.986E14 m^3/s^2. The correct units must've skipped my mind when writing the question. Dec 6, 2021 at 17:21
• if you have gravitational parameter, radial distance vector, and velocity vector you have semimajor axis by way of the specific orbital energy and orbital eccentricity by way of the eccentricity vector. The hyperbolic versions of Kepler's equations flip signs and substitute hyperbolic trig functions, and the parabolic can be treated with its own special set of equations, or politely ignored. Dec 6, 2021 at 21:39
• I'd recommend Shepperd's method from that link in my answer on that other question. You can also implement Danby's method of solving Kepler's equation from David Eagle's Matlab implementation: mathworks.com/matlabcentral/fileexchange/… This answer should help you with testing these algorithms: space.stackexchange.com/questions/20669/… -- you can generate a bunch of Kepler solutions and use that to test all three kinds of solvers. Dec 6, 2022 at 19:35
• David Eagle's repo also has eci2orb1 and orb2eci functions that convert to/from state vectors and keplerian elements, e.g. in this repo: mathworks.com/matlabcentral/fileexchange/… and he's got a glambert.m file there which is a Gooding Lambert solver as well. His eci2orb1, orb2eci, glambert, kepler1 (Danby) and twobody2 (Shepperd) form a good Keplerian solver toolbox. IMO, implement all of them. Dec 6, 2022 at 19:49

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