Orbital Ellipse from State Vectors (math help)

[HarvesteR]

Hi again,

I'm studying math again :thumbup:
Spent last light reading through Wikipedia, and got reminded of how much I hate math notation

Anyways, what I'm trying to figure out is how to plot an orbital ellipse from initial state vectors...

The trouble I'm having is with some parameters which, for simplicity, I've ignored in my simulation...

For instance, the force of gravity in my sim is a constant -9.81m/s/s, and the primary body (the planet) is completely stationary... Which means the planet has no specified mass...

since I'm applying 1g of force as gravity, I assumed it was safe to use Earth's standard gravitational constant... but I'm not too sure about that...

also, my engine uses Y-up coordinates, so I also assumed it was correct to swap out Z and Y from the original formulae from wikipedia, since it says all it's vectors are in Z-up notation...

Here's what code I have until now... all variables are commented as to what they are (or should be)... the initial state vectors are taken from the static FlightState class that contains the current state of the ship, relative to the center of the planet.

float inc; //inclination
    float ecc; //eccentricity
    float lan; //longitude of ascending node
    float aPe; //argument of periapsis
    float meA; //mean anomaly
    float orP; //orbital period

    float sma; //semi-major axis

    float soe; //specific orbital energy

    Vector3 ram; //specific relative angular momentum vector

    float u = 398600.441f; // standard gravitational constant for earth
    float G = 6.6725985e-11f; // Big G constant

    void Update()
    {
        // let's go
        // I'm assuming that all calculations from Wikipedia are Z-up, so I'm swapping out Z and Y wherever they're present

        Vector3 pos = FlightState.position;
        Vector3 vel = FlightState.velocity;

        ram = Vector3.Cross(pos, vel);

        inc = Mathf.Acos(ram.y / ram.magnitude); //swapped z for y here

        Vector3 eccV = (Mathf.Pow(vel.magnitude, 2f) * pos) / u - (pos / pos.magnitude); //points towards periapsis
        ecc = eccV.magnitude;

        Vector3 lanV = Vector3.Cross(Vector3.up, ram); //lanV points towards the ascending node
        lan = lanV.z >= 0f ? Mathf.Acos(lanV.x / lanV.magnitude) : 2f * Mathf.PI - Mathf.Acos(lanV.x / lanV.magnitude);

        aPe = Mathf.Acos(Vector3.Dot(lanV, eccV) / (lanV.magnitude * eccV.magnitude));

    }

Hopefully you guys can make sense of this... it's probably noobishly wrong... (hence this post)

note that it's also unfinished... I couldn't make sense of the formulae for all the parameters, since my simulation is missing a few, so I'm unsure what to use...

what I intend to get out of this is a full transform to apply to a deformable ellipse... but I'm starting out just by trying to calculate all the needed parameters... later I'll deal with the engine specifics

Many thanks in advance for any help

Cheers

 

[tblaxland]

I think you might find that most of the hard work has already been done for you, here: [ame="http://www.orbithangar.com/searchid.php?ID=3825"]KOST 0.75[/ame] (lots of orbit calcs there); and here: [ame="http://www.orbithangar.com/searchid.php?ID=3609"]Free Orbit MFD 1.2[/ame] (includes stuff that you are interested in like drawing ellipses). The code is generally well commented. If you use it, be sure to read the licences

 

[HarvesteR]

Thanks!! I'll definitely check those out!

The thing is that I can't use any code directly, since I'm not coding for Orbiter... so at best i can read up on the algorithms, and try to translate it into something Unity can understand.

In any case, any help is greatly appreciated :thumbup:

Cheers

 

[Moach]

i believe you might be reinventing the wheel...

in my making of DeltaVee (a miniature orbiter-like sim made for simply "passing" the games-for-cellphones class in college) i did encounter and solve that same problem....

it really isn't as hard as it looks.... but i can't remember how i did it....
so knock yourself out - it's all in spaghettish AS2 (*gasp/cringe*)

but the math should be mostly portable,,,,


just note that a lot of fat has been cut off... a bunch of orbital parameters are simply not needed / ignored due to the highly simplistic nature of the flight model



edit:


found it, here's how i did it:

        vee.y -= grav * flightScale;
	vee.x += wind * flightScale;
	
	
	travel += vee.x * flightScale;
	

	// lets try to make altitude change in function of travel...
	//
	var Rds = alt + rad;
	var orbitD = Math.sqrt((Rds*Rds) + (vee.x*vee.x)) - Rds;
	
	
	vee.y += orbitD;
	alt += (vee.y * flightScale);

it's actually REALLY simple.... just think of the distance travelled as a side of a right triangle, and the altitude (+ planet radius) as the other, then figure the hypotenuse for the next altitude you'd get if you were flying a straight line

finally, subtract your last altitude + planet radius from that and you'll get a vertical velocity adjustment that should guide your ship into an orbital trajectory

simple as that... no muss, no fuss... (tho Kepler tosses a bit in his grave) :thumbup:

 

[HarvesteR]

Hmm interesting...

This shall prove useful for my noob-friendly, orbiter-geek-cringing, autopilot :blink:

I got the trajectory prediction to work by integrating the physics ahead of time... so I managed to plot out an accurate flight path prediction.

this works because all the forces that act on the ship in flight are fed by a function that returns the value of that force at a specified position. Both the ship and the trajectory integrator grab their values from that same function, so it's time-free.

Later on I'll do it for all the other forces too, so I can plot out a path that will account for air drag and winds too

The autopilot will probably work by intercepting the player control feed and doing a few uh... corrections

Cheers

 

[Arrowstar]

I wrote a simple MATLAB function to take state vectors (position and velocity) and turn them into keplerian orbital elements. You can use that information, along with the conic equation, to draw the ellipse. Even though it's in MATLAB, it should be fairly straight forward to follow. Let me know if you have any questions.

function [sma, ecc, inc, longAscNode, ArgPeri, TrueAnom] = getKeplerFromState(rVect,vVect,muCB)
% getKeplerFromState() returns Keplerian orbital elements when provided
% with the state (cartesian position vector, cartesian velocity vector) of
% a spacecraft or celestial body.
%
% INPUTS
% rVect - a 3x1 vector that contains the x,y,z components of the orbiting
% body's current position relative to the central body.  Units: [km]
% vVect - a 3x1 vector that contains the x,y,z components of the orbiting
% body's current velocity vector relative to the central body.  Units:
% [km/sec]
% muCB - the gravitational parameter of the central body.  Units: km^3/s^2
%
%OUTPUTS
% sma - semi-major axis of the orbit.  Units: [km]
% ecc - eccentricity of the orbit.  Units: dimensionless
% inc - inclination angle of the orbit.  Units: radian
% longAscNode - Longitude of ascending node of the orbit.  Units: radian
% ArgPeri - Argument of periapse of the orbit.  Units: radian.
% TrueAnom - Current true anomaly of the spacecraft/body in the orbit.
% Units: radian

r=norm(rVect);
rUnitVect=rVect/r;
v=norm(vVect);

hVect=cross(rVect,vVect);
h=norm(hVect);
hUnitVect=hVect/h;
ThetaUnitVect=cross(hUnitVect,rUnitVect);

Energy=v^2/2 - muCB/r;
sma=-muCB/(2*Energy);

p=h^2/muCB;
ecc=sqrt(-p/sma + 1);

TrueAnom=acos((p/r - 1)/(ecc));
if(dot(rVect,vVect)<0)
    TrueAnom=-TrueAnom;
end

inc=acos(hUnitVect(3));

longAscNode_1=AngleZero2Pi(asin(hUnitVect(1)/sin(inc)));
longAscNode_2=AngleZero2Pi(pi-asin(hUnitVect(1)/sin(inc)));
longAscNode_3=AngleZero2Pi(acos(-hUnitVect(2)/sin(inc)));
longAscNode_4=AngleZero2Pi(-acos(-hUnitVect(2)/sin(inc)));

longAscNodeSet1=round(1000*[longAscNode_1,longAscNode_2])/1000;
longAscNodeSet2=round(1000*[longAscNode_3,longAscNode_4])/1000;

[val,ia,ib]=intersect(longAscNodeSet1,longAscNodeSet2);
longAscNode=longAscNodeSet1(ia);

Theta_1=AngleZero2Pi(asin(rUnitVect(3)/sin(inc)));
Theta_2=AngleZero2Pi(pi-asin(rUnitVect(3)/sin(inc)));
Theta_3=AngleZero2Pi(acos(ThetaUnitVect(3)/sin(inc)));
Theta_4=AngleZero2Pi(-acos(ThetaUnitVect(3)/sin(inc)));

ThetaSet1=round(1000*[Theta_1,Theta_2])/1000;
ThetaSet2=round(1000*[Theta_3,Theta_4])/1000;

[val,ia,ib]=intersect(ThetaSet1,ThetaSet2);
Theta=ThetaSet1(ia);

ArgPeri=Theta-TrueAnom;

 

[HarvesteR]

Oh, wow, gotta love this forum!!

Thanks a million man!! this is perfect!! :tiphat:

So very well commented too!! I stripped down an orbital parameter calculator I found online, but the code didn't have a single comment... couldn't make sense of what each variable was, or what units the guy used...

This is exactly what I needed!

Just one question, are the vectors Y-up or Z-up? (or does it matter?) Wikipedia says most orbital maths are done in Z-up... but my engine is Y-up, so I might have to convert... no problem!

Thanks again!!

Cheers

 

[Arrowstar]

The x,y,z are all in the J2000 frame, so I guess you can think of 'z' is up. Since we're in space, though, 'up' is somewhat arbitrary.

 

[HarvesteR]

Yeah, I meant in relation to the reference frame... In my sim that's just the absolute XYZ axes from the engine

Cheers