Simulating a restricted 3-body system with On-Rails timewarping

[Whirligig Girl]

KSP uses the "Patched Conics" system for orbital mechanics. What this means is that you patch together "Spheres of Influence" (roughly equivalent to real-world Hill Spheres) to simulate only the gravity of the most prevalent gravity source. All trajectories are just conic sections patched together into ellipses. This system is also called "Restricted Two-Body" because you use the Parent Body (Planet) with gravity and the Satellite (Spacecraft) without gravity. The advantage of this system is that the position and velocity of the spacecraft is a function of time. This means you can run time forwards (or backwards!) freely, and it will always be accurate. Compare to N-Body simulations which can not warp on rails, so timewarp introduces a whole lot of lag/inaccuracies.

The question is, can you simulate the two most gravitationally significant bodies at one time, whilst still being able to maintain an on-rails system where location is a function of time.

Why would anyone want this? Perhaps a space simulation game would want to simulate, say, Duna and Ike together. (Or Pluto and Charon). The two are actually a binary system, so the gravity of Ike would significantly affect spacecraft orbiting Duna, especially during a transfer orbit flight.

And a quick bit about how planets would be simulated under this model: It's the same deal as the restricted 3-body system except all bodies have gravity. Ike and Duna both move together, and around the sun. You would not need a barycenter SOI. (Except perhaps as a boundary SOI, for when Kerbol's gravity is actually more influential than Ike's gravity. Only simulating Duna's gravity without factoring in Ike is like the difference between weighing a bowling ball versus a canonball. They have to have their gravity affect the craft at all times when inside Duna's Hill Sphere. So you probably need a sort of barycenter artificially determined area of influence to transition from Duna+Kerbol to Duna+Ike.

To be clear, I'm using Duna and Ike as an example, not as an idea for a mod.

+Rep and also sheer amazement if one of you advanced coder-type persons can manage to build a prototype of this. Just a star and two binary-orbit planets and a spacecraft which can change it's velocity, and it can transfer to on-rails warp. I know that if it's possible, some of you programmer-type people are skilled enough to do it. Think of it as a challenge. Not a contest or a prize, but a challenge.

[EzinX]

The math for 3 body is much much much harder than the math for 2 body. And it doesn't scale. The correct method to use that will work to N-bodies is described here : http://http.developer.nvidia.com/GPUGems3/gpugems3_ch31.html

By using an octree to represent far away clouds of particles, the actual computational load is much less than N^2. It's not even that complicated - the code listings for all the key parts are right there in the article, and I've heard of a half dozen folks who have gotten such a simulation running in their spare time. Also, using proper N-body gravity means it would not be any big deal to have ships running their engines during timewarp so that we don't have to wait all day for ion drive burns.

[Whirligig Girl]

The math for 3 body is much much much harder than the math for 2 body. And it doesn't scale. The correct method to use that will work to N-bodies is described here : http://http.developer.nvidia.com/GPUGems3/gpugems3_ch31.html

By using an octree to represent far away clouds of particles, the actual computational load is much less than N^2. It's not even that complicated - the code listings for all the key parts are right there in the article, and I've heard of a half dozen folks who have gotten such a simulation running in their spare time. Also, using proper N-body gravity means it would not be any big deal to have ships running their engines during timewarp so that we don't have to wait all day for ion drive burns.

What do you mean "doesn't scale"? And this doesn't look like it could actually accellerate time without errors accumulating. Anyway, finding an alternate solution is not the point, the point is whether restricted 3-body on rails would work.

[EzinX]

What do you mean "doesn't scale"? And this doesn't look like it could actually accellerate time without errors accumulating. Anyway, finding an alternate solution is not the point, the point is whether restricted 3-body on rails would work.

By "doesn't scale", I mean just that. You could maybe have 1 planet and a moon considered at the same time, but what about 2 moons? What about 1 moon and 1 asteroid and 1 planet? And so on and so forth. The correct solution is much easier and simpler to program than what you're asking for.

As for errors accumulating, actually, if you use symplectic integration, the errors for energy are guaranteed to remain bounded. What that means is, the actual position of the simulated bodies (aka phase angle) will vary over time, but the energy in their orbits will not. That means that stuff won't fall out of orbit, etc. The current method KSP uses does not preserve phase angle, either, so this drawback is not a detriment. (I don't know the exact reason it loses phase angle, but I know that if you put a satellite into perfect geosync it will drift anyway)

[Whirligig Girl]

By "doesn't scale", I mean just that. You could maybe have 1 planet and a moon considered at the same time, but what about 2 moons? What about 1 moon and 1 asteroid and 1 planet? And so on and so forth. The correct solution is much easier and simpler to program than what you're asking for.

As for errors accumulating, actually, if you use symplectic integration, the errors for energy are guaranteed to remain bounded. What that means is, the actual position of the simulated bodies (aka phase angle) will vary over time, but the energy in their orbits will not. That means that stuff won't fall out of orbit, etc. The current method KSP uses does not preserve phase angle, either, so this drawback is not a detriment. (I don't know the exact reason it loses phase angle, but I know that if you put a satellite into perfect geosync it will drift anyway)

Well I don't want 1 planet, 1 moon, and 1 asteroid. And I don't want to feel the gravity of two moons. I only want to deal with 2 sources of gravity at a time. This is not about solving the actual problem, it's about making a restricted 3-body on-rails and off-rails simulation. It's a challenge, like building a ternary computer or even a simple robot that unplugs itself.

And what happens in the N-body solution here when you go up to very high timewarps? 1,000,000x timewarp? KSP can handle it. (As long as you don't change SOIs)

I think the key of this is finding the recursion time of a system. The recursion time of a restricted 2-body system is the time it takes to complete one orbit. But for a 3-body system, the recursion time will obviously be much longer, but eventually (As long as an orbit is stable) there will be a recursion.

Edited December 5, 2014 by GregroxMun

[EzinX]

And what happens in the N-body solution here when you go up to very high timewarps? 1,000,000x timewarp? KSP can handle it. (As long as you don't change SOIs)

To handle higher timewarps, you have to start to sacrifice positional accuracy to a larger degree, which you don't want to do with planets. So, at least for Sol, baked solutions exist that work at any time warp like http://en.wikipedia.org/wiki/VSOP_%28planets%29#VSOP87 . To do what you would like, you would need to create an analytical solution for the Kerbin system for all the major planets and moons that was accurate for several thousand years, just like VSOP87. This is essentially a gigantic equation that will tell you where any given body is at any given time, specific to the major planets in Sol.

A solution like this is what you're really asking for : you would probably dig up the methods and code used to created VSOP87, then simulate the bodies of KSP over time using a powerful computer and then generate the equations from that empirical data.

[Whirligig Girl]

No, all I'm asking for is a star and a binary set of planets. Not the Kerbolar system.

If I wanted to make a whole solar system, then for the planets I'd just have them in normal orbits. And at any rate, I'm more concerned with this spacecraft.

Edited December 5, 2014 by GregroxMun

[Mattasmack]

On-rails timewarping as you describe it requires an analytic solution for the three-body problem. In general the gravitational three-body problem has no analytic solution. The only three-body systems that do have analytic solutions have very particular restrictions. For example, in addition to the third body having negligible mass, there are analytical solutions if the two massive bodies move around each other in circular orbits and all motion is restricted to the plane of that orbit. In more general cases, no. The motion is usually chaotic; there is no recursion time.

Also, regarding restricted N-body simulations of the Kerbal solar system, yes, they can run at 1,000,000x time warp and faster. In the context of KSP the game (i.e., the simulation only need run for ~100 years rather than geologic time periods and only need be accurate enough to produce reasonable-looking behavior in the game), it is even easy.

[K^2]

Indeed. In order to get on-the-rails solution, you must first find constants of motion. There are 6 for a general orbit. In 2-body problem, these are trivially related to the 6 orbital elements. Furthermore, Energy and Angular Momentum account for 3 of the 6, which leaves a time and two angle references as your other 3.

A 3-body restricted problem is not conservative. So energy is not a constant of motion. Neither is reference time. When you can't use time as a coordinate, you know you are in trouble. If the two massive bodies are further restricted to a circular orbit, then you can go into a rotating frame of reference, where energy is conserved. But you still have only one component of angular momentum. So that leaves two constants of motion and two coordinates to disentangle. If orbits are not circular, then you are completely screwed.

I'm not entirely sure why planar motion is a requirement. Intuitively, it feels like cylindrical symmetry here should let you generalize to any restricted 3-body where the two massive bodies are on a fixed circular orbit with respect to each other. I'm going to play with Hamilton-Jacobi for this setup in cylindrical coordinates and see if I can come up with anything half-useful.

Edit: Took a look at your link. It looks to me like 2D case is considered for simplicity. Same methodology should be applicable for general 3D motion. It's just going to be way scarier. But this would suggest that at least something like Kerbin-Mun system could be done as restricted 3-body on rails, giving these two their own, combined SOI. (Minmus would still need its own SOI.)

Edited December 5, 2014 by K^2

[Mattasmack]

It may be that the planar case was considered there for simplicity, but I've seen it mentioned in multiple places as a restriction so I felt safe in claiming that. But this is getting outside of my experience, so I would not be terribly surprised to be shown to be wrong.

[K^2]

I've just realized, as I was writing it all down, that centrifugal potential breaks the cylindrical symmetry in the rotating frame. So it probably needs to be planar, after all. At least, I see no obvious way of separating this crocodile of an equation.

[EzinX]

I've just realized, as I was writing it all down, that centrifugal potential breaks the cylindrical symmetry in the rotating frame. So it probably needs to be planar, after all. At least, I see no obvious way of separating this crocodile of an equation.

So why not just go ahead and do it with numerical integration for the N-body, using an octree to group together and average farther away bodies? GPU-gems makes it look so easy.

[Mattasmack]

I've just realized, as I was writing it all down, that centrifugal potential breaks the cylindrical symmetry in the rotating frame. So it probably needs to be planar, after all. At least, I see no obvious way of separating this crocodile of an equation.

Ah, yeah, that would do it.

So why not just go ahead and do it with numerical integration for the N-body, using an octree to group together and average farther away bodies? GPU-gems makes it look so easy.

Because the context of this thread is GregoxMun asking about analytic solutions to a restricted three-body problem.

Edited December 6, 2014 by Mattasmack

[YNM]

Restricted three body problems are hard to a fair extent (Least are, harder than restricted two-body problems or harder than two body problems). As far as I can tell the only easy solution is for some special points where you can predict it'll stay there relative to the two bodies (ie. Lagrangian points). Other than that, you need to calculate the forces everytime, probably nearing infinite amount of calculations ; but as we can't reach that, you need to simplify it down, maybe calculating it every tens of milisecond. I know you're asking for a trajectory, but when you start to talk about N-body problem (restricted three-body problem is already one of it - the lowest, easiest of all) you should realize most movements will be chaotic, even in infinite time. No real predictions exist unless you already run the simulation to all those time - even it's not an unusual thing to see the true end is different compared to what's simulated. If you ever wonder why probes (or other things) mostly do mid-course correction, it's one of the reason.

Edited December 6, 2014 by YNM
adding some clarity

[K^2]

Also, while not a huge problem with 3 bodies, a general N-body 1/R potential problem is a numerical nightmare. Don't trust any game-oriented literature on the subject. You can run a plausible-looking galaxy with millions of stars using GPU, but you wouldn't be able to simulate trajectory of a comet well enough to say for sure if it will hit your planet within a few years with just a handful of significant sources of gravity in the system. And for a KSP-like game, ability to predict precisely where you'll end up is the key.

[SciMan]

Everybody keeps coming up with "computational load" as a negative factor.

Let me tell you, it's a complete and total NON-issue.

For evidence, I present Orbiter. Yes, the space sim that (partially?) inspired KSP.

It uses the VSOP87 solution that was mentioned several posts above mine, coupled with the ELP2000 solution for the Earth-Moon system.

That covers all 8 "official" planets, Earth's Moon, the two moons of Mars, 4 moons of Jupiter, and a bunch around Saturn and Neptune as well.

It does all that, and it has time warp that allows engine thrust, atmospheric flight, and control inputs. And the time warp even tracks rotation.

And it does that all at up to 600 frames per second. I could easily keep the simulation pegged at 600fps on a 7-8 year old computer, even at full time warp.

And it can simulate non-spherical gravity, gravity gradient torque, Lagrange points and a bunch of other really tiny details that most people don't even know exist.

All of that means that the simulation is accurate and stable enough that the limiting factor is the accuracy of the burns that you perform, not the accuracy of the simulation.

A fully autopilot flown mission in Orbiter will be able to choose which square meter of it's target planet it lands on, and tell you how much fuel you need down to a single cubic centimeter or better. And then it can do it again, as many times as you want or have fuel for.

Bottom line? If this idea is ever implemented, the load of the N-body physics calculations won't be what kills it. Assuming anything "kills" it, of course. I doubt anything except a Unity bug could do such a thing.

Of course, anyone attempting to make this system in KSP will probably need to re-work twenty other things just to get them all to agree on where the ship is going (maneuver nodes, I'm looking at you).

[EzinX]

Also, while not a huge problem with 3 bodies, a general N-body 1/R potential problem is a numerical nightmare. Don't trust any game-oriented literature on the subject. You can run a plausible-looking galaxy with millions of stars using GPU, but you wouldn't be able to simulate trajectory of a comet well enough to say for sure if it will hit your planet within a few years with just a handful of significant sources of gravity in the system. And for a KSP-like game, ability to predict precisely where you'll end up is the key.

So, for a KSP like game, you want to be able to specify a rendezvous with, say, that comet, and perform a series of burns to set up the intercept. With true N-body gravity, everything affects everything to a degree, and you cannot serially run the simulation beyond the serial processing speed of your GPU.

Well, we can't make that intercept perfectly in real life, either. Corrective burns have to be done. So, what you *could* do is the following :

Plan the trajectory to reach the intercept position (a closest approach to the target) by using long timesteps and high order symplectic integration over a brief period of realtime. Essentially, you now have a rough trajectory that using a lower precision integration method than your game's real simulation uses (the longer timesteps cut the precision), ends up somewhere roughly where you want to be.

If the player approves the flight plan, you would perform corrective burns automatically during the flight if the relative differencee between the planned trajectory and the actual trajectory were tiny enough to correct this way.

Another approach is one even you probably haven't heard of, K^2. There is a method to simulate a serial trajectory in parallel, despite trajectories being a series of mathematically dependent states. PM me if you'd like an explanation of how to use the method. Using that method, and given access to a monster GPU, you could run at at full simulation precision at millions of times speedup and predict a trajectory perfectly.

[K^2]

long timesteps and high order symplectic integration over a brief period of realtime.

And there is no symplectic method for 1/R potential. So you are down to implicit integration over collocation points, or some variation on that, which tends to be extremely time consuming. It also can't be broken into parallel computations, because you are solving a non-linear optimization problem at each point. Or we'd all be doing that.

[EzinX]

What's this 1/R potential? At any given point along a spacecraft's trajectory, it's subject to forces from gravity and from it's own engine. Those forces always sum to an acceleration vector. That vector varies over time. You have a starting position, and you want to know the new velocity and new position. Boom, you just invented Euler's method.

You then use a fancier method to compensate for some of the changes of acceleration over the timestep and you end up sampling at several past states and you end up with 4th order symplectic. Mathematically proving it or studying effects of the integrator and it's error over time is a whole field of study, but, regardless, you have something that works pretty well.

I take it this 1/R problem is some change in energy from close approaches to planets? How big of a change are we talking about here?

With a game, you can cheat in several ways that scientists running simulations cannot. For one thing, you can simply declare that the planets obey a pre-baked solution like VSOP87. Even the longest KSP campaign is realistically not more than a century anyway, so some of the numerical integration methods are probably good enough that you won't have them crashing into each other or drifting too far from their orbits. You only really care about N-body effects on a spacecraft and so your numerical simulation does not have to be infinitely accurate.

[K^2]

Do you actually know what "symplectic integrator" means? I have no idea how you've managed to pick up that term and not understand something about potentials.

In the context of gravitational interaction, and without getting into details of symplectic geometry, symplectic integrator would be the one that keeps energy within roughly the correct value. It can change due to numerical errors, but these errors don't explode over time.

For example, if you are simulating a spring, Forward Euler is not symplectic. But then pretty much any second order method is. Which is why nearly every half-decent game simulation is done using Velocity Verlets. It's enough to prevent your collisions from "exploding" the system.

Gravity has no symplectic method. Any integration method will end up with energy being lost or generated in the system. And that very quickly causes the system to diverge from true solution in chaotic and unpredictable way. There is no fix for it. You can keep padding your numerical accuracy and reducing the step size, using implicit methods for every step, but you will eventually have to give up. The reason we can do better with planets is because we use classical perturbation theory to improve our computations. These are analytical methods that let you take into account the most significant deviations from Kepler. So when we get down to numerical integration, we aren't computing the trajectory, but just small deviations from it.

If you have skill to actually write a numerical integrator, I suggest you try this. JPL's Horizons will let you get NASA's solutions for trajectories of major bodies in Solar System. You can get it to dump literal XYZ coordinates with respect to your chosen origin. (I recommend using Solar System's barycenter, because you'll be in accelerated frame otherwise.) Grab the data for the Sun, all planets, major moons, and anything else you might consider a significant source of gravity. Then get trajectory of some comet or asteroid. (67P/Churyumov–Gerasimenko, for example?) Try to take its initial position and velocity, then use your favorite integration method, and see how close you get to the actual trajectory. Don't forget to compare kinetic energy at different points.

After you have tried this, I think we'll end up having a completely different conversation about N-body dynamics.

Edited December 7, 2014 by K^2