eggrobin

[WIP][1.1.3] Principia - version Cantor, released 2016-08-13 - N-Body and Extended Body Gravitation

631 posts in this topic

I was asked to make a dev thread for this mod (which is currently very much a work in progress, and can be rather crashy), so here it is.

The source code is available on GitHub under the MIT license.

Abstract

This mod aims to replace KSP's unstable physics integration with a higher-order symplectic integrator, adding N-body Newtonian gravitation in the process. Once N-body gravitation with point masses is made playable (including faithful trajectory prediction over long periods and reasonably accurate dynamically-updated trajectory predictions over a few years for vessels under thrust), further interesting physical effects can be added.

Acronyms

Spoiler
  • ACS: Attitude Control System
  • AMD: Advanced Micro Devices
  • AMP: Accelerated Massive Parallelism
  • ARG: Ada Rapporteur Group
  • API: Application Programming Interface
  • BIPM: Bureau International des Poids et Mesures
  • CR3BP: Circular-Restricted 3-Body Problem
  • CGPM: Conférence Générale des Poids et Mesures
  • CIPM: Comité International des Poids et Mesures
  • CIL: Common Intermediate Language
  • CLI: Common Language Infrastructure
  • CLR: Common Language Runtime
  • CLS: Common Language Specification
  • CPU: Central Processing Unit
  • CRTP: Curiously Recurring Template Pattern
  • CTS: Common Type System
  • DFT: Discrete Fourier Transform
  • DLL: Dynamic-Link Library
  • ER3BP: Elliptic-Restricted 3-Body Problem
  • ETHZ: Eidgenössische Technische Hochschule Zürich
  • EVA: Extra-Vehicular Activity
  • FAQ: Frequently Asked Questions
  • FAR: Ferram Aerospace Research
  • FFT: Fast Fourier Transform
  • FMM: Fast Multipole Method
  • FPU: Floating-Point Unit
  • FSAL: First Same As Last
  • GPGPU: General-Purpose Computing on Graphics Processing Units
  • GPU: Graphics Processing Unit
  • GR: General Relativity
  • ICC: Intel C++ Compiler
  • IDE: Integrated Development Environment
  • IEEE: Institute of Electrical and Electronics Engineers
  • IMCCE: Institut de Mécanique Céleste et de Calcul des Éphémérides
  • IRC: Internet Relay Chat
  • JPL: Jet Propulsion Laboratory
  • KSO: Kerbin-Synchronous Orbit
  • KSP: Kerbal Space Program
  • LKO: Low Kerbin Orbit
  • MET: Mission Elapsed Time
  • MIT: Massachusetts Institute of Technology
  • NASA: National Aeronautics and Space Administration
  • OBE: Overcome By Events
  • ODE: Ordinary Differential Equation
  • OP: Original Post
  • PDF: Portable Document Format
  • PFCE: PlanetFactory Creator Edition
  • PGO: Profile-Guided Optimization
  • PSA: Public Service Announcement
  • RCS: Reaction Control System
  • RHS: Right-Hand Side
  • RK: Runge-Kutta
  • RSS: Real Solar System
  • RT2: RemoteTech 2
  • SAS: Stability Augmentation System
  • SI: Système International
  • SIMD: Single Instruction, Multiple Data
  • SOI: Sphere Of Influence
  • SPRK: Symplectic Partitioned Runge-Kutta
  • SSE: Streaming SIMD Extensions
  • VB: Visual Basic
  • VS2013: Visual Studio 2013
  • WIP: Work In Progress

Notations

Spoiler

 

B: number of massive bodies.

V: number of massless bodies.

 

Terminology

Spoiler

 

The main integration refers to the integration of celestial dynamics in the absence of thrust, both for the propagation of the vessels and for the predicted trajectories of passive vessels and planets. It also includes the real-time integration of vessels under thrust and, when this is added, the integration and trajectory prediction for planned burns under timewarp (e.g. ion engines).

The dynamic integration, or integration for vessels under thrust, refers to the predicted trajectory of a vessel assuming thrust is cut immediately (what happens is stock KSP when you fire your engines). This integration has to be very fast, as the initial conditions change quickly and prediction over several years may be required for planetary transfers &c.

Fictitious forces in any reference frame are those proportional to the mass of the object upon which they act. This includes gravity. This convention only really makes sense in GR, which I do not currently intend to implement, but it's more convenient.

The geometric acceleration is the acceleration due to fictitious forces.

The proper acceleration is the acceleration not due to fictitious forces.

 

Documentation and download

Principia is not built in the same way as most other mods (most of it is native code, with a thin C# interface layer), so that compatibility issues may arise depending on the platform.

Binaries are available for Windows and Linux (build on Ubuntu, your mileage may vary when using these on other distributions). Note that the mod only works with 64-bit versions of KSP. We do not currently support OS X (macOS, Mac OS X, or whatever its current name is).

Please read the wiki page detailing the concepts carefully; a lot of things behave differently when using principia, in particular map view and flight planning. Many aspects of the mod are unfinished, and this is detailed on that page.

In addition, please read the FAQ and bug reporting guide. Note that bug reporting procedures differ significantly from those of other mods because of our custom logging system, our use of fatal failures, and our journalling system (which, if activated, allows us to replay the events that led to a crash for later debugging).

A tutorial is available on the topic of going to the Mun in various ways, applying the concepts described above.

Since Principia makes all bodies attract each other according to Newtonian gravitation, the stability of the solar system is a concern; when using stock KSP, Principia modifies the Jool system so that it is stable. When using a customized system (e.g. with Kopernicus), you need to make sure that the system is stable.

If you wish to ask questions or report bugs, you may want to do so on IRC (on the #principia IRC channel on EsperNet). Your concerns may be addressed more quickly there than on the fora (bear in mind however that there isn't always someone around who can answer your question right away; if no-one is around, you may want to stay on the channel for a while).

Download the binary (Ubuntu and Windows) here or, if you don't trust our binary, build the mod from the Cantor release.

Status

Dates are ISO 8601 extended format.

2016-08-13

Cantor is out.

This version is mostly the same as بوژگانی (see the change log), the intent is to make sure it is as stable as بوژگانی was: this is the first broadly available release, with a link to binaries outside the IRC channel (see the OP).

In the meantime, progress has been ongoing (note that this is about a future release, Cantor contains none of the features below).

Regarding our internal libraries, we have clarified the way we handle time (which is always an extremely confusing thing), with our Instant type aligned with Temps Terrestre (TT), and support for date literals in TT, Temps Atomique International (TAI), Coordinated Universal Time (UTC), and UT1 (based on the rotation of the Earth). This allows us to specify dates more cleanly in our tests.

We tested the difference between the orbits of Mercury predicted by our ephemeris and by JPL's over a century, and we find the expected error in perihelion precession due to general relativity.

On the side of flashier features, we seem to be well underway to having working axial tilt, by tilting universe the current main body (with terrain) is aligned with the axis used by KSP, and rotating the other ScaledSpace models appropriately; in a future release of Principia, the Jupiter, Saturn, and Uranus systems should look much nicer, with the planet properly aligned with its satellites.

2016-06-26

بوژگانی is out.

A change log can be found here.

2016-05-20

Burnside is out.

A change log can be found here.

2016-05-05

Буняковский is out.

A change log can be found here.

2016-02-24

There is a bug in Buffon (2016022220-Buffon-0-g0455d0cb0e4b0b584a84caf40520cf7993903e0c) that causes a crash when starting a new save with RSS (loading an existing save works fine). The hotfix "2016022220-Buffon-12-g1298cffba22d90119bd59e9933a9ef261922a423" (let's call it "Buffon + 12") resolves that, so if you run into this issue just go back to the IRC channel to get the latest build, it will be Buffon + 12.

There are no other changes between Buffon and Buffon + 12.

2016-02-22

Buffon is out, you can get it by asking on IRC (#principia on EsperNet), as usual.

Changes:

  • The integrators now limit the number of steps they perform, and terminate if their step size vanishes. This avoids issues where the plugin would hang when the trajectory would accidentally get very close to the centre of a celestial body or spend a long time in a low orbit.

  • A use-after-free bug has been fixed which caused a variety of crashes (#872, #881, #889, #896) when the historical trajectory was shortened in a way that would cause it to start after the beginning of the flight plan.

  • The version identifier of the plugin is now displayed in the UI to make it is easier to assert what version is running.

  • A verbosity option has been added to the journalling which makes it easier for us to reproduce crashes.

The first two items above are illustrated by the following two reports.

On 2/12/2016 at 5:19 AM, nilof said:

Great work. Every time I've tried to do that, my flight planner crashes the game because I plotted a trajectory that intersects earth. :D

On 2/12/2016 at 11:29 AM, maccollo said:

One thing that will crash it is if your trajectory history becomes less than the starting point of your flight plan.

For more details see all 19 pull requests between Brouwer and Buffon.

2016-02-11

Some clarifications regarding reference frames and flight planning.

2016-02-09

Brouwer is out, you can get it by asking on IRC (#principia on EsperNet), as usual.

User-facing features:

  • The whole Frenet trihedron is now displayed in the correct reference frame when "fix navball in plotting frame" is selected.
  • The initial state (and thus the evolution) of the system is now deterministic even when not using RSS.
  • Tidally locked bodies no longer spin back and forth madly (on the other hand, they may not be tidally locked if their mean period differs from their Jacobi osculating period).
  • When using stock, the Jool system is modified, cancelling the apocalypse. Specifically, we make the inner Jool system nonresonant, since we have been unable to replicate the results (Manley, priv. comm.) according to which some interpretations of the orbital elements yielded a stable Laplace resonance, despite systematic searches of the Jacobi osculating elements. In addition, at Scott Manley (@illectro)'s and @UmbralRaptor's suggestion, we put Bop in a surprisingly stable, though highly precessing, retrograde orbit. The modified system is stable for upwards of a century.
  • Flight planning has been implemented.

Modder-facing changes:

  • When a Cartesian initial state cfg is not given, the KSP orbital elements are interpreted in a hierarchical osculating Jacobi fashion; for instance, the orbital elements of Jool are the osculating elements at game start of the orbit of the barycentre of the Jool system around the barycentre of the (sun, moho, eve, gilly, kerbin, mun, minmus, duna, ike, dres) system; the elements of Laythe are the osculating elements at game start of the orbit of Laythe around Jool; the elements of Vall are the osculating elements at game start of the orbit of Vall around the Laythe-Jool barycentre.

Optimizations:

  • The Windows build now uses profile-guided optimization (we estimate that this improves performance by ~20%); in theory this could be extended to other platforms.
  • The evaluation of the Чебышёв series has been significantly optimized.
  • @sarbian made trajectory rendering faster (as he pointed out, there is still lots of room for improvement).

Other features:

  • Everything that crosses the interface can now be journalled if the right flag is set, allowing us to replay the C++ side of a session; this is useful for tracking down tricky bugs, and it enables profile-guided optimization.

Highlights of miscellaneous library changes; beware, this gets technical:

For more details, see all 195 pull requests between Bourbaki and Brouwer.

2015-08-16

Bourbaki is out, you can get it by asking on IRC (#principia on EsperNet), as usual.

Please bear in mind that the channel operators may be away from keyboard, so you should wait until you're noticed (it also helps to say the name of the channel operators since that pings them).

As of 2015-08-16T15:38Z, the build for Win32 is ready, we're waiting for Norgg and armed_troop for the Linux and Macintosh builds respectively.

Note that you should read the FAQ before installing principia.

Rough changelog:

 

  • Reimplemented integrators: the symplectic Runge-Kutta-Nyström integrator was reimplemented more cleanly, an embedded explicit Runge-Kutta-Nyström integrator was implemented. Abstractions for differential equations were created.
  • Ephemeris: the celestial bodies are integrated on their own, with their own (much larger) timestep (45 min);
    their trajectories are then fitted using чебышёв series, yielding continuous trajectories. This means that when there are no vessels (including asteroids, see the FAQ), timewarp at very high speed (6'000'000x was tested in RSS) is smooth.
  • The predictions are computed using an adaptive step size, so they're faster and less fiddly (you still get a tolerance setting, but it doesn't need as much attention as the step size setting; the step size will shorten near periapsis and lengthen near apoapsis on its own). The histories are still in fixed steps of 10 s, that will likely change in Brouwer, since it is one of the biggest performance costs now.
  • There is an initial configuration for Real Solar System: the planets will start in the right places as given by the JPL HORIZONS service, and they are given gravity models using the freshest data available (Vesta's model is from Dawn data, some Cassini data gets used). Please note the RSS-specific recommendations in the FAQ.
    A side effect of that is that the moon becomes far more accurate: since the motion of the moon is very much a 3-body problem, it cannot be accurately represented in RSS alone. In particular, real-life eclipses can be observed in principia + RSS (please note the unexpected inaccuracy mentioned in the FAQ).
  • This initial configuration also includes J2 for the Sun, the planets, the Moon, and Vesta, so the resulting effects are felt (precession of Earth orbits, the possibility of heliosynchronous orbits, etc.).

 

Bourbaki is save-compatible with Borel, for RSS users, please note that unless you reset the plugin, the new initial state and gravity model configuration files will not be taken into account.

2015-05-07

The new version, Borel, is out (you can get it on IRC, as previously; the same caveats apply).

Rough changelog:

 

  • ported to 1.0.x (that only took a couple of lines);
  • custom navball settings, so that the navball can be fixed in the reference frame in which the trajectory is plotted; the IVA navball is unaffected;
  • when using the custom navball, the prograde/retrograde vectors are now in the correct reference frame, consistent with what is shown is map view; the rest of the Frenet trihedron (the radial and normal vectors) are unaffected at the moment and will be fixed in another version;
  • fixed a few bugs (the tetrydsbug, the melonbug, the thutmosebug);
  • less memory-hungry;
  • added a setting to clip histories;
  • added a toolbar button to show/hide the UI;
  • the UI now acknowledges the F2 key;
  • other things that I forgot.

 

Moreover, Norgg and armed_troop have made good progress on Linux 64-bit and Macintosh 32-bit builds respectively, so if you're on these platforms you can go on IRC and ask for a build or for build instructions.

With Windows 32-bit, Linux 64-bit and Macintosh 32-bit, we should be covering most users (I don't think it's worth doing a Linux 32-bit build).

2015-03-22

Serialization has been implemented, and rudimentary predictions have been added.

The predictions are currently using the same integration method (McLachlan and Atela's 1992 optimal 5th order method), with the same splitting of the Hamiltonian (kinetic + potential energy), this is somewhat usable but unacceptably slow.

I am currently implementing as many symplectic integrators as I can find in order to compare them, and I will also be comparing various splittings (as as nonsymplectic methods in use for computing ephemerides). I will probably write a post about these at some point (probably an advanced one, rather than an elementary introduction, this is quite a complicated topic).

It seems that there is some interest in builds of Principia, no matter how unstable or slow they may be.

If you want to try out the current version of Principia (Bolzano) for the Windows 32-bit version of KSP, go the IRC channel linked below and ask me for a build (my nickname there is egg, or sometimes something of the form egg|<status>|egg).

CAVEAT:
The current version is not stable, it can corrupt or otherwise destroy saves, it has been known to crash, and predictions are horrendously slow. #principia IRC channel on EsperNet.

2015-01-22

The refactoring of the physics bubble has been done (and some tests have been added), as well as some cleanup in the C# adapter (including some performance improvements, a lot of the performance issues occur on the C# side since the implementation is a bit careless there).

More refactoring occurred (use of a not_null template for non-null pointers).

The next step on the path towards playability is persistence: the state of the plugin must be persisted so that we remember the states of the planets, the histories of the vessel trajectories, etc. (currently loading a save or reverting will result in either a crash or a loss of consistency due to the planets having moved around).

Persistence of our types will be achieved using protocol buffers, stored in config nodes: we don't want to use KSP's config format directly, since we have lots of highly structured data and operate far from KSP's API it would be inefficient and clumsy.

We'll probably also use protobufs for caching when implementing trajectory predictions later on.

Some work has already been done in that direction.

Here are some more screenshots, showing a flight from Kerbin to the Kerbin-Mun L4 Lagrange point (the reference frame for each screenshot is indicated in the "Traces of various description" window).

Javascript is disabled. View full album

antedated 2014-12-13

Support for intrinsic acceleration has been added. Refactoring will be needed, but it works, as demonstrated by the following plane change manœuvre:

4n4ax6el.png

2014-11-07

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

2014-10-28

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

2014-09-30

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. :P

2014-07-04

I wrote the second explanatory post, on Hamiltonian mechanics (PDF).

Prerequisites are chapters 4, 8, and sections 11-4 and 11-5 of Richard Feynman's Lectures on Physics.

The next post will be on symplecticity and the leapfrog integrator.

2014-06-27

The integrator (still a 5th order SPRK, same as in the C# prototype, the Saha & Tremaine stuff isn't here yet) has been implemented, tested and benchmarked (using a Windows port of google/benchmark) in C++, as well as the first level of abstraction above it, principia::physics::NBodySystem (header, body, test, header for test data, test data, test of test data (which is also a nice example of the use of principia::quantities and principia::geometry, benchmark (using the test data)).

A significant change of plans has occured: As Unity uses the .NET Framework v2.0, it is not possible to use plugins compiled for the Framework v4.5. This means C++/CLI projects need to be compiled using the platform toolset v90 (from VS2008), rather than VS2013's v120. Of course v90 does not support C++11 nor C++14, so it is not an option for this code that strongly depends on C++11 features.

As a result, we will keep using the C++11 codebase, compiling it to a native DLL, and using P/Invoke to call it from a C# adapter. The C# adapter will perform no calculations, only transmit data from KSP to the native plugin and apply the changes/render the trajectories returned by the plugin. A simple test P/Invoke plugin can be seen on the repository (header, body, C# adapter) and works in KSP.

This means there will be separate x86 and AMD64 builds for each target operating system (Microsoft Windows and 57 varieties of Unix), possibly more if we decide to do specific optimisations for Intel chips (though ICC is not cheap).

Moreover, this will allow a switch to clang in the near future, so that we can have saner error messages and better compliance with the standard.

2014-05-12

I now have a collaborator on Principia, https://github.com/pleroy (my father). This means the code gets reviewed and that development is faster.

We have decided to completely switch to C++/CLI due to better test tools and useful language features. Most of the code will actually be written in standard C++, with the implementation inlined in header files, so that they are seen by the eventual C++/CLI managed assembly (If we were to use C++/CLI everywhere, we would need managed boundaries between the assemblies, which is quite inconvenient).

As an added advantage, using standard C++ enables putting performance-critical parts into a native assembly if needed at a later time, without requiring a significant rework of the code.

We have started switching to gmock, gtest and glog for mocking, testing and logging. These tools are more convenient and powerful than the Visual Studio testing framework and are open source, so that users of Visual Studio Express (which does not support Microsoft's testing framework) will be able to build and run tests if they want to.

Here is the first test to be converted to gtest: https://github.com/mockingbirdnest/Principia/blob/master/geometry/sign_test.cpp.

2014-04-03

The Quantities library (C++) is pretty much done and tested.

I am currently porting the C# Geometry assembly mentioned previously (now called CTSGeometry) to C++ in order to enable its use with these quantities.

2014-03-15

The Geometry assembly seems to be mostly feature complete, so I'll start writing tests tomorrow (Saturday) after changing a few access modifiers, replacing the ugly switch statements in Permutation by something nicer, and a few optimisations.

See here for an overview of its features.

2014-03-04

I have investigated the feasibility of using other languages for KSP plugins:

The following languages can be used on their own to write a plugin: C#(unsurprisingly), F#, C++/CLI (with no calls to native code).

VB.NET can be used, but wrappers in one of the above languages are needed due to its case-insensitivity.

Unity does not support mixed-mode DLLs, so calls to native code have to be made using DllImports in one of the above languages.

I have started doing some refactoring since the sphaghettiness of the code was getting on my nerves.

I have written strongly typed wrappers for the numerous reference frames spawned by KSP (direct vs. indirect, rotating vs. inertial, world vs. body-centric, etc.) as I have had numerous bugs due to a misplaced xzy, rotation, translation, scaling, inertial force etc.

Of course, since KSP has the brilliant idea of using both direct and indirect reference frames, I needed distinct types for vectors, bivectors and trivectors (basically I had to strongly type Grassmann algebras; there can be no cross product, only wedge products, commutators and left and right actions of alternating forms---where ncnhkgd.png is identified with pynpdhk.png through the inner product, and ogat7w6.png is identified with qcv7ksx.png).

I do not think I will implement strong typing for physical quantities yet (though I'd like to), since C# generics are not powerful enough for that; I would need C++ templates. I'll do that when I rewrite things in C++/CLI.

The next step is to implement my own versor, since Unity's Quaternion is in single-precision float and KSP's QuaternionD is broken.

The rest should be more straightforward refactoring.

2014-02-24

It turns out I have trouble properly setting the position when off-rails (finding out where the reference frame actually is is hard).

However, there are some nicer news: I seem to have found the source of the phantom acceleration bias (not the ones arising from rotation, the bias that is removed by timewarping). The floating origin sometimes floats several km away from the ship, so that's probably just floating-point inaccuracies (the usual Kraken). If that hypothesis turns out to be true, this particular acceleration bias will be easy to fix, just reset the floating origin often enough.

2014-02-19

I have successfully implemented the integration of proper acceleration for the active vessel when off-rails.

Further experimentation on proper acceleration has led me to the following conclusions:

  • There is a bias in proper acceleration coming from some improperly initialised variable in KSP. Indeed, when loading a vessel in LKO, I observe a strong bias in proper acceleration (~20 mm s-2). This bias is observed independently of the way proper acceleration is computed (differentiating position twice, differentiating any of the numerous velocities, etc.) and geometric accelerations have been checked from first principles (the difference in geometric acceleration depending on the point of the vessel it is computed at is negligible). The bias is reduced to below 1 mm s-2 when warping then coming out of warp. It should be noted that the strong bias is not seen in Vessel.perturbation, but Vessel.perturbation consistently has a bias of 4 mm s-2. As I have attempted to compute the proper acceleration in many different ways and all were consistent with each other and inconsistent with Vessel.perturbation, I assume Vessel.perturbation is mostly nonsense.
  • Accelerations below 1 mm s-2 are biased, unavoidable, unusable in stock, and should be clamped to 0. The acceleration from low-thrust engines will have to be fetched by hand.
  • It had previously been mentioned that spinning the ship accelerates it. If spinning the ship with angular velocity É produces a phantom acceleration a, then spinning it with angular velocity -É produces a phantom acceleration -a. It seems a is either prograde or retrograde.

 

2014-02-17

I have finally managed to work on this a bit over the weekend.

Extensive experimentation has failed to yield a better velocity, so we are stuck with computing the proper acceleration from the rather mysterious Vessel.obt_velocity for now (for some reason in whatever reference frame this velocity is in, the geometric accelerations are as expected, except for the Coriolis acceleration which is halved :confused:).

This yields roughly the same results as Vessel.perturbation without the moving average.

The same experiments have further shown that a naïve low-pass filter will not be enough to remove Unity's silliness. For instance, one finds the proper acceleration is strongly correlated with the rotation rate of the ship. Smarter post-processing will be needed.

I shall now attempt to actually integrate whatever proper acceleration I have managed to grab.

I also wrote the first explanatory post, which describes ODEs and Runge-Kutta methods (PDF).

I have tried not only to explain how Runge-Kutta methods work (this can easily be found on Wikipedia) but why they have this structure.

Hopefully you will find it interesting.

2014-02-08

I fixed a bug or two since 2014-02-05, added a few TODOs, but have not cleaned up the code to any measurable extent. I am, however, publishing it on GitHub (under the MIT license).

Caveat Compilator: As previously stated, this is just a proof of concept with a bunch of traces. Pressing the wrong buttons in the wrong order will result in lots of NullReferenceExceptions and off-rails crafts will behave weirdly. Using HyperEdit to set up initial conditions, then switching to Newtonian physics and fooling around with reference frames can be entertaining though.

You might learn something from looking at the Integrators part (which admittedly contains only one integrator), the rest was quickly hacked together, has no tests, is badly structured &c.

The 'Simulator' project probably won't compile, it's leftover from my Alternis Kerbol simulations. It's not needed.

2014-02-05

This is currently hardly more than a proof of concept. The simulation only affects on-rails vessels, thrust is ignored (thus you can't actually play while the simulation is running) and the integrator is too slow. However, it shows a few interesting things:

 

  • Near-unit-roundoff (using IEEE 754 binary64) errors can be achieved while sustaining 100_000x timewarp with B = 17, V < 5. Handwavy asymptotic calculations indicate that for V = 100, this integrator would slow down to about 20_000x timewarp. Using better integrators (and saner tolerances), performance could easily be increased by a couple of orders of magnitude, making the simulation usable for RSS as well as stock.
  • Trajectory plotting in rotating reference frames, or reference frames centered around a body other than the primary, can be useful even for near-Keplerian orbits: it allows for visualisation of RT2 satellite coverage, easier landing prediction, etc.
  • The stock integrator is terrible, and will have to be overriden at all times while in space (when near a Lagrange point, a fraction of a second under stock integration will turn a Kerbin capture into a Mun collision).
  • If you think floating SOIs (or worse yet, SOIs with repulsive gravity) are a half-decent model for Lagrange points, you don't really understand how Lagrange points work.

 

Expect a slow development cycle, due to a combination of laziness, exams, and this actually being a complicated project.

I'll put the source on GitHub as soon as I can clean things up a bit (there are quite a few traces that were too hastily implemented for my taste). This might take a few days, see above.

Nonetheless, here are some pictures that illustrate the fun trajectories in the N-body problem. The second and third pictures have a misplaced trajectory, this bug has been fixed. Some pictures have three strange red, green and blue lines, which indicate the origin and axes of world space.

Javascript is disabled. View full album

Considerations on numerical integration

The current prototype uses a straightforward fifth order SPRK method. The coefficients used are from (Atela & McLachlan 1992). The increments are accumulated using Kahan summation. The integration is performed with a constant 10 second timestep (numerical experiments suggest that the error is close to an unit roundoff with a 5 second timestep for LKO).

Regarding implementation details, the integrator is written in C# (compiled with VS2013) and uses double (IEEE 754 binary64) for all computations.

Scott Manley suggested using hybrid symplectic integrators (Chambers 1999) in order to speed things up. It should be noted that the concerns that the Chambers paper attempt to alleviate

Chambers 1999 said:

The fixed time-step inherent in symplectic algorithms makes dealing with close encounters particularly difficult. Ideally one would like the step-size to decrease during an encounter in order to preserve the accuracy of the overall integration. However, changing the step-size of a symplectic integrator introduces an error with each change. If close encounters do not occur too often, this technique can still be used when the results of the integration are interpreted in a statistical sense (Levison & Duncan 1994). However, such an integrator cannot be relied upon to reproduce the true orbital evolution of a particular body (Michel & Valsecchi 1997).

might not be critical for gameplay purposes (the main integration can probably be done with a very small timestep compared to the ones commonly used in astronomy. Moreover, timestep reduction occurs de facto when getting out of warp, as the game cannot conceivably be played with timesteps of several seconds).

It could become relevant for trajectory predictions in vessels under thrust, which has to be done in a few hundred milliseconds over several years. For those predictions, the timestep would have to be increased to durations comparable to those used in the paper, and similar concerns would arise.

The results from the papers quoted by Chambers, namely (Saha & Tremaine 1994) and (Holman & Wisdom 1991), will probably have to be used in order to make the main integrator faster. Experiments will have to be carried out in order to find out whether a higher order (> 2) integrator is worth using, and how the added cost of Keplerian evolution (which requires numerically solving the Kepler equation for all bodies, which is very costly due to the trigonometric functions) affects performance at acceptable time steps.

Open questions for the interested reader

Regarding KSP modding

The method I currently use for drawing trajectories (LineRenderers, object in layer 9, transform.parent = null, useWorldSpace=true, updated every time the GUI is drawn) has some issues: when rotating the camera in map mode, the lines lag behind. Does this have a known solution? (I guess so, RT2 doesn't have this problem). And indeed, the RT2 code contained the solution. Thanks, Cilph!

Is there a way to directly set the position and velocity of the planets and vessels in world coordinates? OBE, the API does not make sense in highly original ways.

Can anybody think of a way to mod in axial tilt (even a hacky one; remember that I'm the one moving the planets around anyway)? OBE (done in RSS, we actually have axial tilt, it's just that all planets rotate around parallel axes).

Regarding aerodynamics

Could somebody help me with the implementation of orbital decay drag when I get to implementing that?

ferram4 told me he'd help. Thanks! :)

Would it be conceivable to predict the results of aerobraking/aerocapture/ballistic reentry with FAR (within some error margins, to be displayed on the predicted trajectory)? How slow would that be?

Answered by ferram4.

This might end up happening.

Regarding astrophysics and astrodynamics

The uniformly rotating reference frame around the barycenter with respect to which the two massive bodies are fixed is very useful when looking at trajectories in the CR3BP (Lagrange 1772). What would a good reference frame for the ER3BP (for bodies with highly eccentric orbits, e.g., Eeloo in stock, Pluto in RSS) be? the uniformly rotating one around the barycenter that follows the mean anomalies, the nonuniformly rotating one that follows the true anomalies, or something else entirely?

In the rotating-pulsating reference frame that fixes both bodies, the equilibrium points are fixed. Is there a way to draw things in that reference frame on the in-game map that's not too confusing? Can it still be useful to draw the potential for the parts of the equations of motion that derive from a potential, given that the independent variable is true anomaly rather than time?

 

Should long-term trajectory predictions for vessels under thrust be inaccurate, is there a good way of estimating the uncertainty (for instance, would predicting a few perturbed trajectories and plotting them all give a representative idea of what might happen)? ferram4 suggested a few interesting things.

Would drawing the gravitational potential in the current orbital plane (together with the centrifugal potential if orbits are drawn is a rotating reference frame) be useful? Answered by ferram4 with visualisation suggestions. The answer is 'yes'.

How should stationkeeping be implemented? In particular, how should the user set the orbit they want to maintain?

Regarding the KSP fora

Do we have LATEX support? If that isn't the case, is it conceivable to have it in the foreseeable future?

Further modding

Once the mod gets to a functional state (taking thrust into account), I shall attempt to solve the problem of trajectory prediction for vessels under thrust. In the process, the smarter integrators mentioned above will be implemented, and the most efficient one will be kept. This should significantly increase computation speed.

This will also enable instantaneous-burn maneuver nodes (the kind of maneuver node we have now).

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

Speed will be further increased by rewriting the numerical integration in (unmanaged) C++, interfaced with C# through C++/CLI. The unmanaged C++ (compiled with clang) will use 80-bit extended precision 'long double' (64-bit mantissa instead of 53-bit) for added precision.

Once the mod becomes playable with N-body gravitation with point masses, the following effects should be relatively easy to add:

  • Conservation of angular momentum through timewarp and outside timewarp (allowing for spin-stabilisation and thus fixed RT2 antennae if Cilph decides to implement that);
  • Thrust through timewarp, for ion engines and the like, together with maneuver nodes for non-instantaneous burns;
  • Orbital decay drag, with ferram4's help;
  • Stationkeeping, to deal with orbital decay and unstable weakly bound orbits;
  • Gravitational moment (quadrupole) for planets, enabling such amusing things as Sun-synchronous orbits;
  • Insert crazy idea here...

 

Acknowledgements

I would like to thank (in alphabetical order of forum usernames)

  • Anatid for his attempt at documenting the KSP API,
  • armed_troop for his help in making the codebase clang-compatible.
  • Ezriilc, Forecaster, khyperia, Payo and sirkut for their HyperEdit plugin, which is really useful for setting up initial conditions,
  • Scott Manley (illectro) for his help with numerical integration of the N-body problem,
  • Matt Roesle (Mattasmack) for writing an integrator which was far better than the one I used, thus prompting me to write my current SPRK integrator,
  • NovaSilisko for his Alternis Kerbol mod, which piqued my interest in the simulation of the N-body problem,
  • The entire KSP modding community (especially the subset which writes clean, well-commented code) for providing code examples from which one can learn the subtleties of dealing with the KSP API.

 

If I forgot you in the above list, please complain! :)

Bibliography

 

 

Edited by eggrobin
Cantor is out.
46 people like this

Share this post


Link to post
Share on other sites

Very impressive and ambitious -- will be a very different game if orbits are precessing and don't stay the same while you're not looking!!

Share this post


Link to post
Share on other sites

This is very cool, I'm looking forward to seeing where this goes. I suppose I should try and answer this question:

Would it be conceivable to predict the results of aerobraking/aerocapture/ballistic reentry with FAR (within some error margins, to be displayed on the predicted trajectory)? How slow would that be?

It should be able to predict, but there are serious caveats to this. If the vehicle enters in an stable aerodynamic configuration the results should be predicted fairly accurately, though there might be some errors due to initial transient effects on atmospheric entry. However, if the vehicle is unstable, the results will be completely inaccurate; unstable aerodynamic configurations tend to follow highly chaotic motion, which makes them essentially unpredictable due to how much minor errors can affect the trajectory. This means that there needs to be a way to determine whether the predicted trajectory is that of a stable or unstable vehicle; I suspect the best way is to search for sudden changes in the orientation of the vehicle and color-code everything after that to designate the degree of chaos in the trajectory. If that also accounts for stability induced by creating a large amount of spin-stabilization, it would provide a fairly accurate way to estimate the effects of aerobraking.

As for the time to do this... it will likely be slow, since you'd have to simulate it part-by-part for best results. A faster (but less accurate) solution would be to calculate out the vehicle's properties as best as you can at the hypersonic limit and use that, which works for Mach numbers > ~5 or so. This will work for most aerobraking in stock KSP and will certainly work for RSS.

Also, as a guess at an answer for this one:

The method I currently use for drawing trajectories (LineRenderers, object in layer 9, transform.parent = null, useWorldSpace=true, updated every time the GUI is drawn) has some issues: when rotating the camera in map mode, the lines lag behind. Does this have a known solution? (I guess so, RT2 doesn't have this problem).

Are you setting it in Update or FixedUpdate? Make sure you use the former, I've seen effects lagging behind when using FixedUpdate. If that's not it, try seeing if you can do things in localSpace.

Share this post


Link to post
Share on other sites

Thank you for creating an extremely well-put-together post. It would be so easy to just say "I'm doing N-body stuff, here is a screenshot," and then the thread would immediately fall to pieces.

Also this is totally mind-blowing excellence.

Share this post


Link to post
Share on other sites

Well I'll be damned. I didn't think it was possible! Best of luck!

Share this post


Link to post
Share on other sites

Everyone: thanks!

This is very cool, I'm looking forward to seeing where this goes.

Thanks!

Are you setting it in Update or FixedUpdate? Make sure you use the former, I've seen effects lagging behind when using FixedUpdate. If that's not it, try seeing if you can do things in localSpace.

I'm not sure when it happens (I have no idea how Unity actually works). I'm registering a callback to draw the GUI, and I draw the trajectories at the same time. I'll try moving this to Update(). It seems this is due to the origin of ScaledSpace moving (it's always near the camera, so that LocalToScaledSpace(0) != 0) and the trajectories staying in the same position relative to the ScaledSpace origin, so switching to localSpace (whatever that is) might fix it.

Thank you for creating an extremely well-put-together post. It would be so easy to just say "I'm doing N-body stuff, here is a screenshot," and then the thread would immediately fall to pieces.

You're welcome!

I figured if I make a dev thread, I might as well make one from which I can gather useful information for development. This seems to be working.

Share this post


Link to post
Share on other sites

Well it certainly seems like you've got a grasp on the problem. I don't think I've ever seen a post with a bibliography before...

Very interesting work, and I hope it goes well for you. Best of luck!

Share this post


Link to post
Share on other sites

The question no-one here is asking:

Will it work with Kraig's Planetfactory?

Share this post


Link to post
Share on other sites
Thank you for creating an extremely well-put-together post. It would be so easy to just say "I'm doing N-body stuff, here is a screenshot," and then the thread would immediately fall to pieces.

Yeah I also appreciate the level of detail both of things to expect and what not to expect from this project. Helps keep the endless questions at bay. Not to say there won't be gone completely though. XD

Share this post


Link to post
Share on other sites

It's like FAR, but it affects gravity this time! Woo!

Share this post


Link to post
Share on other sites

This sounds mind-bogglingly awesome; together with FAR, this will put KSP on a par with Orbiter in terms of realism.

I can only second Headhunter09's post; I have a sneaky suspicion that you're deep in academia, and this might become your master's or even PhD thesis.

Share this post


Link to post
Share on other sites

What would this do to long term satelites?

Would we have to control each one every so often to do station keeping burns?

Share this post


Link to post
Share on other sites
What would this do to long term satelites?

Would we have to control each one every so often to do station keeping burns?

I think this woudl require some kind of automatic station keeping mechanism to keep things playable... it probably woudl be now fun to periodically manually visit all these 20 satellites on orbit to adjust their trajectory.

And one thing came to my mind. If we are going full realism, how about peturbations due to planets not being perfect sphericam bodies?

Share this post


Link to post
Share on other sites
The question no-one here is asking:

Will it work with Kraig's Planetfactory?

Maybe. It should, but then I haven't looked too closely at the implementation details and I have heard reports of incompatibilities between PFCE systems and FAR. It will be some time before things are playable on my side though, so I really can't say what PlanetFactory will look like by then.

I can only second Headhunter09's post; I have a sneaky suspicion that you're deep in academia,

What gave it away? :P

Undergraduate student in Mathematics at the ETHZ.

and this might become your master's or even PhD thesis.

I don't think this really counts as pure math. :)

What would this do to long term satelites?

Would we have to control each one every so often to do station keeping burns?

I think this woudl require some kind of automatic station keeping mechanism to keep things playable... it probably woudl be now fun to periodically manually visit all these 20 satellites on orbit to adjust their trajectory.

Yes, I forgot to add that to the 'further modding' list. It turns out however that orbits far away from the critical points of the gravitational potential (Lagrange points in the CR3BP) are pretty stable on their own. Predictions should help for the more complicated trajectories.

Once orbital decay drag is modeled, automated stationkeeping will definitely be needed though.

And one thing came to my mind. If we are going full realism, how about peturbations due to planets not being perfect sphericam bodies?

Already in the 'further modding' list:

Gravitational moment (quadrupole) for planets, enabling such amusing things as Sun-synchronous orbits;

See 'Gravitational Moment' in Eric Weisstein's World of Physics (I will add this to the bibliography).

Edited by eggrobin

Share this post


Link to post
Share on other sites

Although there is just only proof of concept available this is absoluteley mindblowing . The structure of the OP is very detailed and shows you are really deep in working on that.I will stay tuned, can't wait when first stuff will be available :D

Share this post


Link to post
Share on other sites

Awesome! I have been privately working on something similar. May scrap it now :( How are you handling the fact that stock trajectories of celestial bodies do not follow full N-Body physics? Most prominently, Kerbin is fixed in the Mun-Kerbin system, where really, both Mun and Kerbin should orbit around the common COM. That makes the Mun disturb objects in LKO much more than it should; there should be just tidal forces, but with a fixed Kerbin, the full gravity of the Mun perturbs the orbit, not just the differential to Kerbin's COM. Are you setting the massive bodies free?

Share this post


Link to post
Share on other sites

If you do this I'll have to learn the math for weak stability boundary transfers.

Share this post


Link to post
Share on other sites
Speed will be further increased by rewriting the numerical integration in (unmanaged) C++, interfaced with C# through C++/CLI. The unmanaged C++ (compiled with clang) will use 80-bit extended precision 'long double' (64-bit mantissa instead of 53-bit) for added precision.

If speed is the aim, SSE extensions could be a faster alternative to x87 floating point but would require sacrificing extended precision.

Share this post


Link to post
Share on other sites

HOLY S***! YES!!!! I hope this could be implemented at some point in the future. As far as computation speed goes, what are you running this on and is it at all playable?

Edited by Jobin

Share this post


Link to post
Share on other sites

Very good development, I've attempted something like that myself but run out of spare time so it was never finished. I just want to comment on one thing.

Speed will be further increased by rewriting the numerical integration in (unmanaged) C++, interfaced with C# through C++/CLI. The unmanaged C++ (compiled with clang) will use 80-bit extended precision 'long double' (64-bit mantissa instead of 53-bit) for added precision.

My experiments have shown that using long double is non-starter in terms of performance. I used doubles and written a host of vector madd (multiply+add) functions using SSE2/3 in C++ and C++/CLI as "glue" interface. You can theoretically make them even faster if you utilize AVX commands, but that would severely limit users who would actualy be able to run it as these commands are relatively new, while SSE2/3-capable CPUs are widespread. Since I used RK5(4) integration which is not much except a ton of madd's and got very good performance out of it. However there are few things to consider:

1. mixed-mode C++/CLI does not conform to CLR specification and as such it is not given that it would work under Unity - I haven't actually tried that and was focusing on optimizing algorithm performance.

2. Unity is cross-platform and runs on MacOS and Linux. Again I'm not sure that mixed-mode would actually work on these platforms.

3. To make matters even worse, there is 64-bit version of KSP for Linux so you would have to ensure you provide 64bit binaries for that platform.

Edited by asmi

Share this post


Link to post
Share on other sites

You have my support and interest in this N-body endeavor.

Edited by Cydonian Monk

Share this post


Link to post
Share on other sites
Are you setting the massive bodies free?

Yes. For fine timesteps, this actually is easier than propagating them on their Keplerian rails, because solving the Kepler equation (which involves a few iterations of Newton's method; the derivative of Kepler's function has a trigonometric function, which is a few hundred clock cycles compared with 10 or so for a division/square root) is pretty slow.

It should be noted that even if you want to keep things in 2-body orbits, it is a good idea to really use the solution of the 2-body problem (2 bodies orbiting a barycenter) rather than the pretty poor Keplerian approximation (one body orbiting another).

Share this post


Link to post
Share on other sites
Very good development, I've attempted something like that myself but run out of spare time so it was never finished. I just want to comment on one thing.

My experiments have shown that using long double is non-starter in terms of performance. I used doubles and written a host of vector madd (multiply+add) functions using SSE2/3 in C++ and C++/CLI as "glue" interface. You can theoretically make them even faster if you utilize AVX commands, but that would severely limit users who would actualy be able to run it as these commands are relatively new, while SSE2/3-capable CPUs are widespread. Since I used RK5(4) integration which is not much except a ton of madd's and got very good performance out of it. However there are few things to consider:

1. mixed-mode C++/CLI does not conform to CLR specification and as such it is not given that it would work under Unity - I haven't actually tried that and was focusing on optimizing algorithm performance.

2. Unity is cross-platform and runs on MacOS and Linux. Again I'm not sure that mixed-mode would actually work on these platforms.

3. To make matters even worse, there is 64-bit version of KSP for Linux so you would have to ensure you provide 64bit binaries for that platform.

The original mapping mod was IIRC using mixed-mode C++/CLI successfully, so I doubt that'd be a barrier, and whilst the KSP game itself is on multiple platforms, there's nothing to say the OP *has* to provide working ports to MacOS or Linux (or indeed his primary could be Linux/MacOS, and it could be win32 he doesn't have to support), unless I've misunderstood the mod developer rules?

Edited by timmers_uk
unbalanced braces

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!


Register a new account

Sign in

Already have an account? Sign in here.


Sign In Now