www.algodoo.com

Insipired by Kerbal Space Program's map I'm trying to plot the orbit of an moving object subjected to the gravitational force. As you can see in the scene (attachment), the orbit is pretty much accurate, i'm using this equation to represent the orbit and I'm getting all the necessary data through the readable [x,y] and vel converting then to polar coordinates.

I've been having some problems with the mass, when I started the Algodoo again the mass value would've changed and that's why I've defined it, but it shouldn't be like this.

My main issue with the scene is on the axis of the conic. As you can see, when the object gets near the apoapsis or periapsis of the ellipsis, the predicted orbit wobbles. I've been monitoring the numbers and I think the problem is here:

Code: Select all

parameter := correct_sgn(k, 0) * math.acos(correct_magnitude(((p / radius) - 1) / ecc));

Specially with the

Code: Select all

correct_magnitude(((p / radius) - 1) / ecc)

We know the arccos(x) function only accepts -1 < x < 1, I saw that this ((p / radius) - 1) / ecc was getting values > 1 or < -1 near those two spots, so I limited it with the correct_magnitude(x) function to not bug the scene and see what's happening. Why it is getting those values? Floating point error? I was suspecting about how I calculate the p or the ecc. First the ecc varies in some range but it was supposed to be constant, but by manually inserting the eccentricity value the wobble didn't disapear. As p is somewhat constant, i've changed its value, but the problem with it is that though it removes the wobble, it makes the orbit not to represent the real orbit. So now I have no clues anymore...Any help from these thyme guys over here?

Sorry if I've misspell something.

Edit: I also would like to tell you that my main goal is have multiple objects (Earth and Moon, for example) and make it like a game. I believe I have two ways of doing it, analitically with Patched Conic Approximation(just like KSP) or numerically by solving x and y ODE, but this is a future project...

Anyone

No one knows...

It's hard to analyse your a bit chaotic code... /-)

First what I see - try to use math.atan2 function instead of manually find angle as atan(x/y).

How do I add comments? I've always felt the necessity to do that. I may split the code to make it less chaotic.

There is a good way - self explaining code

BTW, I've got the idea, I think I could help you.

Kilinich wrote:There is a good way - self explaining code

BTW, I've got the idea, I think I could help you.

I can't make a self explaining code. There are too many variable names, and long ones. If i'd add those names the code would be even more chaotic. I'll make some improvements on some variables name and will upload it soon.

Edit:

First thing to say is that θ₀ is the angle of the tilt of the orbit, it comes from solving r=p/(1+e cos(θ-θ₀)) for θ₀ getting θ₀=θ-arccos((p/r-1)/e) and i've called the arccos((p/r-1)/e) as parameter;

r_dot_ish tells if the radius is increasing or decreasing, this correct_sgn(r_dot_ish, 0) is necessary because the parameter needs to change sign depending on the object location (above or below the orbit's axis)

It's harder than I though )
I think there should be another way with no arccos calculation...

The problem is definitely in floating-point calculation, in your code it's became very noticeable, somewhere up to 0.001
Than happened because of many (x^0.5)^2 operation and additional variables.
Try to optimize calculation to a single function, simplify it (maybe using "mathematica" or similar software).

BTW: You could use clamp(x, min, max) for value limitation...

Since the sum of the energy is a constant, i tested epsilon = -2.9 fits quite well,
the calculation of ecc simplifies a bit, also the reduced mass can be set to a fixed value.
Be also aware of the modification of the two-body-problem-transforming to a single-body-problem
by usage of reduced mass and radius-vectors, and maybe herein lies the solution to your problem;
it could be that the two systems are not correctly transformed into each other.
Here specially L=I*theta_dot might be a place to search.

cant afford more time right now.... hope that helped a bit.
cu