Project Kepler Orbit Simulation Toolkit

[cjp]

KOST is the Kepler Orbit Simulation Toolkit. But it won't cost you anything .

This is a continuation of this thread:
http://orbiter-forum.com/showthread.php?t=4518
Once I found out that the thing I wanted did not exist, I decided to make it myself. I now have a 'zeroth' version, which doesn't do very much, but I'm releasing it anyway to allow for API discussions before going on with the real stuff. Because of this release, I think "Addon Development" is a more appropriate place for the discussion.

KOST is available on [ame="http://www.orbithangar.com/searchid.php?ID=3825"]OrbitHangar[/ame].

To give an idea what it will be: KOST will be a small open source toolkit for doing orbital calculations. It will be very useful for making MFDs and things like custom planet DLLs. While it can be compiled to use Orbiter's data types natively, it doesn't depend on the Orbiter API, so it can also be used for standalone applications.

 

[cjp]

I'm having trouble with orbital elements calculations (please help me). It's probably a minus sign error.

I'm currently developing a remake of Orbit MFD as a test case for KOST. For this, I need a function to translate orbital state vectors (position+velocity) into orbital elements (inclination, longitude of ascending node, etc.).

I'm using appendix C in Orbiter.pdf as a reference. Unfortunately, the coordinate system used in the appendix seems to be 90 degrees rotated w.r.t. the coordinates of the state vectors I'm getting from Orbiter. So I applied the following translation to all formulas to correct for this:
x -> x
y -> z
z -> -y

I'm getting some things right, but not everything. In the attached image, you can see my Orbit MFD on the left, and the built-in Orbit MFD on the right. The debug string displays the state vectors.

The first thing that goes wrong is the LAN (although in the image it's correct, as I had a second minus sign bug which compensated the still existing bug). So I tried to calculate everything by hand:
_____________________________________________
mu = G * Mearth = 398600 = 1 000 000 * 0.3986
1/mu = 10^-6 * 2.5088

r = (-5584109, -3916727, 720329)
= 1 000 000 * (-5.58, -3.92, 0.72)

Rad = |r| = 6.86M (OK)

v = -3365, 2664, -6478
= 1 000 * (-3.37, 2.66, -6.48)

Vel = |v| = 7.77k (OK)
_____________________________________________
h = r x v
hx = ry*vz - rz*vy = 1 000 000 000 * 23.45
hy = rz*vx - rx*vz = 1 000 000 000 * -38.60
hz = rx*vy - ry*vx = 1 000 000 000 * -28.05

n = -y x h = (-hz, 0, hx) = 1 000 000 000 * (28.05, 0, 23.45)

e = (1/mu) * [(|v|^2 - mu/|r|)*r - (r.v)*v]
|v|^2 - mu/|r| = 1 000 000 * 60.37
r.v = 1 000 000 000 * (18.80 - 10.43 - 4.67) = 1 000 000 000 * 3.7088
(|v|^2 - mu/|r|)*r = 10^12 * (-336.9, -236.7, -43.5)
(r.v)*v = 10^12 * ( -12.5, 9.9, -24.0)
________________________________ +
10^12 * (-349.4, -226.8, -67.5)
e = 2.5088 * 10^6 * (-349.4, -226.8, -67.5)
= 1 000 000 * (-876.6, -569.0, -169.3)
_____________________________________________
Inc = arccos(-hy / |h|) = arccos(0.7260) = 43.45 degrees (OK)

LAN = arccos(nx / |n|) = arccos(0.7672) = 39.89 degrees (BUG)
check: nz < 0: not the case

All in SI units.

Does anybody see where it goes wrong?

 

[tblaxland]

Orbiter uses a different coordinate system to that used in Appendix C. In Orbiter, you have a left-handed co-ordinate system, x points to the vernal equinox and y points ecliptic north and z is determined by the left-hand rule. You do not need to do any coordinate conversion, just remember:

h=r x v
n=y x h

I haven't double checked it with calculations for your particular state vectors. I did verify it by left hand rule (don't you wish you had two left hands sometimes!) so you should be OK with that.

 

[mjessick]

What I do is translate Orbiter's left handed coordinate sets on input and output from my code. For the usual orbit mechanics right handed frames with Z along the rotational axis and X in the known direction, all you do to convert is swap the Y and Z components. (Oddly, the algorithm is the same going either way, which is nice )

So try it as you had it above, but without the minus sign.

 

[cjp]

In Orbiter, you have a left-handed co-ordinate system

Wait a minute... Appendix C says: "We assume a right-handed system". So, there isn't just a rotation, there's also a mirroring of the coordinate system!!!

This leaves me wondering what would be the standard mathematically speaking. Wikipedia says:
"The standard orientation, [..] is called right-handed or positive."

I think it's best to follow the most usual standards in KOST as much as possible, and convert the state vectors from Orbiter's unusual system to the standard system before giving them to KOST. Then I'll see if that works.

You do not need to do any coordinate conversion, just remember:
h=r x v
n=y x h

I used to have that, but it was inconsistent with my z -> -y translation, so I changed it. BTW, don't you think the other places where coordinates are mentioned explicitly should also be changed? For instance:

Inc = arccos(hz / |h|)
LAN = arccos(nx / |n|)
and all those rules like if(ny < 0) then ...

What would be the correct transformation from Appendix C to Orbiter? Maybe this:
x -> x
y -> z
z -> y

That mirrors vectors through the y=z plane, so it changes handed-ness.

 

[Hielor]

I think it's best to follow the most usual standards in KOST as much as possible, and convert the state vectors from Orbiter's unusual system to the standard system before giving them to KOST.

Yup, martin picked the left-handed coordinate system because it was least used, and therefore most likely to cause problems!

Kidding, obviously. I'm pretty sure the reason that Orbiter uses left-handed coordinates is because DirectX uses a left-handed coordinate system for rendering.

Why? I don't know. If you want, I'll see if I can find the right person on the DirectX team to ask when I get to Redmond in January. Watch, the explanation will be something like "well, we had drawn it on the whiteboard as a right-handed system, but the camera was accidentally held upside-down for the photo of it we used later in the spec."

 

[cjp]

Why? I don't know. If you want, I'll see if I can find the right person on the DirectX team to ask when I get to Redmond in January. Watch, the explanation will be something like "well, we had drawn it on the whiteboard as a right-handed system, but the camera was accidentally held upside-down for the photo of it we used later in the spec."

I'm sure it has something to do with the difference between mathematics (where the y-axis goes up) and computer graphics (where the y-axis goes down). In CG, the y-axis goes down because in a CRT monitor, the electron beam scans the screen downward. The electron beam goes to the right and downward because that's the standard reading direction in western languages.

But of course I'm interested in the real story. Whatever the real reason is, it is another example of what I see as the dirtiness of DirectX, compared to the (right-handed) OpenGL.

 

[tblaxland]

I think it's best to follow the most usual standards in KOST as much as possible, and convert the state vectors from Orbiter's unusual system to the standard system before giving them to KOST. Then I'll see if that works.

Fair enough. Right-handed is definitely more common.

What would be the correct transformation from Appendix C to Orbiter? Maybe this:
x -> x
y -> z
z -> y

Correct. Once you have done that transformation you can do the calcs as per Appendix C.

 

[cjp]

Correct. Once you have done that transformation you can do the calcs as per Appendix C.

It seems to work now. Thanks!

Inc is correct
LAN is correct
LPe is 360 degrees too much, but I can correct for that
MnL is still insane, so that's the next item to work on.

I'll soon have exactly one mode working (FRM=ECL, NT, MOD=no visualiation, DST=radius). For the rest only visualisation is non-trivial.


-----Posted Added-----


I now have all parameters correct (= same as built-in Orbit MFD) for elliptical orbits. I have support for DST=Alt mode (still have to test its accuracy for non-spherical bodies) and partially for setting a target.

I'm having a new issue: hyperbolic orbits. How can I calculate PeT for them? For elliptical orbits I simply used MnA and T, but T is not available, and MnA does not contain a useful value (Orbit MFD says MnA=0.0). Still PeT is a very useful value e.g. when arriving at your destination planet, so I really need it.
(Edit: already found a solution that works)

 

[tblaxland]

I'm having a new issue: hyperbolic orbits. How can I calculate PeT for them? For elliptical orbits I simply used MnA and T, but T is not available, and MnA does not contain a useful value (Orbit MFD says MnA=0.0). Still PeT is a very useful value e.g. when arriving at your destination planet, so I really need it.
(Edit: already found a solution that works)

Hmm, seems some posts just got eaten, I'll try again:

Can you post your solution? I found this one:
http://books.google.com/books?id=2I...&hl=en&sa=X&oi=book_result&resnum=1&ct=result

 

[cjp]

Hmm, seems some posts just got eaten, I'll try again:

Can you post your solution? I found this one:
http://books.google.com/books?id=2I...&hl=en&sa=X&oi=book_result&resnum=1&ct=result

I can't read your Google book. It gives some error on the book being unavailable or me reaching my viewing limit.

I have something that seems to work, but I'm not 100% sure it is correct. Can you confirm this:

From kost_elements.cpp:

        params->SMi =
            sqrt(elements->a*elements->a* (elements->e*elements->e - 1.0));

        params->EcA = acosh((1.0 - absr/elements->a) / elements->e);


        /*Copy sign from sin(TrA)*/
        if(sin(params->TrA) * params->EcA < 0.0)
            params->EcA = -params->EcA;

        params->MnA = elements->e * sinh(params->EcA) - params->EcA;

        params->T = 
            M_TWOPI * sqrt(fabs(elements->a*elements->a*elements->a/mu));

        kostReal tPe = params->MnA * params->T / M_TWOPI; /*Time since last Pe*/

        params->PeT = -tPe;

 

[tblaxland]

I can't read your Google book. It gives some error on the book being unavailable or me reaching my viewing limit.

I have something that seems to work, but I'm not 100% sure it is correct. Can you confirm this:

Looks the same except for the calculation of the hyperbolic anomaly, EcA. What is "absr"?

 

[cjp]

Looks the same except for the calculation of the hyperbolic anomaly, EcA. What is "absr"?

absr is the length of the position vector.

What is the alternative way of calculating EcA?

 

[tblaxland]

absr is the length of the position vector.

What is the alternative way of calculating EcA?

What you have looks OK for EcA then.

An alternative that doesn't require absr, derived from eq 4.155 of the above link:

params->EcA = 2.0 * atanh( tan(params->TrA / 2.0) / sqrt( (1.0 + elements->e) / (elements->e - 1.0) ) )

 

[cjp]

What you have looks OK for EcA then.

An alternative that doesn't require absr, derived from eq 4.155 of the above link:

params->EcA = 2.0 * atanh( tan(params->TrA / 2.0) / sqrt( (1.0 + elements->e) / (elements->e - 1.0) ) )

Then I think I'll keep mine. I already have absr, and I don't like calculating an atanh, a tan and a sqrt instead of just an acosh, if I don't have any clue whether it will be more accurate.

 

[tblaxland]

Then I think I'll keep mine. I already have absr, and I don't like calculating an atanh, a tan and a sqrt instead of just an acosh, if I don't have any clue whether it will be more accurate.

No more accurate, just different (EDIT: different method, that is, same answer).

 

[cjp]

I attached a new version of KOST to the first post in this thread. Basically, it's just the same code as used in Free Orbit MFD, but with different default settings, and some modifications to make it work with gcc -ansi -pedantic.

There is one big change: I added a lot of documentation.

So, what do you think? Is it worth placing this on Orbiter Hangar? Is the documentation clear enough? Are you missing certain features? Comments on the existing features?

 

[tblaxland]

So, what do you think? Is it worth placing this on Orbiter Hangar? Is the documentation clear enough? Are you missing certain features? Comments on the existing features?

I'll have a look over the weekend.

One question, does KOST include the ability to output true anomaly as a function of time? No big deal if it doesn't, because I was planning to write this function for another project anyway. Perhaps it could be added to KOST if you don't already have it.

 

[Enjo]

cjp: I'll definitely use it when I have too, but I for example don't have time to evaluate just now, so don't let the low response level stop you from releasing/developing it. It's not for the majority of end users, but for minority of developers mainly (same applies for FreeOrbit btw.).

 

[tblaxland]

I'll have a look over the weekend.

One question, does KOST include the ability to output true anomaly as a function of time? No big deal if it doesn't, because I was planning to write this function for another project anyway. Perhaps it could be added to KOST if you don't already have it.

Sorry it took me a bit longer than expected to have a look at this.

This is a very useful toolkit, well done. :speakcool:

Since true anomaly as a function of time is not included I'll write up a function that suits my purpose. I think this would be a useful inclusion in KOST. It may be a couple weeks before I get to this, though.

One thing that is on my mind is that because it involves recursion, I think I should make it so that each iteration can be calculated on separate time steps so as to not bog down Orbiter. Also, since each iteration increases accuracy, an arbitrary level of accuracy could be achieved provided you are prepared to wait enough time steps for it. As I understand it, you want to keep KOST as generic as possible, so I'll have to think on the best way to achieve this without getting any Orbiter code in there - perhaps if KOST just includes the "inner" part of the iteration and the calling function handles the recursion.