Calculating E or M from a State Vector

[kepler]

Hi,

Does anyone knows how to calculate the eccentric anomaly E, or the mean anomaly M, given a state vector?

Regards,

Kepler

 

[kepler]

Work, work...

Hi again,

After picking up a paper and a pencil I was able to calculate what I wanted from scratch. Now, I have a doubt: I've created a routine that gives the orbital elements of a body from the state vectors ( position and velocity in 3 axis ). I tryied to use the system to the Moon and Earth; I get always wrong results.

Aren't the laws of physics the same for the interaction of two bodies? - I only corrected the GM parameter.

Any ideas?

Kepler

 

[tony873004]

mu = 398600.44
mu = mu * Mass 'mass of the central body of your system expressed in Earth masses (ie Earth = 1, sun = 332945.65)
rmu1 = 1 / mu

hi = Y * zd - z * yd ' x,y,z are positions, xd, yz, zd are velocities
hj = z * xd - X * zd
hk = X * yd - Y * xd
svH = Sqr(hi * hi + hj * hj + hk * hk)
r = Sqr(X * X + Y * Y + z * z)
v2 = (xd * xd + yd * yd + zd * zd) - mu / r
rv = X * xd + Y * yd + z * zd
ei = rmu1 * (v2 * X - rv * xd)
ej = rmu1 * (v2 * Y - rv * yd)
ek = rmu1 * (v2 * z - rv * zd)
ec = Sqr(ei * ei + ej * ej + ek * ek) 'eccencitrity
vn = Sqr(hj * hj + hi * hi)
apg = Acos((-ei * hj + ej * hi) / (vn * ec))
apg = apg * 180# / Pi ' perapsis in degrees instead of radians
ai = Acos(Abs(hk / svH)) 'inclination
ai = ai * 180# / Pi ' convert to degrees
anl = Acos(-hj / vn)
If (hi <= 0#) Then anl = Pi - anl
anl = anl * 180# / Pi ' assending node in degrees
sma = svH * svH * rmu1 / (1 - ec * ec) ' semi-major axis

 

[Jenab]

Eccentric and mean anomalies from a state vector.

Does anyone knows how to calculate the eccentric anomaly E, or the mean anomaly M, given a state vector?

Quantities in the procedure below will default to MKS units and radians for the angles. Exceptions will be noted.

Initial data.
Heliocentric position in ecliptic coordinates: x, y, z.
Sun-relative velocity in ecliptic coordinates: Vx, Vy, Vz.
An epoch (i.e., a date/time of validity) for the position & velocity: t.

(Note: The epoch will be entered in decimal days.)

Handy constants and conversion factors.

AU = 1.49597870691E+11 meters

GMsun = 1.32712440018E+20 m^3 sec^-2

DPY = 365.256898326 days

The semimajor axis (a).

r = [ x^2 + y^2 + z^2 ]^0.5

V = [ Vx^2 + Vy^2 + Vz^2 ]^0.5

a = 1 / { 2/r - V^2 / GMsun }

Angular momentum per unit mass (h).

hx = Y Vz - Z Vy

hy = Z Vx - X Vz

hz = X Vy - Y Vx

h = [ hx^2 + hy^2 + hz^2 ]^0.5

The eccentricity (e).

p = h^2 / GMsun

q = x Vx + y Vy + z Vz

There are two ways to get the eccentricity. The first is...

e = [ 1 - p/a ]^0.5

The second way to get the eccentricity is...

ex = 1 - r/a

ey = q / { a GMsun }^0.5

e = [ ex^2 + ey^2 ]^0.5

The two values for (e) should check each other.

The eccentric anomaly (u).

u = arctan2(ey,ex)

Note: The arctan2 function is the two-argument version of the arctangent, which adjusts the result for the actual quadrant of ex and ey.

The two-argument arctangent function: arctan2(ey,ex).

If ex is NOT zero, then

BEGIN

| u = arctan(ey/ex)

| If (ex < 0), then increase u by pi radians.

| If (ex > 0) and (ey < 0), then increase u by 2 pi radians.

END ELSE (if ex = 0) BEGIN

| If (ey < 0), then u = -pi / 2

| If (ey = 0), then u = 0

| If (ey > 0), then u = +pi / 2

END

arctan2 = u

The mean anomaly (M).

M = u - e sin u

The inclination (i).

i = pi/2 - arcsin (hz/h)

or, in other words...

i = arctan { [ ( h^2 / hz^2 ) - 1 ]^0.5 }

However,

If (hz < 0), then i (revised) = pi radians - i

The longitude of the ascending node (L).

L = arctan2(hx , -hy)

The true anomaly (F).

Fx = h^2 / (r GMsun) - 1

Fy = h q / (r GMsun)

F = arctan2(Fy , Fx)

The argument of the perihelion (w).

xx = ( x cos L + y sin L ) / r

IF (i = 0), THEN
yy = ( y cos L - x sin L ) / r
ELSE (if i is not zero)
yy = z / (r sin i)

w = arctan2(yy , xx) - F

Adjust (w), if necessary, to the interval [0, 2 pi) by adding or subtracting whole multiples of 2 pi radians.

The orbital period (P).

P = DPY (a/AU)^1.5

(P results in days.)

The mean motion (m).

m = 2 pi / P

(m results in radians per day.)

The time of perihelion passage (T).

T = t - M/m

(T results in decimal days.)

Jerry Abbott

 

[kepler]

The vectors

Hi Tony,

I've tested the routine that you gave here, and the results aren't correct for the Earth-Moon system by some misterious cause. I tried with AU as unit of distance, and with meters. It gives wrong results both ways; but just to be sure, here are the value of the vectores for the moon at 1/1/2000, 12h00, in AU units:

x = -0.0019490056;
y = -0.0018384464;
z = 0.0002424016;
vx = 0.0003717380;
vy = -0.0004221159;
vz = -0.0000066657;

Maybe I wrote something wrong in the routine.

Regards, and thanks,

Kepler

 

[Jenab]

I've tested the routine that you [Tony] gave here, and the results aren't correct for the Earth-Moon system by some misterious cause.

Oops. The procedure I gave above was for orbits around the sun. You would need to change GMsun to GMearth for it to work with low Earth orbits. I used MKS units for most of the quantities, so most of it should work if you change GM. The period calculation still won't work, though, until you change it to some non-canonical MKS form. But you don't need to go that far down the procedure just to get the mean and eccentric anomalies.

GMearth = 3.9885E+14 m^3 sec^-2

Jerry Abbott

 

[tony873004]

Hi Tony,

I've tested the routine that you gave here, and the results aren't correct for the Earth-Moon system by some misterious cause. I tried with AU as unit of distance, and with meters. It gives wrong results both ways; but just to be sure, here are the value of the vectores for the moon at 1/1/2000, 12h00, in

Actually, it wants kilometers and km/s. I inputted the following values for the moon's earth-relative position on Jan 1, 2002 12:00:00:

X = -226545.429682968: y = 286349.964075857: z = 22167.1721237376
xd = -0.83215986700533: yd = -0.67880264763759: zd = 0.073458217887696

and it returned the following:

sma= 390528.587811472
ecc= 6.42167877853417E-02
perapsis = 51.5661090865652
inclination = 5.26463171719322
ascending node 87.1361489019904

Which seems in the right ball park anyway.

The above input was generated by JPL Horizons system for the Moon for Jan 1, 2002 at 12:00. Asking for elements instead of vectors gives the following (angles in radians):

JDCT Epoch Julian Date, Coordinate Time
EC Eccentricity, e
QR Periapsis distance, q (km)
IN Inclination w.r.t xy-plane, i (degrees)
OM Longitude of Ascending Node, OMEGA, (degrees)
W Argument of Perifocus, w (degrees)
Tp Time of periapsis (Julian day number)
N Mean motion, n (degrees/sec)
MA Mean anomaly, M (degrees)
TA True anomaly, nu (degrees)
A Semi-major axis, a (km)
AD Apoapsis distance (km)
PR Orbital period (sec)

2452276.000000000 = A.D. 2002-Jan-01 12:00:00.0000 (CT)
EC= 5.152445531449781E-02 QR= 3.653672052491503E+05 IN= 5.264631717193200E+00
OM= 8.713614890199040E+01 W = 5.396714814738498E+01 Tp= 2452276.865982342046
N = 1.522267598936786E-04 MA= 3.486102607231271E+02 TA= 3.473660261921209E+02
A = 3.852152090756327E+05 AD= 4.050632129021150E+05 PR= 2.364893007323014E+06

I purposely left out my code for mean anomoly since my routine gets it wrong. I'm anxious to see if Jeneb's code works for this. I tried his code and got a wrong value for inclination, but I haven't checked my work for typos yet. Something to do after work I guess...

Also, I'm not sure how much difference it makes, but I believe the small differences between Horizons output and my code's output is because my code treats the object orbiting Earth as a massless test particle. I'll try this later, but I believe entering the Mass of Earth + Mass of Moon for the variable mass will give me a more correct answer for SMA and Ecc. Not sure about the others.

 

[kepler]

I'm a little bit confused...

Hi Tony,

it seems the problem might be solved according to what you describe; but I'm a little bit confused: which set of routines did you use?(there are two here in this thread) And what value did you use for GM (Mu) and rmu1?

Regards,

Kepler

 

[tony873004]

Hi, Kepler, I used the first set of routines, the one I posted. The first 2 lines of this code are:

mu = 398600.44
mu = mu * Mass 'mass of the central body of your system expressed in Earth masses (ie Earth = 1, sun = 332945.65)

I used 1 for mass to represent Earth's mass. So, 398600.44 * 1 = 398600.44 is my value for mu when computing the Moon's orbit around Earth.

I've been having good luck with Jenab's method. I might even prefer it to the one I posted, once I debug my slight translation a little more... silly stuff like making sure that my computer knows what Pi is before using it in code.

 

[kepler]

...oopsss...

Dear Tony,

the problem was that I had initial data incorrect - the routine does work. Complemented with mine, everything's ok - but it doesn't corresponds to the JPL data. Now I have a doubt ( the final one I think ): in the system that you computed, with that Mu value ( 398600.44 ), what are the units for velocity if for the position were to be AU?

And just another thing: where, in JPL Horizons, do we choose to output position/velocity vector and orbital elements?

Do you have any idea about this?

Regards,

Kepler

 

[tony873004]

Dear Tony,

the problem was that I had initial data incorrect - the routine does work. Complemented with mine, everything's ok. Now I have a doubt ( the final one I think ): in the system that you computed, with that Mu value ( 398600.44 ), what are the units for velocity if for the position were to be AU?

And just another thing: where, in JPL Horizons, do we choose to output position/velocity vector and orbital elements?

Do you have any idea about this?

Regards,

Kepler

If you want to input your distances in AU, you can do one of two things. The easiest is to just add a few lines of code like:
AU= 149000000 ' look up a more accurate value for au
x=x*AU
y=y*AU
z=z*AU

Or you could modify mu to work with your desired units. Try this: First, use the existing mu to compute SMA. Then do a hand calculation to convert this SMA into your desired units. Then change x,y,z and xd, yd, zd to your desired units. Then keep altering mu until you zero-in on the correct SMA for the new units. It'll probably be something like current mu * au or / au. Try that first.

Horizons works like this:

Create an e-mail to: horizons@ssd.jpl.nasa.gov
Put the word job in the subject line

copy and paste the following text. Modify to suit your needs:

!$$SOF
EMAIL_ADDR=''
START_TIME = '2005-Jan-01 00:00:00'
STOP_TIME = '2005-Jan-01 00:00:01'
TABLE_TYPE = 'Elements'
REF_PLANE = 'Ecliptic'
CENTER = '@0399'
COMMAND='301'
!$$EOF

Note that the start time and stop time are 1 second apart
Table type = "Elements" or "Vectors" I think that's the answer to your question.
center = "@" + ID number of object that is the center of the system
command= ID number of the object that is orbiting the center object.

The numbering scheme works like this:
010 Sun
199 Mercury
299 Venus
399 Earth
...
999 Pluto
301 Earth's Moon
401 Mars' first moon
402 Mars' second moon
501 Jupiter's first moon (Io)
...etc...


Then they send you back something that looks like this for elements:
$$SOE
2452276.000000000 = A.D. 2002-Jan-01 12:00:00.0000 (CT)
EC= 5.152445531449781E-02 QR= 3.653672052491503E+05 IN= 5.264631717193200E+00
OM= 8.713614890199040E+01 W = 5.396714814738498E+01 Tp= 2452276.865982342046
N = 1.522267598936786E-04 MA= 3.486102607231271E+02 TA= 3.473660261921209E+02
A = 3.852152090756327E+05 AD= 4.050632129021150E+05 PR= 2.364893007323014E+06
$$EOE

and very similar if you requested vectors. The rest of their return e-mail will explain the abbreviations to you.

If you send them a blank e-mail with the word batch-long in the subject line, they e-mail you instructions on how to use Horizons.
Also, see this web site: http://ssd.jpl.nasa.gov/

And here's one more website you might find interesting:
http://www.geocities.com/SiliconValley/2902/orbit.htm

 

[BobG]

Do you know how to find your True Anomaly from a state vector?

The Eccentric Anomaly treats all orbits as circular with a radius equal to the semi-major axis and looks at the point on the circle directly above the satellite's location on the ellipse. The second difference is that it's looking at that point from the geometric center of the ellipse instead of the true focus.

There's an equation for that, but you can just look at the ellipse and figure it out. The difference between the geometric center and the true focus is the linear eccentricity (ae or semi-major axis times the eccentricity). You add the cosine of true anomaly times the radius (actual radius) to find the point on the semi-major axis directly below the point you want to be pointing at. You'll have the hypotenuse of a triangle - a line pointing from the eccentric anomaly towards a point on a fictitous circle directly above the satellite's location (key point is that it's equal in length to the semi-major axis) - and the adjacent side of a triangle (the linear eccentricity plus r times cosine of true anomaly). You can find the angle of eccentric anomaly.

Mean anomaly is just E - e sin E.

 

[Jenab]

The eccentric anomaly is hard to describe in words. I just now drew and uploaded a picture.

http://www.jabpage.org/images/anomaly.jpg

Jerry Abbott

 

[tony873004]

The eccentric anomaly is hard to describe in words. I just now drew and uploaded a picture.

http://www.jabpage.org/images/anomaly.jpg

Jerry Abbott

Nice drawing.
What you described as true anomoly, I always thought was mean anomoly. Any chance of putting that in your drawing too?

 

[Jenab]

Nice drawing.
What you described as true anomoly, I always thought was mean anomoly. Any chance of putting that in your drawing too?

The true anomaly is the angle, subtended at the sun, measured from the perihelion, in the plane of the orbit, in the direction of motion, to where the planet truly is.

The mean anomaly is also an angle, subtended at the sun, measured from the perihelion, in the planet of the orbit, in the direction of motion...but there's a difference in where the angle is measured to.

Imagine that the circle circumscribing the ellipse in that drawing were moved to the right (by a distance ae) so that instead of being centered at the center of the ellipse, it was centered on the sun. Now imagine that a "phantom planet" orbits the sun in that circular orbit with the constant angular speed appropriate for a circular orbit of that radius.

The period of the phantom planet is equal to the period of the real planet because the radius of the phantom planet's orbit is equal to the semimajor axis of the real planet's orbit.

But whereas the phantom planet moves around the sun with a constant angular speed, the real planet moves around the sun faster near perihelion and slower near aphelion.

The mean anomaly is measured from the perihelion to where the phantom planet is. Thus the mean anomaly has a constant rate of increase with time. The true and eccentric anomalies increase at variable rates with time.

The true, mean, and eccentric anomalies all have the same values at perihelion (zero) and at aphelion (pi radians). But everywhere else they are all different.

Jerry Abbott

 

[CharlesP]

I want to understand this

Some questions please:
mu = 398600.44
Is this GM for Earth?
mu = mu * Mass 'mass of the central body of your system expressed in Earth masses (ie Earth = 1, sun = 332945.65)
rmu1 = 1 / mu
Is h the orbital angular momentum per unit mass vector components?
hi = Y * zd - z * yd ' x,y,z are positions, xd, ->yd, zd are velocities (did I fix a typo?)
hj = z * xd - X * zd
hk = X * yd - Y * xd
svH = Sqr(hi * hi + hj * hj + hk * hk)
r = Sqr(X * X + Y * Y + z * z)
Is this the total energy per unit mass?
v2 = (xd * xd + yd * yd + zd * zd) - mu / r
This looks like r dot v:
rv = X * xd + Y * yd + z * zd
I have no idea what is going on here or why. Help!
ei = rmu1 * (v2 * X - rv * xd)
ej = rmu1 * (v2 * Y - rv * yd)
ek = rmu1 * (v2 * z - rv * zd)
ec = Sqr(ei * ei + ej * ej + ek * ek) 'eccencitrity
vn = Sqr(hj * hj + hi * hi)
apg = Acos((-ei * hj + ej * hi) / (vn * ec))
apg = apg * 180# / Pi ' perapsis in degrees instead of radians
ai = Acos(Abs(hk / svH)) 'inclination
ai = ai * 180# / Pi ' convert to degrees
anl = Acos(-hj / vn)
If (hi <= 0#) Then anl = Pi - anl
anl = anl * 180# / Pi ' assending node in degrees
sma = svH * svH * rmu1 / (1 - ec * ec) ' semi-major axis

 

[tony873004]

I guess I should attempt to answer your questions since that code is from my post. It is basically the code I use for my program Gravity Simulator ( www.gravitysimulator.com ).

I pretty much hijacked, modified, and translated into Visual Basic, that code from the javascript "view source" of the following web page:

http://orbitsimulator.com/sheela/vec2ele.htm

So although I'll be a little useful in answering your questions, I won't be complete.

Is this GM for Earth?
This is one of the parts I modified.

Yes. G is gravitational constant and M is mass expressed in Earth masses. G is usually memorized to be 6.672 x 10^-11 N m^2 kg^-2.
But actually G comes in many different flavors depending on your input units and your desired output units. Consider the following example:

The formula for figuring out acceleration due to gravity at a given distance from a mass is a = GM/d^2

If you'd like your input units to be expressed in Earth Masses for mass, Earth Radii for radius, and your output unit to be meter / second^2, then G does not equal 6.672 x 10^-11 since we're diverging from the MKS system of units. In this case G = 9.81.

9.81 * 1 / 1^2 = 9.81 meter / second ^2 which is the value for acceleration at Earth's surface.

Try this for the Moon. G = 9.81, mass = ~0.0123, d = ~0.25 (off the top of my head, so this is approximate). 9.81 * 0.0123 / 0.25^2 = 1.93 meters / second^2 (1/5 - 1/6 Earth gravity), approximately correct.

398600.44 is the number that the code needs to use for G for the formulas to work. And it assumes that mass will be expressed in Earth Masses, and that D will be expressed in kilometers, which is why it is different from what a physics student would memorize for the value of G.

Is h the orbital angular momentum per unit mass vector components?
hmmm... I used to know what h was, and that doesn't sound familiar. But I could be wrong.

(did I fix a typo?)
Yes you did!

Is this the total energy per unit mass?
Again, doesn't sound familiar, but remember, I hijacked this code.

I have no idea what is going on here or why. Help!
Me either. Jeneb or BobG might. They're good with this stuff.

And, I don't trust the answers I get for the anomolies. I think one of them is wrong, but I don't remember which.

Also, Sheila Belur, author of the webpage I stole the code from is very good at answering her e-mail. She's a very nice person who would probably be willing to help you understand. I host her webpages for her on my orbitsimulator.com domain now.

 

[BobG]

Some questions please:
mu = 398600.44
Is this GM for Earth?
mu = mu * Mass 'mass of the central body of your system expressed in Earth masses (ie Earth = 1, sun = 332945.65)
rmu1 = 1 / mu
Is h the orbital angular momentum per unit mass vector components?
hi = Y * zd - z * yd ' x,y,z are positions, xd, ->yd, zd are velocities (did I fix a typo?)
hj = z * xd - X * zd
hk = X * yd - Y * xd
svH = Sqr(hi * hi + hj * hj + hk * hk)
r = Sqr(X * X + Y * Y + z * z)
Is this the total energy per unit mass?
v2 = (xd * xd + yd * yd + zd * zd) - mu / r
This looks like r dot v:
rv = X * xd + Y * yd + z * zd
I have no idea what is going on here or why. Help!
ei = rmu1 * (v2 * X - rv * xd)
ej = rmu1 * (v2 * Y - rv * yd)
ek = rmu1 * (v2 * z - rv * zd)
ec = Sqr(ei * ei + ej * ej + ek * ek) 'eccencitrity
vn = Sqr(hj * hj + hi * hi)
apg = Acos((-ei * hj + ej * hi) / (vn * ec))
apg = apg * 180# / Pi ' perapsis in degrees instead of radians
ai = Acos(Abs(hk / svH)) 'inclination
ai = ai * 180# / Pi ' convert to degrees
anl = Acos(-hj / vn)
If (hi <= 0#) Then anl = Pi - anl
anl = anl * 180# / Pi ' assending node in degrees
sma = svH * svH * rmu1 / (1 - ec * ec) ' semi-major axis

[tex]\mu[/tex] is the geocentric gravitational constant. It's equal to the universal gravitational constant times the Earth's mass. For any other body, you can compare the ratio between the other body's mass and the Earth's mass.

(Specific) angular momentum (per unit of mass) (h) is the cross product of the position vector and the velocity vector. So, your equations for hi, hj, and hk are correct. (A subroutine that does cross products and dot products would come in pretty handy, since just about all your elements are obtained via vector operations).

The rest is to find perigee and eccentricity. Your eccentricity vector is a confusing process. First, you compare centrifugal force to centripedal force by finding the Laplace (Runge Lenz) vector.

[tex]\dot {\vec {r}} \times \vec {h} - \frac{\mu}{r}\vec {r}[/tex]

If you divide [tex]\mu[/tex] out of the result, you get your eccentricity vector in a more usable magnitude (It's norm equals the eccentricity of the orbit). It points towards perigee. Measuring the angle between your eccentricity vector and the position vector gets you your true anomaly.

Your equations try to save some work, but don't give as much insight as to how you got there as the actual vectors do.

Is this the total energy per unit mass?
v2 = (xd * xd + yd * yd + zd * zd) - mu / r

No. Your specific energy (per unit of mass) is the sum of your kinetic energy and potential energy with the mass divided out.

[tex]E = \frac {v^2}{2} - \frac{\mu}{r}[/tex]

(And this is important, because the semi-major axis is equal to:

[tex]a = - \frac{\mu}{2E}[/tex]

The equations listed work, but, like I said, it's hard to really picture what's happening with them (for me, anyway, but I like working with vectors).

Two additional notes about the angular momentum vector.

The vector serves as the orbital pole. It's perpendicular to the orbit, so you can use it to find the angle between the orbital plane and the equatorial plane. (Actually, you use a dot product to find the angle between the angular momentum vector and the Z-axis of your geocentric coordinate system).

It's magnitude is equal to the area of a parallelogram with the radius as one side and the velocity as one side. If you divide the magnitude by two, you have the area that your orbit sweeps out per second. Recalling Kepler's second law, this is a pretty key value. If you calculate the total area of your ellipse and divide by the amount of area swept out per second, you're left with the number of seconds required to complete one orbit (Kepler's third law).

 

[CharlesP]

"...
hj = z * xd - X * zd
hk = X * yd - Y * xd
svH = Sqr(hi * hi + hj * hj + hk * hk)
r = Sqr(X * X + Y * Y + z * z)
v2 = (xd * xd + yd * yd + zd * zd) - mu / r
rv = X * xd + Y * yd + z * zd
ei = rmu1 * (v2 * X - rv * xd)
ej = rmu1 * (v2 * Y - rv * yd)
ek = rmu1 * (v2 * z - rv * zd)
ec = Sqr(ei * ei + ej * ej + ek * ek) 'eccentricity
Got this far ok.

vn = Sqr(hj * hj + hi * hi) ' component in XY plane
apg = Acos((-ei * hj + ej * hi) / (vn * ec))
apg = apg * 180# / Pi ' perapsis in degrees instead of radians
argument of perigee? I can't see these vectors working in the formula right even though I know what they mean. This is the only part that I can't get.
I'll keep trying.

ai = Acos(Abs(hk / svH)) 'inclination
ai = ai * 180# / Pi ' convert to degrees ' I think this is the inclination of the orbital plane to the XY plane.

anl = Acos(-hj / vn)
If (hi <= 0#) Then anl = Pi - anl
anl = anl * 180# / Pi ' ascending node in degrees This looks right.

sma = svH * svH * rmu1 / (1 - ec * ec) ' semi-major axis Got this one ok.