[WIP][1.8.1, 1.9.1, 1.10.1, 1.11.0–2, 1.12.2–5] Principia—version ‎‎Kronecker, released 2024-11-01—n-Body and Extended Body Gravitation

[undercoveryankee]

My point was that L-points are not very useful for communications satellites, which explains why we don't now have and never have had any communications satellite at real-life Lagrange points. Think of this: imagine you could put comms satellites in a stable orbit at the Kerbin-Mun Trojan points. These points lie nearly on the Mun's orbit, but 60° ahead of or behind the Mun. So, even with satellites at both points, there would still be a 60° region opposite Kerbin that is invisible to the satellites. The current RT solution of putting satellites on the Mun's orbit, but just outside its SOI works better because it covers more of the Mun's surface.

Lagrange points are only relevant for a particular two body system and only stationary within the rotating reference frame of that system. A Jool-Sun L4 point, for instance, would rotate around the Sun (or rather around the Jool-Sun barycenter) with Jool. It would not be stationary with respect to any of the other planets, so you couldn't use it as a reliable relay. RT2 has a couple of late-tech dishes that can reach from one side of Eeloo's orbit to the other, so nothing inside the kerbal solar system is out of their range.

SF writers like to talk about putting large colonies at Lagrange points or parking asteroids to mine there. Basically, since the real world has enough gravity sources to require station keeping, they're the places where station keeping is cheaper.

[CaptainKipard]

My point was that L-points are not very useful for communications satellites, which explains why we don't now have and never have had any communications satellite at real-life Lagrange points.

I kind of think that what explains it best is our lack of need for complete coverage. NASA rover drivers are absolutely fine with Earth-Mars conjunction blackouts

Think of this: imagine you could put comms satellites in a stable orbit at the Kerbin-Mun Trojan points. These points lie nearly on the Mun's orbit, but 60° ahead of or behind the Mun. So, even with satellites at both points, there would still be a 60° region opposite Kerbin that is invisible to the satellites. The current RT solution of putting satellites on the Mun's orbit, but just outside its SOI works better because it covers more of the Mun's surface.

Lagrange points are only relevant for a particular two body system and only stationary within the rotating reference frame of that system. A Jool-Sun L4 point, for instance, would rotate around the Sun (or rather around the Jool-Sun barycenter) with Jool. It would not be stationary with respect to any of the other planets, so you couldn't use it as a reliable relay. RT2 has a couple of late-tech dishes that can reach from one side of Eeloo's orbit to the other, so nothing inside the kerbal solar system is out of their range.

All of this is completely irrelevant. You're still missing the point. You can't direct a signal through the sun. AFAIK even RT2 doesn't allow it.

Edited July 30, 2014 by Cpt. Kipard

[Mr Shifty]

All of this is completely irrelevant. You're still missing the point. You can't direct a signal through the sun. AFAIK even RT2 doesn't allow it.

Ok yes. You're right, for Kerbin-Duna conjunctions (which occur once every 200 days and last for about 5 minutes) a Trojan relay would be useful. Most of the other planets are inclined enough relative to the kerbal system ecliptic that conjunctions will hardly ever be a problem. (Possible exception: Moho, since the sun takes up a larger chunk of its sky, and its synodic period relative to Kerbin is only 33 days.)

[CaptainKipard]

Ok yes. You're right, for Kerbin-Duna conjunctions (which occur once every 200 days and last for about 5 minutes)

I didn't realise that. This gives me an idea.

[Mr Shifty]

I didn't realise that. This gives me an idea.

Wow, you opened a can of worms there.

[eggrobin]

Stuff about Lagrange points

For some reason talk about N-body physics around here always end up revolving about Lagrange points, and more specifically to follow halo or Lissajoux orbits around them.

As said by Mr Shifty, Lagrange points are nothing more than the singular points of the potential of the CR3BP in the barycentric rotating reference frame in which that potential is time-independent (it's even worse for the ER3BP, since you need a rotating-pulsating reference frame). They are nice landmarks when considering 3-body trajectories (for instance some relatively stable trajectories revolve around them in Lissajoux or halo orbits), but they're not interesting per se.

Such orbits can be useful for observation (of the Sun and Earth for Sun-Earth L1, of the universe from Earth's shadow for Sun-Earth L2), but these concerns are sadly irrelevant in KSP.

There are however many other things that can be done with N-body gravitation, e.g., weak stability boundary trajectories (which also make ion engines far more interesting). I'm not sure why "orbiting" Lagrange points always steals the limelight.

It is rather fortunate that engines in timewarp (and thus sensible ion engines) and N-body gravitation, which go so well together, have the same basic requirements

contributing

Math is really not the problem, and is relatively easy to implement. Most of the work consists, as usual, in making up appropriate abstractions and implementing those. This API design work requires extensive discussions with pleroy, so contributing would be complicated. As things take shape it should become possible to contribute more conveniently (in any case, have fun reading the code, and feel free to ask any questions you have here or on IRC!).

I don't think I document that anywhere, but the styleguide used is the Google C++ Style Guide with a couple of amendments.

[CaptainKipard]

I'm not sure why "orbiting" Lagrange points always steals the limelight.

Maybe because hardly anyone has your experience.

Also partly because people repeat it. I only found out about them in this forum.

[Captain_Party]

Is this still in development?

[drtedastro]

Is this still in development?

Seems to be stuck, somewhere around L4 or L5 would be my guess...

[NathanKell]

Heh.

Or rather the fact that it's exam season.

[eggrobin]

Absolutely, I discussed some changes to NBodySystem and Body, and a new class, Trajectory, in order to have the requisite abstractions to keep several alternative trajectories (predicted, actual, manoeuvre nodes, etc.) with pleroy last weekend. Actual coding on my part has halted due to exams, and will resume in September. pleroy might find some time to implement these API changes, which I could then review.

Heh.

Or rather the fact that it's exam season.

NathanKell's cunning knows no bounds (and is therefore not compact).

[MeateaW]

as an aside; since you are going to be pinvoking an external library; to make things easiest on users; you could do what I did for my library.

Basically; KSP lazy loads the libraries; so you can execute code that chooses which precise library to load.

http://forum.kerbalspaceprogram.com/threads/87533-64bit-KSP-External-library-question?p=1302106#post1302106

I am not going to claim it is the best way to do it (it is certainly pretty ugly - and requires a lot of duplication when you have lots of functions in your library); but it certainly makes the effort of picking which version to install trivial for the user.

The only reason I would consider not doing this method; is if your library compiles to many megabytes, in which case to save download space you may not want to include multiple copies. But I suspect that this is not the case. (no embedded difficult to compress resources - unless there are some kind of large pre-computed lookup tables that can speed up the mathematics?)

[MeateaW]

And thinking about it even more; I have realised the ultimate way to do it, that simplifies even my code, and avoids unnecesary duplication.

On initialization of the plugin; copy the relevant library to a common name, and use the standard pinvoke calls. Due to the lazy-loading; you should be able to copy (or link?) the to the required file during plugin initialization.

[NERVAfan]

What is the current status of this?

[eggrobin]

What is the current status of this?

Things have resumed, I've started writing the C++ side of the actual plugin (as opposed to just the underlying abstractions and libraries).

[Guest]

Looking forward to it. Physicists rejoice! I'm downloading this the moment you have a playable version.

[Jsk2003]

Eggrobin, I absolutely love what you've done with this! I wish I knew where to begin going from multivariable calc and a small understanding of C to understanding symplectic integrators and coding them so I could help.

Curious, what are some of the roadblocks/problems?

[dreadicon]

Just wanted to mention that a multithreading plugin is in the works by ferram. I could see that being useful potentially. Or not, I noted the heavy use of C code, so this probably is effectively already in it's own thread apart from KSP anyhow. Just thought I would mention it anyhow.

[fibonatic]

The ability to look at orbits in rotating reference frames might actually make commsat management easier.

Yes, the same Matt Roesle (Mattasmack) also did a simulation of the Alternis system (he's mentioned in the 'acknowledgements' section).

Scott Manley remarked that the system might be much more stable if the orbital parameters were read as barycentric, rather than Jool-centric. Experiments will be carried out in due time.

I was playing with simulations of the solar system of KSP. With this I also tried to get Vall to stay. For this I initially tried to use barycentric orbital parameters, but it did not help. I think that the masses of the neighboring moons are just to massive compared to the mass of Jool. For instance in a real world example of 1:2:4 resonance with Jupiter's moons Ganymede, Europa and Io, their mass ratios are in the order of magnitude of 10000, while those of Jool's moons have ratios in the order of magnitude of 100.

[eggrobin]

I was playing with simulations of the solar system of KSP. With this I also tried to get Vall to stay. For this I initially tried to use barycentric orbital parameters, but it did not help. I think that the masses of the neighboring moons are just to massive compared to the mass of Jool. For instance in a real world example of 1:2:4 resonance with Jupiter's moons Ganymede, Europa and Io, their mass ratios are in the order of magnitude of 10000, while those of Jool's moons have ratios in the order of magnitude of 100.

This is interesting, since it contrasts with Scott's conclusions on the (pre-Pol) Jool system, cf. [OP::Further modding]. What kind of integrator did you use, and could you share your derivation and results for the initial conditions?

[...] a small understanding of C [...]

[...] I noted the heavy use of C code

C++11, not C (it's a good idea to specify 11, since C++11 is so different from C++03). C++11 and C are very different beasts.

so this probably is effectively already in it's own thread apart from KSP anyhow.

All the computationally intensive stuff will be done in C++.

[fibonatic]

This is interesting, since it contrasts with Scott's conclusions on the (pre-Pol) Jool system, cf. [OP::Further modding]. What kind of integrator did you use, and could you share your derivation and results for the initial conditions?

I implemented Störmer-Verlet method: Xi+1 = 2 * Xi - Xi-1 + A(Xi) * h^2

with a step size as small as 10 seconds in MATLAB.

I do have some about whether I implemented the barry-centric initial conditions correctly. I am not 100% sure if I even converted the orbital elements correctly in to positions and velocities.

But here it goes, I assumed that the semi-major axis is measured from the center of mass of the it and its primary body, so Jool in this case. This means that the actual distance between the moons and Jool is bigger than their semi-major axes, namely by a factor 1+m/M, where M is the mass of Jool and m is the mass of a moon. I incorporated this an effective gravitational parameters of Jool for each moon: GMeff=G*M^3/(M+m)^2.

I even tried to set their orbital periods to exactly 1:2:4 ratios with:

SMA_vall = SMA_laythe * (2 * (M + m_laythe) / (M + m_vall))^(2/3)

SMA_tylo = SMA_laythe * (4 * (M + m_laythe) / (M + m_tylo))^(2/3)

I positioned all moons relative to (0, 0, 0) after which I placed Jool in such a position that the center of mass of all moons and Jool is positioned at (0, 0, 0). When calculating the velocities I also used the effective gravitational parameters. When defining the volcity of Jool I did something similar as with its position, namely I set total momentum of the Jool system to zero. After all this I added the position and velocity of "Jool" relative to the sun to each body of the Jool system.

For the inital two positions I used:

X1 = X - V * h/2 + A(X) * h^2 / 8

X2 = X + V * h/2 + A(X) * h^2 / 8

The actual function which I used for the initial conditions in MATLAB can be found here.

Edited September 22, 2014 by fibonatic
added method of initial two positions

[eggrobin]

I implemented Störmer-Verlet method: Xi+1 = 2 * Xi - Xi-1 + A(Xi) * h^2

with a step size as small as 10 seconds in MATLAB.

Sounds good.

I do have some about whether I implemented the barry-centric initial conditions correctly.\

Remark: it's written barycentre (or barycenter in American), it comes from the Greek βάÃÂο (weight).

I am not 100% sure if I even converted the orbital elements correctly in to positions and velocities.

There is a typo in the comments labeling the columns of your table of elements, column 5 is the argument of periapsis, rather than the longitude of periapsis. It seems it is used as the argument of periapsis, so there's no mistake here.

There seems to be a mistake here:

        angle = A(i,5)*pi/180 + A(i,7); % Total angle away from ascending node

A(i,7) is the mean anomaly here, which is not an angle (it's more of a normalised area).

what you want is the corresponding true anomaly, which is directly computed from the eccentric anomaly, which in turn has to be computed iteratively from the mean anomaly.

Moreover, in the computation of the position,

pos(:,i) = A(i,2) * (1 - A(i,3)^2) / (1 + A(i,3) * cos(A(i,7))) * er

you should be taking the cosine of the true anomaly (angle) in the denominator here.

The velocity looks wrong too.

You could use the formulae from here, steps 1 trough 7 to have a correct conversion from orbital elements to Cartesian coordinates (in the computation of h, μ = G(M +m)).

Once you have computed these, compute as you said an initial state for Jool that cancels out the momentum and position of its satellites.

It's been a ridiculously long time since I've done one of these:

Status update

Trajectory, an abstraction for a tree of trajectories that can be forked into child alternatives (useful i.e. for predictions, manoeuvre nodes, but also simply computations at different precisions), has been written (part of the physics namespace)

The plugin has returned, first as a simple orbit-freezing plugin, then as a plugin that actually integrates the motion of vessels (only on-rails for the moment). Is does not plot yet, so it is still behind the C# prototype in terms of features, but it is significantly faster.

As was previously mentioned, all calculations are done in a native assembly, called by P/Invoke via a (cdecl) interface.

It should be noted that the plugin uses glog everywhere for logging (since logging from native code is needed).

Of course, this means Principia will have its own log files (in <KSP directory>/glog/Principia, with a log of the last session in <KSP directory>/stderr.log.

Another interesting aspect is that there will be no silent exceptions as in managed plugins, instead, the game will either segfault (when dereferencing an invalid pointer) or abort (SIGABRT) if a glog CHECK macro fails (on Windows the latter yields the "KSP.exe has stopped working" message).

I have been informed by Sarbian et alii that this will result in Fun support.

[dreadicon]

I'm still following this; the update news is quite exciting! It looks like this may actually help the game's computational load thanks to running in a separate thread (unless it does not and I am mistaken?).

I would try to contribute, but both the math and the C++ is over my head, i'm afraid. I do know a little about Planets in KSP as I have been working on understanding the Kopernicus mod's implementation so I can help with it's development. In order to control the planets, I'm not sure how deep you would have to go. what may be the easiest solution is to completely override them such that you are controlling all physics besides collision/ship based ones. Possibly even those. The problem I foresee with this is that KSP is very, very picky about Kerbin itself. Apparently parts of the game are hard-coded to look for a KSP planet object with the name of 'Kerbin', and the game kinda implodes if it cant find it. Once I finish my many, many side projects, I will start puzzling out a way to tell the KSP planetary physics systems as a whole to take a hike; but that'd be a while from now....

Also, as to global coordinates, I can foresee that being potentially difficult; a few versions back, Squad changed how the coordinates work for the realtime/rails scene rendering and physics making the active vessel the 'center of the universe', specifically to avoid floating point errors. Not sure what impact that has on your endeavors, just thought I would mention it.

There is a mod which predicts atmospheric trajectories recently released which predicts in-atmosphere, including for FAR; they actually borrowed the fixed body prediction code from wither this or the N-Body C# mod: http://forum.kerbalspaceprogram.com/threads/93685-0-24-2-Trajectories-v0-2-3-%282014-09-24%29-atmospheric-predictions-with-FAR-NEAR

Can't wait for a usable release; will help test at that time! It looks like you are really close; I'd ask how much longer, but such things are forbidden in the modding world So I shall just wish you Good luck!

Edited September 30, 2014 by dreadicon

[eggrobin]

control the planets

The C# prototype already did that, and the current C++ plugin does that now!

the easiest solution is to completely override them such that you are controlling all physics besides collision/ship based ones

This is the plan (with additional fixes to Krakensbane etc. in order to enforce conservation laws as well as we can).

Also, as to global coordinates, I can foresee that being potentially difficult; a few versions back, Squad changed how the coordinates work for the realtime/rails scene rendering and physics making the active vessel the 'center of the universe', specifically to avoid floating point errors. Not sure what impact that has on your endeavors, just thought I would mention it.

Yes, I have started to gain a better understanding of that insanity. A couple of these open questions are no longer relevant actually.

There is a mod which predicts atmospheric trajectories recently released which predicts in-atmosphere, including for FAR; they actually borrowed the fixed body prediction code from wither this or the N-Body C# mod: http://forum.kerbalspaceprogram.com/threads/93685-0-24-2-Trajectories-v0-2-3-%282014-09-24%29-atmospheric-predictions-with-FAR-NEAR

I saw that, it's neat.

[fibonatic]

There is a typo in the comments labeling the columns of your table of elements, column 5 is the argument of periapsis, rather than the longitude of periapsis. It seems it is used as the argument of periapsis, so there's no mistake here.

I copied to names of the elements from this google document, which uses that name for that column, without thinking. The top name actually does state "Argument of periapsis", however its abbreviation, LPE, seems to suggest longitude of periapsis.

A(i,7) is the mean anomaly here, which is not an angle (it's more of a normalised area).

You are totally correct, I often switch those two (true/mean), however since most eccentricities are low (or even zero) and most starting mean anomalies are close to zero or PI, then the true anomaly will not be that different or in some cases equal to the mean anomaly. I have implemented the correct anomaly into my script now.

Moreover, in the computation of the position,

pos(:,i) = A(i,2) * (1 - A(i,3)^2) / (1 + A(i,3) * cos(A(i,7))) * er

you should be taking the cosine of the true anomaly (angle) in the denominator here.

Are you still referring to the mistake between true and mean anomaly, since cos(A(i,7)) is in the denominator?

The velocity looks wrong too.

I used these equations I derived a while ago, which I think are correct. These equations can be simplified down to this:

Again I used the mean anomaly instead of the true anomaly, however I did this consistently, so only the starting true anomalies of a few celestial bodies where off, but the other orbital parameters should be correctly translated.

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, however I do use a different definition for μ.

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.

Edited October 3, 2014 by fibonatic