calculating change in eccentricity

[jedidia]

I got a quick and dirty math question, to test out something:

When I apply a certain amount of delta-v at periapsis, how do I calclate the change of eccentricity?

 

[Zatnikitelman]

If you can calculate your new apoapsis, then you can plug the values into the equation here relating apoapsis, periapsis and eccentricity http://en.wikipedia.org/wiki/Orbital_eccentricity#Calculation.
You could also use your delta-v and perhaps plug it into the speed at apoapsis/periapsis equation: http://www.rclsoftware.org.uk/gravel/orbits.html

 

[jedidia]

thanks. That brought me along a fair bit, but the predictions of my little calculator aren't accurate enough... got to do with the periapsis moving quite a bit of course, since the delta-v doesn't get aplied all at once. I'll have to split it up a bit...

What I'd need now is a possibility to calculate gamma (angle between radius and prograde direction) given sma, eccentricity, periapsis and time to periapsis. Let's see what I can find...

 

[tblaxland]

What I'd need now is a possibility to calculate gamma (angle between radius and prograde direction) given sma, eccentricity, periapsis and time to periapsis. Let's see what I can find...

This will do the trick: http://en.wikipedia.org/wiki/Elliptic_orbit#Flight_path_angle

 

[jedidia]

*sigh* me and teh math... I'm in trouble!

Maybe I should just write down the whole problem, because the longer I look at it the more I get the impression that I'm trying to mount the horse from the rear end...

Basically, I'm trying to program a little calculator to help me along with eject burns on low thrust. The auto-burn mode of IEATMFD burns a third orbit around periapsis, which works quite well. The goal of my calculator is to calculate how long it will take until orbit ejection is reached (i.e. eccentricity gets bigger than 1).

I have a basic version functioning, but it's got severe accuracy problems, because of above mentioned lack of PeR-modification. Allthough the predictions done by the calculator currently are already of some help, I can't get rid of the feeling that it might be alot more help with the right formulae.

Currently I have an iteration that goes through the orbits, aplies the delta-v burned in one third orbit at periapsis, calculates the new eccentricity and orbital period and does that until eccentricity goes above 1. The major improvement I want to add is to predict the modification of the periapsis.

However, I can't see how I should go about doing that, since I need PeR to calculate eccentricity and SmA, so I need some other basis to calculate the deviation of the PeR.

So what I'd need is a formula that lets me calculate the new PeR based on the old PeR, eccentricity, delta-v aplied and time of burn, but I can't find any. Anyone got a hint for me?

 

[tblaxland]

So, you starting thrusting T/6 seconds before periapsis and continue thrusting along the velocity vector until T/6 seconds after periapsis, where T is the orbital period?

The only way I can see to solve that is to resolve the orbit parameters at small intervals during the burn (every 1 s?). Whether or not that is easier than just numerically integrating the full trajectory, I'll leave for you to decide.

If you do decide to go down that route, you could use [ame="http://www.orbithangar.com/searchid.php?ID=3825"]KOST[/ame] to propagate your orbit parameters for you between each time step. Then, at each time step, convert your orbit parameters to state vectors, increase the magnitude of the velocity vector by however much dV you are predicting the engine to provide in that time step, and then convert the state vectors back to orbit parameters.

A temporary note of caution on KOST: it has some small bugs which cjp and I are working on ATM. These relate only to some very specific cases of equatorial orbits, orbits with eccentricities very close to 1, and hyperbolic orbits far away from periapsis. That said, there is nothing to stop you using KOST functions as they are since the function definitions all remain the same. You can just plug in our bug fixes when they are ready.

 

[jedidia]

Have spent a day overworking the autoburn funtion of IEATMFD (it now tries to keep at least the longitude of the periapsis constant, seems to work fine so far), and two days in bed because of my sinuses... argh!

The only way I can see to solve that is to resolve the orbit parameters at small intervals during the burn (every 1 s?). Whether or not that is easier than just numerically integrating the full trajectory, I'll leave for you to decide.

well, it is easier becuase I don't know jack about integrations... :lol:

If you do decide to go down that route, you could use KOST
KOST
to propagate your orbit parameters for you between each time step.

oh, cooool! That might make it a bit easier for me, thanks a lot! I'll see what I come up with when I'm on my feet again.

 

[tblaxland]

oh, cooool! That might make it a bit easier for me, thanks a lot! I'll see what I come up with when I'm on my feet again.

You might like to check out the source code for [ame="http://www.orbithangar.com/searchid.php?ID=3824"]State Vector MFD V0.9.3[/ame] too (the StateVectorMFD::Update function does all the work with regard to KOST).

Another note of caution: KOST uses a right hand coordinate system so you need to convert to/from Orbiter's left handed system. That is covered in the State Vector MFD code too:

	// Convert from left-hand coords to right-hand
	std::swap(state1.pos.y,state1.pos.z);
	std::swap(state1.vel.y,state1.vel.z);

 

[jedidia]

I have kind of a different problem now, but I thought it might be best to put it into my old thread.

I have a rather weird phenomenon here, and I can't find the mistake. Basically, I'm trying to calculatethe time it takes from the current point in orbit to a specified point in orbit. The thing seems to work most of the time, but in certain cases, it fails, because the mean anomaly of my targeted point comes out smaller than the mean anomaly of my origin, although the true anomaly of the targeted point is way bigger. Here's one of the cases where that happens:

The Orbit has an eccentricity of 0.0812, and a semi major axis of 7.636e6 m.

The point the ship is currently at has a true anomaly of 2.6496 (all angles in radians). I calculated its mean anomaly by going over eccentric anomaly, the formula looks as follows and should be right (at least it returns the same result as is contained in the ORBITPARAMS at thatt point (TA stands for true anomaly):

Eo = arccos((e + cos(TA)) / (1 + e + cos(TA) ))

and

Mo = Eo - e * sin(Eo)

giving me a eccentric anomaly of 2.61 and a mean anomaly of 2.57 (matching my orbit MFD at that point)

I then do the same, with the exact same formula, for my targeted point with a true anomaly of 6.06 (only slightly before periapsis), and I get an eccentric anomaly of 0.206 and a mean anomaly of 0.19

I checked the operations in my code several times, the chance that I overlooked errors in the code are very slim. So how is it possible that I get a much smaller mean anomaly for a point so much further along the orbit??

 

[tblaxland]

Without running the numbers my self, I would say that you are taking the result of the inverse cosine in the incorrect quadrant. The sine of the true, eccentric and mean anomalies should always be the same for an elliptic orbit. So for true anomaly = 6.06, eccentric anomaly = 2 * PI - 0.206 = 6.077 and mean anomaly = 2 * PI - 0.19 = 6.093.

 

[jedidia]

The quadrants! darn, they haven't been mentioned in any of the example problems I came upon, so I kind of assumed that they could be neglected, at least when calculating with radians...

So, even when calculating with radians, I have to split a circle into 4 pieces of 90 degrees each?

 

[tblaxland]

So, even when calculating with radians, I have to split a circle into 4 pieces of 90 degrees each?

No, not quadrants, just semicircles. It is possible to have a true anomaly > PI/2 and eccentricity and mean anomaly < PI/2, but the inverse cosine function in the above formula can handle that. You just need to know whether the true anomaly is moving away from periapsis or moving toward it because taking the cosine of the true anomaly conceals that information.

 

[jedidia]

ok, here we go! thanks a lot!

 

[Miner1]

There's a better formula for converting between eccentric anomaly and true anomaly.

tan(E/2) = sqrt((1-e)/(1+e))*tan(TA/2)

Since this equation uses half-angle tangents, you don't need to do quadrant checks. The inverse tangent will give you a half-angle in the range [-pi/2, pi/2], so that when you multiply by 2 you get the full range [-pi,pi].

Hope that helps.

 

[jedidia]

Thanks for the formula. It doesn't make that much difference currently, since setting up my program to work with negative angles would take a bit of work now, and I all just got it working yesterday (not the whole program of course, just the particular element that you don't have to manualy adjust your LPe anymore...).

 

[jedidia]

increase the magnitude of the velocity vector by however much dV you are predicting the engine to provide in that time step

I got pretty much everything working that didn't need KOST, now I'm trying to write the calculator for precisely predicting time to eject, deviation of the periapsis etc. For that I need to do exactly what you said above. I familiarized myself a bit with the KOST functions, and everything seems to be pretty clear, just one thing I don't understand yet (since I never worked with state vectors before, and looking at the source code of state vector mfd doesn't help that much):

If I understood the concept of state vectors correctly, it's simply a three-dimensional vector describing position and velocities along the X/Y/Z axis of a three dimensional cartesian coordinate system. Now the first question would be, what's it centered on? (I assume heliocentric?).
the second and tougher question, how would I go about finding which way is prograde? I assume I'd have to convert them first into vectors relative to the body I'm orbiting, then add them up, apply the appropriate DeltaV to it and split it up again. Does that sound about right?

 

[tblaxland]

If I understood the concept of state vectors correctly, it's simply a three-dimensional vector describing position and velocities along the X/Y/Z axis of a three dimensional cartesian coordinate system.

Just to clarify: the state vectors are a set of two 3 dimensional vectors, one for position and one for velocity.

Now the first question would be, what's it centered on? (I assume heliocentric?).

It depends

In State Vector MFD, the state vectors are relative to the current reference body (the reference body is user selectable). The reference body defaults to being the body with the major gravitational influence on the vessel at the time the MFD initialises.

For input into the KOST functions, the state vectors should be relative to the body you are orbiting.

the second and tougher question, how would I go about finding which way is prograde? I assume I'd have to convert them first into vectors relative to the body I'm orbiting, then add them up, apply the appropriate DeltaV to it and split it up again. Does that sound about right?

The prograde direction is parallel to the velocity vector. So, in pseudocode:

kostStateVector state; // state vector
oapiGetRelativeVel(vesselHandle,referenceBodyHandle,state.vel); // get velocity vector
state.vel += deltaVelocityForThisTimeStep * (state.vel / length(state.vel); // add delta-v in the prograde direction

EDIT: The Orbiter API provides the "length" function and the "+=" operator override for you :speakcool:

 

[jedidia]

For input into the KOST functions, the state vectors should be relative to the body you are orbiting.

silly me, I'm giving the function my orbital elements, so it should have been pretty obvious that I'm getting the state vectors back centered on the center of that orbit... That makes things a lot easier.

The prograde direction is parallel to the velocity vector. So, in pseudocode:

this makes things even more easy. I was already writing a function for putting the vectors together and taking them appart again... allthough, if I'm going to try to introduce plane-change in the ejection burn, I'll probably still have to. But I'm not sure yet if I will, let's take this one step at a time. Thanks a lot!

just one question on your bit of pseudo-code, in

(state.vel / length(state.vel)

what does the length stand for?

 

[tblaxland]

just one question on your bit of pseudo-code, in

(state.vel / length(state.vel)

what does the length stand for?

"length" is the length of the vector (also known as magnitude), ie:

sqrt(state.vel.x^2 + state.vel.y^2 + state.vel.z^2)

But happily, the Orbiter API already provides you a function that performs that calculation for you. The function name is "length" and it takes a VECTOR3 as the argument and returns a double.

The Orbiter API also provides a divide ("/") operator overload for the VECTOR3 type, so state.vel / length(state.vel) will return a VECTOR3 that has unit length and points in the direction of the velocity vector.

 

[jedidia]

cool. Thanks!

---------- Post added at 05:07 PM ---------- Previous post was at 09:35 AM ----------

It's really not on topic either, but I just noticed that I have no Idea how to calculate fuel consumption per second (in kg) for a given engine. Neither do I seem to be able to find it on the net...