Open se5a opened 5 years ago
Retrograde orbits being broken is a side affect of trying to flatten a potentially 3d orbit into drawing a 2d image. as such if i > 90 degrees, (indicating a retrograde orbit) the AoP (and TrueAnomaly?) need to be sign flipped.
To help debug kepler orbits I've created a widget to display different lines/distances, positions, and angles. This can be found by:
@UberWaffe Are you any good at this type of maths? a fresh pair of eyes on this would be nice.
@se5a I'll give it a look. But it has been years since if done any serious physics math.
Still trying to get the UI to display properly on Windows. So most likely I'll be stuck with testing it via unit tests.
Cheers, thanks. There are some reasonably comprehensive tests that are failing, but if you manage to get the ui displaying it may help you visualize what's going on, especially with the retrograde orbit stuff.
I've got an old netbook with windows on it, I'll see if I can get that running to have a look at the ui on that, unfortunately it doesn't seem to like running in a virtual machine. Unless there's missing libraries it should "just work"... if you can create an issue for it with a screenshot it would be good.
Do you have a URL link that explains the maths better? I've done some searching but I don't quite get why we are trying to get vector perpendicular to the angular velocity vector and the Z-axis. (nodeVector) https://github.com/Pulsar4xDevs/Pulsar4x/blob/91ce915bd6432fe1ecd7aa072a61d43817e6e5f1/Pulsar4X/Pulsar4X.ECSLib/Helpers/OrbitMath.cs#L49 If I am visualizing this correctly, the angular velocity vector is already representing an instantaneous plane perpendicular to the orbit (and therefore in line with the Z-axis). So the cross product of that and a vector that is only in the Z-axis (the line 49 references) always gives a 0.
Here is the visualisation (as I understand it) from Test 1 (the one that passes, but that appears to be a fluke if I step through the code.)
https://academo.org/demos/3d-vector-plotter/
I know the scale is off, since I didn't convert the Km to AU, but the cross product plane remains the same. The result is in the Z-Axis.
The cross result of this and the hard-coded vector in the Z-axis (0, 0, 1) is not defined. (It is just all zero.) This results in all the division calculation further down that uses nodeVector to result in NaN. (Divide by zero error.)
I'm trying to find the main source of info I used for that, but these are places I've also referenced:
http://orbitsimulator.com/formulas/OrbitalElements.html but more specificaly the source: view-source:http://orbitsimulator.com/formulas/OrbitalElements.html
it's possible that the accepted answer here was my main reference, but I thought I had something else that was a bit clearer than this: https://space.stackexchange.com/questions/1904/how-to-programmatically-calculate-orbital-elements-using-position-velocity-vecto it looks like it might have a similar formula there though, so should give some idea of what it's trying to do.
Also this for some calcs: http://www.bogan.ca/orbits/kepler/orbteqtn.html
also pretty much every wikipedia page on kepler elements. https://en.wikipedia.org/wiki/Longitude_of_the_ascending_node it's a bit late tonight for me to dig further and understand the above, I'll take a look in the morning. about all I can glean right now is that it's meant to be used to find the Longitude of Ascending Node.
Edit: it's possible that this was my prime reference: https://github.com/RazerM/orbital/blob/0.7.0/orbital/utilities.py#L284
@se5a Thanks, I'll dig into the references and update my knowledge.
Still busy investigating, but found a few things. (I've added some tests specifically for the Kepler Elements calculations. The experimental branch can be found here. https://github.com/UberWaffe/Pulsar4x/commit/9b235033b3e13abe28ea84b2f67bbc4b3a2f6e3f) (Be careful of merging the branch. I also played around with updated ImGUI.Net to a newer version to see if that resolves the Windows UI issue. It doesn't. The last two commits can be cherry picked rather safely.)
Firstly, there are rounding errors. If I do the tests, if I demand accuracy beyond one decimal place, then the tests fail. I would suggest we do the calculation in km and use decimal for the data type. It has nearly twice as many digits precision, and when working in m (not even km) can cover distances from -1x10^12 light years to +1x10^12 light years. That should be sufficient even for the multiplications of vectors.
Secondly, for the second test (the first that fails), if I ignore the accuracy issue, then the following elements are out / incorrect: LongdOfAN => Expected = 45.000001531994371 ; Calculated = 0 ArgumentOfPeriapsis => Expected = 3.1407949771775465 ; Calculated = 3.9269908169864447 MeanAnomalyAtEpoch => Expected = 3.2214551857593583 ; Calculated = -3.1415901287309396 TrueAnomalyAtEpoch => Expected = 3.142390330156426 ; Calculated = -3.1415926535889964 Periapsis => Expected = 0.00014095441681280844 ; Calculated = 1.4089617201124753E-13
I cannot tell you why yet (I don't get the maths yet. Rocket science, go figure). But I figured mentioning the results so far might help others in tracking the problem down.
I've been/still am fighting with trying to update virtualbox all morning, so I've not had a chance to look at stuff closely yet, where do you think the rounding error is happening? I did notice that but thought that a single decimal place was too much to be a rounding error. (if we're talking radians that's a massive difference in angle) Try doing the calc in decimal and seeing what happens - the reason I'm passing the sgp into a lot of the orbitmath calculations was so that you could pass them the distances in whatever units of measure you choose. tbh I'm not a huge fan of storing distances in AU, but the effort of changing it and deciding what to store it as wasn't worth it. Not too keen on storing and doing everything in decimals, not keen on the performance loss (that being said, I've not profiled it) anything thats not being done all the time could be converted and calculated in decimals if necessary though (ie conversions between state vectors to kepler elements does not happen often enough to worry about performance, kepler to state vectors does though) I do want to go through at some point and change all the angles to be stored as radians though.
are those expected/calculated values in radians or converted to degrees?
"are those expected/calculated values in radians or converted to degrees?" I had missed converting LongdOfAN to radians before the previous post. So it should be an expected of 0.7853982 radians. Everything else is in radians.
"where do you think the rounding error is happening?" My suspicion is during the multiplication and division calculations. Multiplying very small numbers with one another results in even smaller numbers.
For example, test 1: Velocity is { 0 AU, 6.684587e-9 AU } (Note: 6 digits of precision) https://github.com/Pulsar4xDevs/Pulsar4x/blob/Master/Pulsar4X/Pulsar4X.ECSLib/Helpers/Vector4.cs#L220 So velocity.Length() = Math.Sqrt(LengthSquared()); velocity.Length() = Math.Sqrt((X X) + (Y Y) + (Z Z) + (W W)); Y = 6.684587e-9. (6 digits of precision) Y * Y = 4.4683703360569e-17 (13 digits of precision)
At least the above example is what I suspect is where the problem lies. I still need to step through the code line by line and see where the precision is being eaten up. But see how multiply and divide can quickly make the needed digits of precision scale up.
ah yes, I'd forgotten that velocity is going to be pretty small compared to the distances
Ok I thought the reason that AU was chosen was because using Km would limit the size of the star systems too much. but doing the math now, if we store distances as meters in doubles, we can still have systems that are: 1.2016836379091028E+297 Au in radius or: 1.9001631454016871E+292 light years in radius. am I doing the math right? au = double.MaxValue / 149597870700 ly = au / 63241.077
so not quite 600ly in diameter I'm dumb. scientific numbers always screw with my head. it's like nearly 600 zeros.
if this is the case, lets just change everything to meters.
I must be doing the math wrong somewhere there.
Your math is right. Even with decimals (more precision, less overall range) if you work in meters you can still have systems with positions between -1E+12 light years to +1E+12 light years on both axis, in the system.
I would still strongly suggest decimal and km. The additional digits of precision would handle the large differences (in scale) between position and velocity much better.
Take the following example where your position is close to earth. That is around 1 AU from zero point (the sun). That is around 150 billion meters for your position values. That is already 12 digits of precision needed if you are comparing it to m/s velocity. (The cross products steps.) Double only allows for 15 digits of precision. That is a rather tight margin of 3 digits of precision lee-way.
Decimal would allow for 28 digits. I.e. you only start running into the same tight margin once you are comparing positions > 1 million billion billion meters (> 1 Septillion meters) to velocities in m/s.
My biggest concern with using decimals is the performance hit. (though in all honesty I don't know how big that will be) Second concern is they're a little more annoying to work with than doubles, Math doesn't have trig functions for them. we'd have to find (or write, eucgh) math for them.
if I understand the problem correctly, it's only where we're doing calcs with velocity, the bits of code that need to go fast don't need velocity, Position is calculated without velocity when doing kepler to state positions, where we need accurate math with velocity is when we're trying to convert from a state vector to kepler. and currently iirc that's the only time we're bothering to calculate the velocity (other than maybe for the ui) which happens relatively infrequently. It might be better to cast to decimal for those calculations.
Actualy no, there is other times we're using the velocity, and that's when we're in a hyperbolic trajectory (couldn't figure out the kepler math for that so it's just doing a Newtonian trajectory). I'd prefer to keep that as fast as possible performance wise since they have to be calculated frequently or they'll loose precision. (all the kepler stuff we can get an 'accurate' position for any given datetime)
I'm pretty keen to change the units we're storing everything in now though. why do you suggest km over m? we're going to need to do velocity in the same units as position/distances.
I'll have a look at casting the values in the calcs tomorrow and see what happens.
Did a quick test by calculating the Kepler elements when using Km instead of Au. (No change in calculation, just in what is passed in.) Doesn't seem to help much for the precision issue.
And yes, after taking a look, switching to decimal is a lot of work.
Also, beside the precision thing there is still the other failing tests, which seems to be a different problem.
My biggest problem with the tests is that it was hard to pinpoint what was wrong, and what was wrong because it was relying on something else that was wrong. there's a lot of circular interdependence going on.
Err, I probably should have mentioned, the tests in OrbitTesting2.cs are a bit cleaner than the tests in OrbitTests and attempt to target specific functions while incrementing through a given orbit. Currently it only tests a circular orbit, but there's data there to test elliptical, retrograde, and 3d retrograde. I think I need to look at TestCaseSource Attribute to get it testing the other types of orbit, I'd not gotten that far due to it not even passing the circular one.
I should maybe have not used the Latin symbols in the vars there, but I'd just figured out it could be done so I kind of had a bit of fun with it.
I've implemented the TestCaseSource Attribute in the OrbitTesting2 tests. I've also got it calculating the LonditudeOfPeriapsis from the OrbitMath calc - which takes inclination into account and flips the sign on the ArgumentOfPeriapsis if it's a retrograde orbit, I think it's passing a few more tests. I expect this to cause problems on the last test set though. (its only good for 2d).
I've changed vector4 to vector3.
the fourth dimension if memory served was to make matrix math easier, however we're not doing that at all currently, and if we do end up doing that I can't imagine it'll be often enough to be worth having the fourth dimension around doing nothing.
Also added some vector3 "Precise" functions ie, Vector3.MultiplyPrecise(vector3, vector3) which will cast the given vector3s to decimals and return a decimal tuple.
so far there are:
MultiplyPrecise
DividePrecise
DotPrecise
CrossPrecise
MagnatudePrecise.
Note that Magnatude returns a double as it casts back down after the multiplication to do the sqrt.
if they end up being used usefully I'll look at creating a decimal vector3.
I've not written tests or anything for them though...
I was intending to put them in another namespace. ie pulsar4x.Vectors.Vector3, however I've left them as is so as not to make merging with your stuff difficult. (the current change shouldn't be too bad your end, a find and a replace will take care of most of it, and removing the last item anywhere calling the constructors.
Just found this which looks written fairly strait fowardly: https://downloads.rene-schwarz.com/download/M002-Cartesian_State_Vectors_to_Keplerian_Orbit_Elements.pdf
If we figure out the True Anomaly calculations, (those seem to be out by 180degrees sometimes)
we can drop the EccentricAnomalyFromStateVectors calculations, and calculate EccentricAnomaly from the TrueAnomaly. I think that might solve most of our problems.
Edit: no it's worse than just 180degrees out. I might throw the results into the graphical orbit icon lines widget and maybe I'll get some useful information out of that. maybe it's doing something stupid like calculating it from the reference direction instead of the AoP. Might do the same for the EccentricAnomaly calcs as well, might see something interesting.
I've been working on getting tests for sub-components and function in.
So here is a fun tidbit. The maths do not work out unless you use SI units, so all distances must be in m and velocity in m/s.
Specifically these tests for CalculateAngularMomentum fail unless done in SI units. https://github.com/UberWaffe/Pulsar4x/blob/experimental/OrbitTests2/Pulsar4X/Pulsar4X.Tests/OrbitTests.cs#L186
In hindsight this should have been obvious to me. As once you do the cross product, your SI unit changes. m m = m^2. So km km = km^2. To convert from km^2 to m^2 is a factor of a million, not 1000. So if the maths are done in numbers other than the SI units, you need to also remember that when converting you are working with squared factors and such. Easier to just stick to SI units.
Still doesn't fix all of the tests though. But certainly one of the issues lurking about.
Edit: Also be ware when getting conversions and such from websites. Most have the same rounding issues for large values (floats in javascript). So when setting up tests you cannot really check for precise numbers.
Himnnn thats a good point, are we actualy doing any calcs where that's going to come into play though? I think everything is done in AU at the moment, but I guess there could be something somewhere where converting up and down is being done wrong. but yeah another good reason to change everything to same units.
Check out my latest push, and have a play around with the widget I mentioned above, I think you'll notice some interesting things with eccentric and true anomaly calcs, namely: Some calcs are right (or at least very close) if there's some eccentricity (or maybe it's not having a LoAN/AoP of 0, not sure) Some calcs are right, until the angle gets above Pi, then it starts decreasing. Examine both the Ensuing Calm (Circular Orbit with LoAN, AoP etc 0) and Mercury which is Eccentric and has LoAN and AoP angles. Note if you hover the mouse over the value name in the ui it'll highlight the arc on the map. Examine them while time is running, note you can also open up the time widget and decrease the tickrate fro 0.5 seconds to 0.1 seconds (ie it'll run turns every half second default, change it to 0.1)
Awesome. WIll take a look. Currently setting up tests for the longitude of ascending node (LoAN / Omega) and it fails for most cases. Not sure why yet.
So working up the energy to change all the units, pros and cons of km vs m? radians instead of degrees is also on my list. the number of bugs because I've used degrees is silly.
"pros and cons of km vs m?" The formulas nearly everywhere (in the literature) is in m. So I would go with that. If you go with km, you introduce a lot of conversion issues. (I've been running into hoards of those while trying to setup tests.)
Quick note. I've created some SI unit classes for angles and distances. "https://github.com/UberWaffe/Pulsar4x/tree/experimental/OrbitTests2/Pulsar4X/Pulsar4X.ECSLib/Helpers/SIValues" (A double with an SI unit). It makes it a lot easier to keep track of what unit the value is in (since you can just specify it in whatever, and where it is used it pulls the correct unit value via GetXYZ() ). Example of its use: https://github.com/UberWaffe/Pulsar4x/blob/experimental/OrbitTests2/Pulsar4X/Pulsar4X.ECSLib/Helpers/GameMath.cs#L61 The class knows what it stores it in, and when you create it, it auto converts. When you request the value, it spits it out in whatever you request. So you have a lot less finger trouble.
They might be useful to use instead of just doubles. (Not sure.)
They've been helpful in the tests for me. (Busy playing around with a simple orbit processor (plain circular orbits).)
Yeah I've been considering creating a struct for Radians at least - and maybe set a MaxValue of 2pi and min of -2pi, and check these on add/subtract, mutliply etc, and roll the value around. Maybe.
For distances... I'm tossing that up, if I did I'd name it as what it was stored in (ie Meters), as that it helps with conceptualize how big (MaxValue) you can go and also how small you can go before running into rounding errors.
I've fixed one of the EccentricAnomaly calcs. if it was given a negitive True anomaly it needed to return a negitive result. that one was easy, but got to figure out the others.
Both:
TrueAnomaly(Vector3 eccentVector, Vector3 position, Vector3 velocity)
TrueAnomaly(double sgp, Vector3 position, Vector3 velocity)
Returns an incorrect value in the circular test orbit, but 'normal' orbits seem fine, nor sure which value it's having problems with. the former is probibly breaking if the eccentricityVector is 0.
GetEccentricAnomalyFromStateVectors2(double sgp, double semiMajorAxis, Vector3 position, Vector3 velocity)
has the same problem.
GetEccentricAnomalyFromStateVectors(Vector3 position, double semiMajAxis, double linierEccentricity, double aop)
is only correct in the circular test orbit but sign flipped half the time even in that case.
I think if we can fix any one of the first three above functions we can shuffle the kepler from state vectors around and have stuff 'work'
TrueAnomaly using EccentrictyVector I don't think is having the conversion problem you mentioned, as it's all being done in au currently, and not converting. don't think I've got any tests on the eccentricity vector.. Maybe I need to bite the bullet and just do the unit conversion work.
The unit conversion will only fix the accuracy being eaten up. The incorrect calcs are likely unrelated. I would focus on getting the maths working first before worrying too much about accuracy.
I've not made much progress. Having trouble finding calculations to use as test cases for the code. (The one website that spits out calculation steps doesn't use the same calculation for LoAN.)
Yeah, true.
Might have to test LoAN it by setting up a kelper orbit, then testing the state to kepler conversion against that. (which is how the calculations are being tested in OrbitTesting2.cs - it is testing it in one place, but is failing on some orbits, but passing the circular one. )
I'm currently looking at the eccentricity vector, I suspect that may be wrong. Edit, yeah it's all over the place in our "test" orbit. but seems ok in a real one looking at the graphical widget.
I'm an idiot, of course EccentricityVector calc is going to be 'confused' in a perfectly circular orbit. it's a completely conceptual value in such an orbit. Probably the problem you're having with LoAN as well. however due to other angles being calculated from these, we need to detect when we have a near circular (within some value of precision) orbit and set them to 0. instead of having them flip all over the place with floating point errors.
The problem I am having with LoAN is similar, in that when your inclination is zero (which I think is always the case in 2D?) the LoAN is arbitrary. I've (in my branch) locked the LoAN to 0 in that case.
(I moved the LoAN calc to a function.) https://github.com/UberWaffe/Pulsar4x/blob/experimental/OrbitTests2/Pulsar4X/Pulsar4X.ECSLib/Helpers/OrbitMath.cs#L164
you'll have to check for an inclination of 180 as well (signifying a retrograde orbit)
I've sent you an invite, you can now push directly. Push what you've got. I've fixed one of the TrueAnomaly calcs (though it's still failing a test which is odd, test thinks its wrong by 180 degrees. need to check that before I push)
So I messed with some problems with the eccentricityVector, tried doing it in Km, though the accuracy difference wasn't all that significant. then I realized, it was passing the eccentricity test, but not the retrograde eccentricity test. so yeah, the velocity vector/heading calc is messing up in retrograde orbits.
Actualy, thinking about the checkign for an inclination of 180 on the LoAN migth be wrong... I'm currently looking at the heading. in the widget I've now got it displaying in the correct direction... but where I'm doing that may be in the wrong place. for example, the OrbitMath.GetHeadingFromPeriapsis returns the heading ralitive to the periapsis (and on the same plane if it were an inclined orbit) as such, a body at its periapsis and in an orbit orbit with a periapsis of 0 and heading +90 (up towards the Y axis) is neumericaly going to be the same as a retrograde orbit wrt the heading... ick this is confusing. Best way to think about it is drawing on orbit on a piece of paper, and folding the paper along the LoAN, the angle of the fold is your inclination. when we go to render it, we're taking this information (i should be 0 or 180) and taking a picture of it from top down.... if I'm rendering angles and arcs, I can't draw them raw in a retrograde orbit, I have to flip them... I think. (if I try flip them in the math which I think I'm doing in some cases, it messes up other math)
Glad you are making some progress. I haven't had time this past few days to work on this. Will likely have time this weekend to really dig into it again.
Yeah, slowly. I understand what the problem is here, but... not how to fix it. since my visual tool doesn't work in this situation I can't see where it's going wrong...
Currently the OrbitMath calculations for veclocity vector return a vector that's ralitive/is on the same plane as the AoP. We need a calc that will give us the velocity vector in the parent bodys frame of reference. I think that's the best way to fix it. OrbitProcessor sort of tries to do this in a 2d way, but it's incorrect, i think. also we're going to have to go the other way too, velocity vector used to get kepler elements needs to be translated from a world Cartesian vector to it's body ralative vector.
Ok I think I've got it, what I need to do to convert from orbit reference to world coordinate vector is:
I could probably add the AoP in at the heading function, but having it referenced from the AoP rather than the LoAN seems more 'correct'. I guess worldspace polar coordinates don't even make sense without including the LoAN and inclination.
Graphically I should be matrix rotating everything around the LoAN by the inclination. then if I ever want to 3dize it all I need to do is 3dize the camera.
Made some very minor progress, got a matrix3d class in now, and cherry picked some of your LoAN work. I may have gotten distracted after buying Kenshi.
I'll pull your latest code and check. I didn't make any progress over the weekend.
after messing with it more, I'm kind of stuck, I thouhgt I had the whole thing figured out but... I'm thinking around in circles. everything apears to work/pass tests in normal orbits, but retrograde/inclined orbits fail right at the beginning where eccentricity vector magnatude doesn't equal the original eccentricity, which leads me to suspect that I'm still screwing up the velocity vector in such cases. (starting with a kepler orbit, and translating into state vectors).
Maybe this is what I'm trying to do here: https://www.orbiter-forum.com/showthread.php?p=311667&postcount=10 he's using a quaternion, to do what if I'm understanding him correctly what I was trying to do with a matrix rotation.
I've solved it.
Most the tests are passing now in the Orbit2 tests.
I was on the right track with rotating inclination on the axis of the LoAN.
needed to do MatrixRotateZ(-LoAN) MatrixRotateX(-inclination)
that seemed to do the trick. though the - inclination doesn't seem right to me, I feel like that should be a positive, ~~I think my matrix rotation func is correct on XRotations, a 90degre X rotatation should change a 0,0,1 to a 0,1,0 right?~~ edit: nope it should be 0,-1,0 since right hand rule thumb should point positive* x not negative. duh.
I still need to do a bit of tidy up and check some funky stuff with the AoP calculations, there's one test there that is not passing, but I can probily get rid of that calculation. however there's some functions there that are passing, but in the ui are sometimes flicking to odd angles or NaN now and then. There's another test there that I've not looked at which is basicaly doing all the previously passing tests, I think it's just a case of 360degrees vs -360degrees not being equal. And then there's some rendering fixes I need to do to bring the graphics side into line with the math/physics fixes.
Ok here's a question: Should I normalize TrueAnomaly, AoP, LoAN, functions etc etc to return a positive angle?
What would be the benefit to normalizing? Wouldn't you then have to include some additional direction field?
mostly consistency, shouldn't need a direction field unless we're talking about angular momentum ,in which case it wouldn't make sense to normalize to a positive. (or at all).
Well I crossed a bunch of things off the list. all the OrbitTesting2 fuzz tests are passing, though I think there's still a problem with the AoP calc in some cases, which leads to a jump in position in non newtonion translations to new orbits sometimes.
If you want to go back to messing with the cargo/storage stuff feel free, I can continue tweeking and improving this slowly, its a whole lot less broken than it was.
I think I may have fixed the bug that was popping up with retrograde orbits being broken. maybe. I never managed to define exactly when this was popping up, but I recently found and fixed a bug where inclination was getting calculated via atan2(r, z) instead of atan2(z, r) z should be zero for the game, but most the orbit calcs can do 3d. the above was giving the orbits a bit of inclination when the inclination should have been zero, which was screwing stuff up.
Most the bugs with this are now fixed, I think the remaining bugs are all ui side now. There are still a few failing tests, though they may be problems with the tests themselfs, needs checking.
Outright bugs in the Movement and Orbital code:
Retrograde (clockwise) orbits are severely broken.- graphical representation of the orbit is prograde (this could be fixed easly but not sure with how things will change due to other retrograde orbit problems, might be a better way to do it)Velocity Vector is incorrectPrograde orbitsFixed.Retrograde orbits still broken.Needs tests.Prograde orbitsFixed.Retrograde orbits still broken.Needs tests.Hyperbolic (newtonion move) orbit ship icon's heading is not correctFixed in 91ce915bd6432fe1ecd7aa072a61d43817e6e5f1219 Translation to new orbits sometimes jumps to new position.
Prograde orbitsFixed??? needs confirmingTranslation widget/position incorrect exit coordinates. (should be easy)Fixed in commit 1682537c19544275213ac264c9bba5c329f4cf6a mouse world coordinates were incorrect.See #223 for Todos/nice to have/lower priority stuff: