eggrobin

Let me ask a very silly question. How slow would CFD on the meshes have to be in order to be good enough for our purposes? I have no experience in aerodynamics (and generally next to no experience in solving PDEs numerically) so I don't have a clue what the order of magnitude is. I suspect it's way too slow, otherwise you'd have done it. EDIT: some idle Googling seems to show that this is not quite outside the realm of feasibility with lattice Boltzmann methods using GPGPU (the Wikipedia page says they don't work well for high Mach numbers though). Apologies for treating your field like http://xkcd.com/793/, I'm afraid this also applies to mathematicians.

Well, it is easy to check that the floating origin isn't with the ship: just trace the ship's worldpos. One easily sees that the origin tends to be a few km away until you timewarp (it is properly reset by timewarp). I'm not having any issues with Krakensbane at the moment. It seems the planets aren't quite where I put them though and this (or something else) is causing problems when trying to set the vessel's position when off-rails. This is driving me ever so slightly insane.

Good points. There's also setOffset or whatever it's called, which would enable me to manually set the origin to the ship's position as often as possible (since all the physics takes place within a sphere, the only sensible location for the origin is at the center of this sphere). Anyway, my main concern at the moment is setting the position correctly when off-rails. What I'm doing at the moment differs significantly from Unity's calculations (i'm diverging from Unity by 1 m after 10 s, so something's wrong).

If that is the case (and experiments show it could be the case) it's pretty silly. They cut physics off at 2.5 km from the ship due in part to floating point inaccuracies, but they allow the ship to be 6 km away from the origin. Frankly, there is no good reason for the ship not to always be at the origin.

Well, I don't really feel like interfacing with FORTRAN, I want to add some effects beyond point-mass gravitation (thrust, drag etc.) and writing the integrators is the fun part. It seems it is the easy part too (it's just maths, physics and numerics; the rest involves writing untestable code to interface with an undocumented, buggy and ill-designed API). Speaking of which, Status update: 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.

Welcome to the fora! Caveat lector: I am not an astrophysicist nor a numerical analyst, I am an undergraduate student in mathematics whose only experience in numerics stems from taking a couple of courses on numerical mathematics, so the statements below reflect my incompetence. I fail to understand how unit conversion would speed up the calculation. They are currently done in an inertial reference frame whose origin lies near the sun. The units are those used by KSP (m for q, m/s for v). Should I switch to furlongs and furlongs per fortnight? The reference frame being nonrotating, it is necessary to rotate the result of the integration to use it in KSP, whose reference frame sometimes rotates. My experiments on simulating the Alternis Kerbol system with a satellite in LKO seemed to show that it performed well as far as errors in H, p and L are concerned (the results were similar for reasonably larger timesteps; fifth order convergence was observed). In the current prototype, the timestep is such that the error should be no more a few unit roundoffs anyway. I have not implemented the smarter integrators yet, but this will be the next step after fixing a few bugs with the moving around of off-rails vessels and a bit of refactoring (the last thing I want is my codebase turning into a mess of spaghetti code before it is even playable). I'll start by implementing the Saha & Tremaine stuff (individual timesteps and separation of the Hamiltonian in H = (T + VSun) + Vinteraction rather than H = T + V), generalising it to higher orders, and then I'll see whether it is good enough at the timesteps I will use. If the Saha & Tremaine integrator is not good enough (I suspect this will only be a problem for the dynamic integration, the main integration is almost solved by the current SPRK integrator already, so I'll be able to use comparatively small timesteps with Saha & Tremaine) I'll look into the Chambers paper. As for the choice close-encounter integrator, I think I'll go with the usual method of implementing a few integrators and seeing how they behave in-game.

Status update: 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. The direction of a does not seem to bear any relation to ̉ۡ. EDIT: It seems a is either prograde or retrograde. The number of astrophysicists here seems to have doubled recently.

What's a PMI? EDIT: Oh, it's just the off-diagonal entries of the inertia symmetric bilinear form. Engineers have weird naming conventions when it comes to mathematics. My original intention for off-rails was to derive the proper acceleration of the vessel from first principles by talking to the engines as you say, but as discussed earlier, there are lots of sources of proper acceleration: besides engines, you have decouplers, collisions/getting out and pushing, engine exhaust, explosions etc. So computing the acceleration from first principles entails knowing what can exert thrust. This would break compatibility, and would lead me to effectively writing my own physics engine. Since this is definitely beyond the scope of this mod, I have to grab the acceleration from Unity. As Unity is bad at physics (and the KSP layer is pretty lousy too), I will have to postprocess things.

No, it means I'm working on the problem of making this playable, i.e., of making rockets useful. Currently, off-rails thrust is ignored. This means you can't change your orbit. I'm trying to grab whatever acceleration comes from your engines/getting out and pushing/decoupling etc. and to change your trajectory accordingly. Remember that in this mod, I'm the one moving stuff around, not Unity (this also applies to off-rails), so I have to fetch these accelerations in order for them to affect the trajectory. On-rails (timewarp) thrust is a feature I will work on once the basic stuff works, like the rest of the 'further modding' section.

You make some very good points. I guess if I need more significant digits I'll use double-double rather than extended precision then. No, I'm compensating for the geometric accelerations applied by KSP. What I mean is what you see when you spin the ship very fast: the camera (which is centered on the center of gravity) wobbles because Unity doesn't conserve momentum in its interactions between rigid bodies. It seems most of the biases I was measuring when the ship was not rotating very quickly were coming from the reaction wheels though; I'm not sure if that's because I have a joystick or whether it's a deeper bug, but with reaction wheels enabled, a test ship (Mk1 - FL-T200 - LV-T30) picks up angular momentum quickly without input (SAS off), while this phenomenon is much slower with reaction wheels disabled. Ignoring the reaction wheels issue (which may well occur between the chair and the keyboard), it seems a low-pass filter for rounding errors, with a band-stop around the rotation rate and a maybe a high-pass on torque might do the trick. I'll see about what post-processing I want to do later on, this will require quite a bit of logging and data analysis. FFT is hard; let's go integrating.

Status update: 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 ). 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.

This pretty much has to exist since thrust through timewarp is planned. I guess I'll need to have a few options as far as orientation is concerned (thrust through timewarp has to be automated, since you can't control the ship directly during timewarp).

By adding them to the 'further modding' list. There are quite a few things I will want to show on the predicted trajectory using interactive markers (you should be able to click on them to get more details). Examples include aerodynamic events (free molecular flow to fluid flow transition, Mach numbers, instability), stationkeeping burns, long-duration burns (start and end). A maneuver node can be such a marking. When tweaking it, the new trajectory is computed using the dynamic integration (I really don't like this name, can someone come up with something nicer?). As for the interface for tweaking it, I think I'll go with a window and sliders; the current system has ergonomic issues. An interesting question is that of the axes. In stock KSP, we are always in a Keplerian orbit, so it makes sense to see burns as a linear combination of prograde, radial and normal. This will also work for strongly bound orbits. For weakly bound orbits, other axes may be useful. Remember that double precision limits the size of the system. The resolution with a 53-bit mantissa at the aphelion of Pluto is 0.8 mm, so you can't do much bigger. If I were to use extended precision (which is slower, as x87 is a scalar FPU, but I don't know how critical performance will end up being), the resolution at 1 light-year would be 0.5 mm. This would be much nicer for RSS. If performance really isn't an issue, I could try double-double (implemented in software, so an order of magnitude slower), which would give me an error of 0.01 mm at the edge of the observable universe (though for systems that big you'd really want GR, and I don't want to look at the curvature of that rabbit hole). I could talk about IEEE 754 quads, whose error at the edge of the observable universe would be 87 nm, but this is getting rather silly unless you're interested in molecular dynamics at the scale of superclusters. Anyway, I'm not currently in the planet-making business, and I hear PFCE is rather buggy (e.g. compatibility with FAR or lack thereof). I'm not sure if the problem of making a fully functional star, i.e., a star which is a light source, feeds the solar panels, etc. has been solved either.

From the game's point of view, the SOIs are still here and the biomes should still work like they did. It's just that I'm moving stuff around instead of letting the game do its thing, and I'm ignoring SOIs. It's actually really messy to define properly. There are mostly three definitions for the SOI of a body---for concreteness, let's call it Earth (in a system with two massive bodies: the Earth and the Sun), The sphere of attraction is the region in which the force exerted by the Earth is greater than that exerted by the Sun. It is tiny and mostly useless (the Moon is outside of Earth's sphere of attraction with respect to the Sun). The sphere of action, or Laplace sphere (Laplace 1799) is a tad more complicated to define, and finding the definition on the internet takes a while: look at the unperturbed problem of a (massless) test particle orbiting the Earth, and consider the position rEarth of the test particle relative to the center of the Earth. The particle undergoes some acceleration r"Earth = aEarth. Now add the Sun. Call pEarth the resulting change in the second derivative of rEarth (it's the perturbation due to the Sun of the 2-body approximation of the relative motion of the test particle and the Earth). Similarly, define aSun and pSun. The sphere of influence of the Earth is the region where the acceleration to perturbation ratio for the Earth is smaller than the corresponding ratio for the Sun, namely, aEarth/pEarth < aSun/pSun. In this region, the 2-body approximation around B1 is better than the 2-body approximation around B2. The sphere of action is larger than the sphere of attraction. This is the SOI used by KSP. The proper definition for the Hill sphere is slightly more involved. The idea is that when your orbital energy (in the corotating frame of reference for the CR3BP) is low enough (and you're close to the Earth), there is a bounded region around the Earth that you can't leave. If you keep increasing your orbital energy, there's a point where this region will include the Lagrange point L1 and will suddenly expand so that you are able to reach the sun: when you have this orbital energy, conservation of energy does no longer prevent you from escaping the Earth. Just before this happens, the region around the Earth that you can't leave (which extends up to L1) is the Hill sphere; it's the region inside the equipotential through L1 (on the Earth side). If you want to visualise this, look at the image below (the rotational potential of the CR3BP) and imagine having a horizontal surface of water that rises on this weird hill with two holes. Say you are on a boat in the tiny lake around Earth. The area where there is water is the area where your orbital energy allows you to go (by trading kinetic energy for potential energy). You could be in the lake around the sun, but there is a mountain pass at L1 between the two lakes preventing you from going from one to the other. As the water (the orbital energy) rises, there is a moment where the lowest point of that mountain pass (L1) is submerged. From then on there is nothing preventing you from going to the Sun. Just before the lakes merge, the lake around Earth is its Hill sphere. Page 580 of this paper by Ross & Scheeres 2007 shows the shape of the forbidden regions depending on the energy (CJ = -2 H, where H is the energy, because astronomers are weird and count the energy backwards). The Hill sphere is larger than the Laplace sphere. It should be noted that Jupiter has quite a few satellites outside of its Hill sphere, so having an energy below that of L1 is not a required condition for having a stable orbit. The Laplace sphere is mostly useful for determining which patched conic approximation to use (this is why it was defined together with the patched conic approximation by Laplace in his Traité de mécanique céleste). Indeed, none of these are actually spheres (the first two get more spherical as the ratio of the small mass to the larger goes to 0, but the Hill sphere remains pinched at L1).

Is this also the case if you run RSS / Realism Overhaul? The mass distribution should be better with these mods.

There would be quite a bit of stationkeeping to do (maintaining a couple dozen LKO satellites in LKO with orbital decay, keeping stuff around Lagrange points etc.), so this would get quite tedious. Moreover, stationkeeping has to be done quite regularly. If you had to do it yourself (or, even if that were automated by kOS, if kOS did the stationkeeping by loading the vessels), you couldn't timewarp for very long, rendering interplanetary missions impractical. In addition, NeverUnload is not a solution: the floating origin stays where it is and Unity physics is done in single precision, so if the active vessel is on Duna, a vessel performing stationkeeping in LKO would jump 2km at a time. Needless to say, instead of performing a minor correction to its orbit, it would end up in a random orbit (or crash on Kerbin). Finally, stationkeeping is way harder than you seem to think, see thorfinn's post. You can't conceivably write the gradient descent you'd need to minimise this cost function in kOS, not to mention an integrator. You'd need a full-fledged programming language and a few weeks of careful coding (not counting the time to do research, and assuming you have experience in numerical math) in order to write your stationkeeping algorithm. You'd basically be writing a 1000-line plugin that does stationkeeping. While this is fun, and I recommend you read the papers and my eventual implementation, I'm not sure this is a good idea from a gameplay perspective. As a conclusion: you need stationkeeping, you need it to be automated, you can't load the vessel to do stationkeeping and kOS wouldn't do as programming language.

I'm downsampling by a factor of 100 to plot the trajectories in order to improve performance and reduce memory usage (I probably don't need to for such short predictions, but I didn't really think about it when I typed 100; I need that for predictions over a decade). As a result, fast orbits look polygonal (you get 1 point every ks). Eventually, smarter downsampling (velocity-dependent) should be implemented, as well as something more advanced than LineRenderers. In the meantime, I have pushed a new version that should somewhat alleviate this (by reducing the downsampling to 100 seconds from 1000 seconds). The barycentric coordinates are also somewhat broken (you'll occasionally see obvious jumps in the trajectory). I suspect this comes from Unity's Quaternion. I am using Quaternion rather than QuaternionD because the obsolete radian-based methods are not implemented for QuaternionD; I don't trust either---why should I trust something that thinks constantly converting between radians and degrees is a good thing---, I'll just write my own versor.

Good point (also, I had never heard of cyclers, thanks for mentioning that!). I guess I'll end up having a few hardcoded options. Regarding people writing their own stationkeeping, I guess something like kOS, with a proper language and extensive APIs, would be nice... *daydream mode on* contemplates making a plugin that interfaces with Ada, provides some packages for vehicle control, and integrates GNAT with a minimalistic IDE in-game. *daydream mode off* Let's get back to the current crazy idea.

My question is mostly a gameplay one actually. The most sensible reference frame is the rotating-pulsating one, but I'm afraid this might be a tad confusing on the map (it's not an isometry). It is also convenient to use true anomaly, rather than time, as the independent variable, which affects what the potential looks like, but I think this is less confusing than the lack of isometry induced by a pulsating map view. The thing is, you don't deviate from your predicted orbit (the prediction and the actual orbit are computed in the same way). For LKO and the like, you could do what you say with the deviation from the Keplerian orbit, thus compensating for orbital decay. The question is, how do you say you want to stay in a Lissajoux or halo orbit around a Lagrange point? Or a horseshoe orbit? Or another kind of exotic orbit? You can see from the prediction that the orbit will end up being unstable, but saying "I want to make corrections to keep me in that kind of an orbit" is something complicated (this is why astrodynamicists exist ). I don't need no stinkin' physics engine. Actually, I'm trying to prevent Unity from doing too much anyway, so I'll definitely not be loading the vessel. The idea would be to model the vessel as a rigid body (no, not an Unity RigidBody, an actual rigid body), see where the thrusters are and what torques/forces they apply, compute the resources required for the correction, perform the correction by just feeding it to the integrator and subtract the resources from the vessel.

That would be interesting. "Hello, I'm Scott Manley. Today we'll to throw a main-belt asteroid at Eeloo by crashing a Kerbal on EVA into it. Let us first derive a closed-form expression for our trajectory." I don't think asteroids are in 0.24 though, they are in the NASA thing. Interesting conclusions from the poll results: It seems Scott's not here. Could somebody draw his Manley attenion to the remaining open astrophysics question (the one about reference frames for the ER3BP)? Answering it would do too of course. The majority seems to know basic calculus, so I'll use that. I'll post a link to Feynman's Lectures on Physics for the ~40% who don't, there's a quick introduction to calculus which doesn't involve defining the derivative at a point as a linear form on a Banach space. It's really well written (I'm pretty sure that's where I first learned calculus around 2007). There are 6 people who checked neither of the bottom 2 options. Either they can do numerical integration but not elementary physics, or they didn't see that it was a multiple choice poll I think the first post I'll write will be on the integration of general ODEs with RK integrators. In later posts I'll talk about Hamiltonian mechanics, the Kepler problem and the patched conic approximation (it would be nice if people understood why SOIs make sense in their current form, but not as Lagrange point hacks or any such sort of nonsense). I'll then get to SPRK integrators, and eventually to the Saha & Tremaine stuff. I'll try to write the first post over the evenings of the coming week. There will be little modding this week because exams. I'll probably be more productive from Saturday onwards though.

As previously stated (though not in the OP), the unstability of the Jool system can probably be fixed by tweaking the initial conditions (actually, interpreting them differently). This would remove potential collisions. Handling collisions of celestial bodies is currently not on the roadmap.

When I heard about it in Rock, Paper, Shotgun. I got the first version with symmetry. It's interesting that it seems to peak between .17 and .20. The comparatively small number of people from the beginning probably stems from the fact that the game was not very well-known at the time (assuming we have a representative sample).

I can't, it's been a really warm and dry winter here in Switzerland. Though it's currently sort-of-snowing-but-sunny-and-windy-and-it's-+7-degrees-centigrade-so-how-is-it-snowing-anyway. They are. The first-order approximation they have makes sense from a physical standpoint and saves them the expense learning how to do numerical integration: it's going to be several months before I get to the stage where things are playable and where I have been rewritten the integration in C++ for better performance, even though I have taken courses in numerical maths, received advice in numerical integration from Scott Manley and constant help from my fellow modders. It turns out time warp isn't as complicated as you'd think, astrophysicists have already solved this problem for us. See the OP and my upcoming explanatory posts. I do not doubt that we can do timewarp as fast as we need it and with good accuracy. Thoughts: It's Joseph-Louis Lagrange or Giuseppe Lodovico Lagrangia, not LaGrange or La Grange or Geordi La Forge. Yes, "la grange" means "the barn" (or "the grange"), but his name is written in one word without CamelCasing. Similarly, mine is Leroy, not LeRoy, Le Roy, Leroi, LeRoi nor Le Roi The only reason the SOI system works well (ok, it doesn't, the planets should really move around barycenters rather that each other, but let's assume they do) is that it's a first-order approximation of the equations of motion around the planet whose SOI you're in. In order to have this kind of approximation near Lagrange points, you'd have to use the 2-body Hamiltonian anyway, so you'd have to use an integrator (the solutions to the 3-body problem are given by a slowly converging series useless for practical purposes). If you have to do numerical integration, you might as well drop the first-order approximation and do N-body physics. I'll explain this point in further detail when I get to my series of exlpanatory posts. From a gameplay standpoint, SOIs might make sense, but regions around the critical points of the potential where you suddenly choose to change your approximation don't: How big would these Lagrange SOIs be? Around Lagrange points you would be integrating the 3-body problem, but when you get somewhat closer to a planet you'd use the 2-body approximation. I'm pretty sure your Lagrange SOIs shouldn't be spheres. Post-Scriptum: This is also a very good point. Lagrange points are useful landmarks when looking at trajectories in the CR3BP, but trajectories differ from the patched conics approximations in many other places (basically, anywhere except low orbits around either massive body, or in really high orbits around the system). If you had 'Lagrange SOIs' you wouldn't get horseshoe orbits. I should watch this movie someday. I don't know how I will do it yet (the stationkeeping, not watching the movie). The problem is defining which orbit you want to stay in. I'll add this to the list of open problems, there are people on these fora who have experience with that sort of stuff (Arrowstar comes to mind, there are probably many more).I'll get to work on that low-pass filter! (actually, I think I'll start by logging the data and doing DFTs in Mathematica to see what things look like). EDIT: come to think of it, I shouldn't get bogged down in data analysis. I'll try to do the integration without any particular preprocessing and I'll worry about filtering later.

Brilliant! I hadn't thought of looking at it in frequency space (and I missed it in your previous post because it's 4:40; I'll go get some sleep). Trying to grab the most noise-free acceleration is probably a good idea though. I think by computing the change in momentum myself (perhaps sorting the momenta to minimize cancellations) I can avoid the double->float->double conversions that KSP is likely to do (I mean, FlightIntegrator exports the Coriolis force in single precision, but FlightGlobals computes it in double-precison ). You have to be careful though, RCS can fire frequently. No, but it would prevent the "stick an ion engine on the side" technique. Automated stationkeeping (and attitude control) is on the roadmap, so that would just end up using up your ACS fuel very quickly. Real programmers don't use Pascal. This article has its faults, but it does make a good point that in Pascal, there is no escape, whereas you always have unchecked options in Ada ('Unchecked_Access or Ada.Unchecked_Conversion being examples). There are cases where you have to use them, and you're glad they exist. In the context of trying to escape Unity's incompetence, this is amusing... EDIT: RCS thrust probably doesn't flip its sign every other FixedUpdate, so I think I'm safe.

That's the second language I learned (after VB6 ) and the first one in which I wrote algorithms (in Ada 95, if you want to sort something, you write a sorting algorithm. Same goes for most useful data structures, so you can't really avoid it). The fact that my father's job at the time included defining the standard for Ada 2005 might have had something to do with it. No, I mean "under what threshold should I clamp the accelerations and attribute them to Unity's ineptitude"; long-duration burns will be handled by talking to the engines (I can't use Unity physics to compute the acceleration during timewarp and if I could I wouldn't want to). EDIT: It's probably not too hard to integrate the torque due to the engines in timewarp. It turns out that the hard part is doing physics when Unity is trying to do its thing, rather than when we don't use the physics engine. EDIT2: A well-known use of Ada in aerospace would be the Ariane 5 software, which was reused from Ariane 4. Ariane 5 was faster, this led to a Constraint_Error during a type conversion because they had removed a range check due to efficiency concerns. As a result, the rocket blew up. The spacecraft was Cluster. Turns out if you do stupid stuff, even a nice language with strong typing and range checks everywhere won't help you (also, their identifiers have terrible names).