[WIP][1.8.1, 1.9.1, 1.10.1, 1.11.0–2, 1.12.2–5] Principia—version ‎Καραθεοδωρή, released 2024-04-08—n-Body and Extended Body Gravitation

[eggrobin]

To be sure I also checked whether my approach gives the same results for position and velocity components as in document you linked to, by expanding all the vector products. My final results match

Sounds good then, if you make sure you use the correct anomaly.

For my definition of μ I assume that the semi-major axis is given relative to the barycenter of that body and its parent body. Because for orbital elements the fictitious parent body lies at the barycenter, thus I changed μ to such a value that the gravitational force of this fictitious body would be equal to that of the actual parent body. The resulting μ then becomes μ = G * M^3 / (M + m)^2. I do remember seeing μ = G(M +m), however I good not think of a logical argument why I should use it.

Mostly "because the other one is wrong", see e.g. here for circular orbits of two bodies around their barycentre.

You say "The resulting μ then becomes", so you must have an actual definiton for your μ, but you do not give it, so I have no idea where that expression comes from.

EDIT:

If you want to see how that μ ends up there, here, p. 41-49 is a thorough derivation of the whole 2-body problem, ending with an explicit time-dependent expression (which is mostly useless, since it involves series of Bessel functions, but it's nice to see that the time-dependent solution is ugly; it tells you that you should not expect something which is a direct simple computation to be correct).

Note that μ means "the numerator of the 2-body Newtonian potential" here, namely G(m+M), and this is the same μ that appears, e.g., as an inverse square root in the period (denoted Ä, on page 48).

Edited October 3, 2014 by eggrobin

[DerpyHoovesX]

Do you have alpha release(pack) of this mod, to be able to play?

[Guest]

It's not fully functional yet. You'll need to wait a while longer.

[fibonatic]

Mostly "because the other one is wrong", see e.g. here for circular orbits of two bodies around their barycentre.

You say "The resulting μ then becomes", so you must have an actual definiton for your μ, but you do not give it, so I have no idea where that expression comes from.

EDIT:

If you want to see how that μ ends up there, here, p. 41-49 is a thorough derivation of the whole 2-body problem, ending with an explicit time-dependent expression (which is mostly useless, since it involves series of Bessel functions, but it's nice to see that the time-dependent solution is ugly; it tells you that you should not expect something which is a direct simple computation to be correct).

Note that μ means "the numerator of the 2-body Newtonian potential" here, namely G(m+M), and this is the same μ that appears, e.g., as an inverse square root in the period (denoted Ä, on page 48).

I am still not convinced, since none of your sources explain why G(m+M) is used.

This is how I derived my expression for μ:

Suppose you would like to construct elliptical orbits, with eccentricity e, of two body with mass m1 and m2, with a given semi-major axis a1 relative to the barycenter (a2 can be derived from m1 and m2), which initial conditions would you have to give them, assuming that the barycenter would be the origin of the coordinates system?

For simplicity I will limit this problem to 2D (so inclination can be ignored), use an initial true anomaly of zero, and choose the longitude of the ascending node and argument of periapsis in such a way that unit vector of the initial position will be equal to the x-axis.

The initial positions will be simple, namely:

x1 = a1 * (1 - e)
y1 = 0
a2 = a1 * m1 / m2
x2 = -a2 * (1 - e)
y2 = 0

The initial velocities are harder, for this I will initially use μ1 and μ2 as temporary variables, which will have to be solved later.

By using the following equation it can be shown that:

vx1 = 0
vy1 = sqrt(μ1 * (1 + e) / (a1 * (1 - e)))
vx2 = 0
vy2 = -sqrt(μ2 * (1 + e) / (a2 * (1 - e)))

The total momentum should be zero (barycenter is stationary), so

m1 * sqrt(μ1 * (1 + e) / (a1 * (1 - e))) = m2 * sqrt(μ2 * (1 + e) / (a2 * (1 - e)))
m1^2 * μ1 / a1 = m2^2 * μ2 / a2
m1^3 * μ1 = m2^3 * μ2

From these equations it can already become clear that μ1 = μ2 = G * (m1 + m2) is incorrect. To derive the equation, which I mentioned in an previous post, I will use acceleration and the fact that r1 / r2 will be constant.

r1 / r2 = m2 / m1
F = G * m1 * m2 / (r1 + r2)^2
acc1 = G * m2 / (r1 + r2)^2 = G * m2 / (r1 + r1 * m1 / m2)^2 = G * m2 / (r1 * (m1 + m2) / m2)^2 = G * m2^3 / (r1^2 * (m1 + m2)^2)
acc2 = G * m1 / (r1 + r2)^2 = G * m1 / (r2 * m2 / m1 + r2)^2 = G * m1 / (r2 * (m1 + m2) / m1)^2 = G * m1^3 / (r2^2 * (m1 + m2)^2)
μ1 = G * m2^3 / (m1 + m2)^2
μ2 = G * m1^3 / (m1 + m2)^2

When using these equations for μ1 and μ2 it will be clear that the total momentum will be zero.

EDIT:

I also found out that my integration method, also called position Verlet, is actually not symplectic because it does not use velocities (I am still a noob when it comes to symplectic integrators).

Edited October 3, 2014 by fibonatic

[eggrobin]

I am still not convinced, since none of your sources explain why G(m+M) is used

EDIT:

These lecture notes from ETHZ indeed pertain to the problem with a body orbiting another fixed body, so they are irrelevant. Sorry about that. The remarks below are still relevant.

The question is not whether μ = G(m+M). This is the definition of the standard gravitational parameter in the 2-body problem.

If your formula for the velocity is correct for your definition of μ, it is incorrect with the usual notation, and conversely.

I may at some point write a pdf about solving the 2-body problem from first principles, but right now I'm busy writing tests for the existing plugin.

DerpyHoovesX, if you want you can try to compile the plugin, instructions can be found here. They may be incorrect, in which case I would welcome feedback on where they fail.

Edited October 4, 2014 by eggrobin

[fibonatic]

EDIT:

These lecture notes from ETHZ indeed pertain to the problem with a body orbiting another fixed body, so they are irrelevant. Sorry about that. The remarks below are still relevant.

The question is not whether μ = G(m+M). This is the definition of the standard gravitational parameter in the 2-body problem.

If your formula for the velocity is correct for your definition of μ, it is incorrect with the usual notation, and conversely.

I may at some point write a pdf about solving the 2-body problem from first principles, but right now I'm busy writing tests for the existing plugin.

Thanks for replying to my posts. Every thing I know about Kepler orbits is mainly self-taught so I am not very familiar with the usual notation. I mainly base my calculations and equations of Newton's law of universal gravitation, Newton's second law and Kepler's laws.

By the way my posts where about finding a stable configuration of the three heaviest moons of Jool. A related question on stackexchange asks why three Galilean moons around Jupiter, with 4:2:1 resonance, are stable. It does not have any answers, but it does has a comment with a few upvotes which is linking to a document about this. I have only skimmed through this documents and did not see anything about weight ratios. Maybe this is because it is only about Galilean moons, so not the general 4:2:1 resonance situation.

[S1gmoid]

Hey,

Thanks for the hard work, and the scientific explanations. I'm really looking forward to this mod.

[eggrobin]

Status update:

Code has been written to plot the trajectories in nonrotating body-centred reference frames (and barycentric corotating, not reviewed yet). While abstractions will have to be written to do that more cleanly and efficiently, this yields some pretty pictures.

Caveat lector: While some of the following orbits may look very wild, bear in mind they are in non-inertial reference frames. Two of these wilder trajectories are followed by the same trajectory in a more pedestrian reference frame. In the Kerbin-centric non-rotating reference frame, they would all look like two or three elliptic orbits connected by strange segments

where the perturbation by the Mun is to blame.

The last picture below (in the barycentric corotating frame of the Kerbin-Mun system, that is, the unique frame in which the Kerbin-Mun barycentre as well as the line through Kerbin and the Mun are invariant) shows a trajectory that differs strongly from what stock KSP would have yielded: in stock this would have been a circular orbit around the Mun near the edge of its SoI. Instead it is a transfer to an eccentric Kerbin orbit, which is picked up by the Mun again after a while and ends with a crash on the Mun.

Javascript is disabled. View full album

With the addition of plotting, the C++ plugin has caught up with the features of the initial C# prototype. After the abstractions for plotting are written (we will add the remaining interesting reference frames in the process), the next step will be to take care of altering the trajectory when, e.g., engines are used.

It is apparent that the C++ plugin is significantly faster, much more reliable, and easier to debug than the C# prototype (it is also correct; some hasty shortcuts were made in the prototype that resulted in unreliable initial states and other problems).

Edited October 29, 2014 by eggrobin

[TMS]

Heh... our very own interactive Spirograph.

Looking good, eggrobin.

[AndreyATGB]

It's been a few months since I've tried it, at least a couple hundred commits (nice progress BTW) and even then it was displaying trajectories a lot better in the future but still not enough to be usable. Are you closer now to being able to show these relatively accurately for 1-2 years in the future, maybe more if used with RSS? Maybe you can use Kepler simulations for interplanetary sections since they shoulnt differ much from Newtonian ones.

[eggrobin]

It's been a few months since I've tried it, at least a couple hundred commits (nice progress BTW) and even then it was displaying trajectories a lot better in the future but still not enough to be usable. Are you closer now to being able to show these relatively accurately for 1-2 years in the future, maybe more if used with RSS? Maybe you can use Kepler simulations for interplanetary sections since they shoulnt differ much from Newtonian ones.

So, there are really two questions here.

When the vessel is not being perturbed by player-dependent actions (aka firing the engines), displaying predictions and prolonging them is not more costly than moving the vessels around.

When the vessel is being moved around by the player, well, at the moment that's not implemented yet. It seems however that it will be possible to recalculate predictions fast enough. Many possible optimisations have already been discussed in the thread, including FMM, caching the orbits of the celestials (caching would almost certainly massively help and isn't all that ridiculous in terms of memory usage), varying the timestep (we don't care as much about breaking symplecticity for predictions, and we only have one vessel we care about), etc.

[dreadicon]

Really glad to see work continuing on this. I still can't wait to try it out!

[John FX]

Glad to see this has not been forgotten. Hope there will be some useable fruit soon.

[Baughn]

How about threading? Would that be a possible optimisation, or are your calculations natively single-threaded?

[Kerbas_ad_astra]

Your introductory documents are excellent, but could you maybe recompile them so they don't get the curl connecting the "ct" and "st" ligatures? They're visually distracting, and they also disrupt text-searching within the document.

[eggrobin]

How about threading? Would that be a possible optimisation, or are your calculations natively single-threaded?

Parallelism may be possible, as discussed previously (see e.g. suggestions by SSR Kermit).

Your introductory documents are excellent, but could you maybe recompile them so they don't get the curl connecting the "ct" and "st" ligatures? They're visually distracting, and they also disrupt text-searching within the document.

But they're pretty!

If they're really distracting I may reconsider the ct and st "historical" ligatures, but search will never really work (while the fi ligature is searchable, Th isn't for instance).

[Kerbas_ad_astra]

But they're pretty!

If they're really distracting I may reconsider the ct and st "historical" ligatures, but search will never really work (while the fi ligature is searchable, Th isn't for instance).

Well, beauty is in the eye of the beholder. I wouldn't mind so much if there weren't so many words with "ct" and "st" in them -- and since we're talking about symplectically integrating system state vector functions, they stand out. Anyway, it seems that the "st" ligature is searchable (I'd tried searching for a word with "ct" and couldn't find any, so I'd assumed that "st" was the same way), and it seems to be the more common kind, so I guess I can make do if you insist on keeping the darned things.

[AndreyATGB]

Is there some key combination I'm missing? It used to display 2 windows one for selecting frame of reference and I can't really remember what the other did TBH. Now there's one which can start/stop plugin and 2 boxes next to all bodies. None of them seem to draw paths on the map anymore though, it does draw my past so to speak.

Damn building release is so much faster than debug, really impressive how little (almost none) performance impact this has.

Edited November 5, 2014 by AndreyATGB

[eggrobin]

Is there some key combination I'm missing? It used to display 2 windows one for selecting frame of reference and I can't really remember what the other did TBH. Now there's one which can start/stop plugin and 2 boxes next to all bodies. None of them seem to draw paths on the map anymore though, it does draw my past so to speak.

Oh, you built it? This means the setup instructions work then.

The 2-window version you're referring to is the C# prototype, which was last worked on about 6 months ago. In the meantime a lot of work was done to reimplement this in C++ with lots of abstractions (partly deal with the KSP API madness, partly because my C# code had got really messy on its own), and the plugin was reborn (with the actual work done in C++ called through P/Invoke) in the last update.

The C# prototype had 2 windows, they have now been fused into one. Some functionalities from the C# plugin has not been implemented yet, predictions are one of those (but predictions are only really interesting when intrinsic acceleration, a.k.a. rockets, works). Histories are drawn.

Damn building release is so much faster than debug, really impressive how little (almost none) performance impact this has.

Indeed, the C++ plugin is much faster than the C# prototype was (you need to build it on release, otherwise no optimization is done, and the abstractions get costly without inlining).

Status update:

Some work has been done on better abstractions for rendered trajectories, nothing significant yet (in particular there still is something in Õ(n) which could be in Õ(1), runs at each frame, where the constant factor involves evaluating trig functions, and n~10_000 is the number of rendered points, so that rendering trajectories in the barycentric rotating frame is slow), but some pretty pictures.

This shows a vessel near the L₅ point of the Kerbin-Mun system.

In the first image, the trajectory of the vessel is displayed in the nonrotating reference frame centred on the Mun, in the second it is displayed in the nonrotating reference frame centred on Kerbin, and in all subsequent images it is displayed in the barycentric corotating reference frame of the Kerbin-Mun system.

Javascript is disabled. View full album

Work on intrinsic acceleration is ongoing.

Edited November 7, 2014 by eggrobin

[fibonatic]

I am not sure what is meant by fixing the barycenters:

The barycenters will be fixed in order to stabilise the Jool system. Scott Manley (illectro)'s predictions indicate this will make it stable over at least 1_000 years (Pol was not part of those predictions, its stability remains to be determined).

If it is referring to viewing the orbital elements as seen from the barycenter, then according to my simulations this does not fix it. Or is it referring to fixing Jool onto the barycenter and thus Jool would not experience any forces from its moons? If physical parameters would have to be changed I would prefer not to change the difficulty relative to the stock game. Changing the mass of one or more moons would affect this, so I would prefer not to do that, however it is already impossible to land on Jool itself, so increasing its mass might be a possibility. Another possibility might be increasing the inclination of one or more of its moons, however this will affect the difficulty to transfer between them.

The main problem is that Vall has a resonant orbit with Laythe and Tylo and gets kicked out by them (mainly because its mass is an order of magnitude lower). These moons are analog for Io, Europa and Ganymede. However according to this sheet of orbital elements the phase difference does not match the phase difference for 1:2:4 resonance. To fix this I set the mean anomaly (at epoch) from 0.9 to 3.14-pi=0.0016. This in combination with orbital elements interpreted from the barycenter, as stated above, it not yet enough to make them stable.

My next attempt involved increasing the mass of Jool, also because this is an analog for Jupiter, which has relative to the sun and Earth a bigger mass than Jool has to Kerbol and Kerbin, this is namely of by a factor 4. However increasing the mass of Jool by a factor 3 already seems to make Vall stable for 63.4 years (2e9 seconds).

Another option would be increase the inclination of Vall. Due to the resonance Vall should always have close encounters with Laythe and Tylo in one spot in its orbit and with the corrected starting mean anomaly the close encounters with Laythe should be 180° away from those with Tylo. So by placing the ascending node of Vall between those two points of close encounters should allow me to increase the distance of both close encounters when increasing the inclination. However after increasing the inclination up to 90° still do did not give an stable configuration.

So concluding I would suggest that increasing the mass of Jool would be the way to go. If you want to keep the orbital periods the same as stock you could also increase the semi-major axes of all moons of Jool, this should not affect stability.

[eggrobin]

If it is referring to viewing the orbital elements as seen from the barycenter,

It is.

then according to my simulations this does not fix it.

According to Scott's it does (at least without Pol). We'll see which simulation was right when I get to that, at the moment sorting out the abstractions for the plugin and implementing intrinsic acceleration are much greater priorites.

The main problem is that Vall has a resonant orbit with Laythe and Tylo and gets kicked out by them (mainly because its mass is an order of magnitude lower). These moons are analog for Io, Europa and Ganymede. However according to this sheet of orbital elements the phase difference does not match the phase difference for 1:2:4 resonance.

Starting from real-world positions and velocities, the Jupiter system is (strongly) unstable if you assume Jupiter is spherical. You need to take the J₂ coefficient into account.

So concluding I would suggest that increasing the mass of Jool would be the way to go. If you want to keep the orbital periods the same as stock you could also increase the semi-major axes of all moons of Jool, this should not affect stability.

If the system turns out not to be stable, this is definitely a possibility.

[bs1110101]

Does anyone know what happened to the plugin that sorta did n body orbits and could timewarp with ion drives running? I lost all my mods to a reformat and i cant find it anymore, and i really need it with 6.4 scale kerbin and it's 40 minute long burn times for most nuclear ships and 2 hour+ burntimes for anything electrically propelled. Does that plugin still work, or was it killed by one of the updates and never patched, i really think that would be a shame.

[NathanKell]

eggrobin: awesome!

bs1110101: Orbit Manipulator

[Chase]

I am curious are you doing n-body with the planets as well? I am wondering if it would just be easier to keep those on rails.

[AndreyATGB]

I am curious are you doing n-body with the planets as well? I am wondering if it would just be easier to keep those on rails.

Go big or go home. I don't think on rails planets should be a thing anymore. This does n-body for all, Vall gets kicked out in a few days.