Determining the orbital parameters

Hi Gents,
I'm trying to write something to go from kepler parameters to a position and
velocity vector.
This article was a great help
http://www.ucalgary.ca/~sneeuw/lectu...423/kepler.pdf

But I can only get orbits that look real (ie: the Earth is in the correct
place) if I change the sign of the rotations, specifically, change the line
r = R3(-Om)R1(-I)R3(-w)q

into

r = R3(Om)R1(I)R3(i)q

and same for velocity. Then all my planets end up in the right 'place'.

When going back from velocity to position, I again need the sign change, but
I'm having unrelated problems with totally circular orbits.
Obtaining a (semi-major azis) and e (eccentricity) from the position and
velocity vector is easy enough, but if e is zero, then I do not know how to
calculate the eccentric anomaly (E).
An equation is given for distance a body at a distance r
Cos E = (a-r)/(a*e)
As you can see, for circular orbits this is undefined.

I realise that in the circular orbit case the eccentric anomaly is equal to
the mean anomaly, but the above paper calculates the mean anomaly from the
eccentric anomaly! Simply obtaining the angle in the q-plane doesn't work
either, because one of the rotations to create the q-plane is created from
the Eccentric anomaly!
But I waffle on...How can I calculate the eccentric anomaly for circular
orbits?

"Makhno" wrote in message ...
Hi Gents,
I'm trying to write something to go from kepler parameters to a position and
velocity vector.
This article was a great help
http://www.ucalgary.ca/~sneeuw/lectu...423/kepler.pdf

But I can only get orbits that look real (ie: the Earth is in the correct
place) if I change the sign of the rotations, specifically, change the line
r = R3(-Om)R1(-I)R3(-w)q

into

r = R3(Om)R1(I)R3(i)q

and same for velocity. Then all my planets end up in the right 'place'.

When you say they end up not in the right place, are
you referring to a screen plot of the position? Make
sure that things aren't getting confused by the
orientation of the screen axes. For example, in BASIC,
the Y-axis runs from the top-left corner to the
bottom-left corner, increasing Y-values are further
down the screen.


When going back from velocity to position, I again need the sign change, but
I'm having unrelated problems with totally circular orbits.
Obtaining a (semi-major azis) and e (eccentricity) from the position and
velocity vector is easy enough, but if e is zero, then I do not know how to
calculate the eccentric anomaly (E).
An equation is given for distance a body at a distance r
Cos E = (a-r)/(a*e)
As you can see, for circular orbits this is undefined.

I realise that in the circular orbit case the eccentric anomaly is equal to
the mean anomaly, but the above paper calculates the mean anomaly from the
eccentric anomaly! Simply obtaining the angle in the q-plane doesn't work
either, because one of the rotations to create the q-plane is created from
the Eccentric anomaly!
But I waffle on...How can I calculate the eccentric anomaly for circular
orbits?

Exactly which parameters are you given to start with
for your circular orbit? The classical orbit elements
can be tricky to work with in specific cases. For example,
for an equatorial orbit Om, the longitude of the ascending
node is undefined. For a circular orbit, the argument of
periapsis and the true anomaly (along with the mean and
eccentric anomaly) are undefined, since there is no
periapsis. In these cases you have to work from the
longitude of the ascending node. Perhaps the argument of
latitude at epoch will be given (angle between ascending
node and the position vector at epoch)?

When you say they end up not in the right place, are
you referring to a screen plot of the position

The correct place in space. I have tables which say where Earth and Mars
should be. I can only get them there by reversing the sign of the rotations.
It probably is an axis problem, though I cannot see what. My axes are very
consistant. But I am wondering if the author of the paper I mentioned
previously made an error.

Exactly which parameters are you given to start with
for your circular orbit?

Semi major axis, eccentricity, inclination, longitude of ascending node,
longitude of perihelion, mean longitude [at epoch].

The classical orbit elements
can be tricky to work with in specific cases. For example,
for an equatorial orbit Om, the longitude of the ascending
node is undefined. For a circular orbit, the argument of
periapsis and the true anomaly (along with the mean and
eccentric anomaly) are undefined, since there is no
periapsis.

So for an equitorial, circular orbit I lose 3 parameters?
How can the true anomaly be undefined? The body must lie at some point on
the orbit at a given time.

In these cases you have to work from the
longitude of the ascending node. Perhaps the argument of
latitude at epoch will be given (angle between ascending
node and the position vector at epoch)?

The argument of latitude at epoch can be calculated (ref the paper), but
depends on Om.

"Makhno" wrote in message ...
When you say they end up not in the right place, are
you referring to a screen plot of the position

The correct place in space. I have tables which say where Earth and Mars
should be. I can only get them there by reversing the sign of the rotations.
It probably is an axis problem, though I cannot see what. My axes are very
consistant. But I am wondering if the author of the paper I mentioned
previously made an error.

Could happen. I prefer to use a rotation matrix
approach to do coordinate transformations. But you
still need the appropriate angles to plug in.


Exactly which parameters are you given to start with
for your circular orbit?

Semi major axis, eccentricity, inclination, longitude of ascending node,
longitude of perihelion, mean longitude [at epoch].

So, for your circular orbit the longitude of perihelion
will be undefined (no perihelion). If it happens to be an
equitorial orbit, you lose the ascending node, too.

The mean longitude at epoch will help. For the circular
orbit it should be the same as the true longitude at epoch,
which is the angle between the Vernal Equinox and the
radius vector at epoch, measured eastwards to the
ascending node, and then in the orbital plane to the
radius vector.


The classical orbit elements
can be tricky to work with in specific cases. For example,
for an equatorial orbit Om, the longitude of the ascending
node is undefined. For a circular orbit, the argument of
periapsis and the true anomaly (along with the mean and
eccentric anomaly) are undefined, since there is no
periapsis.

So for an equitorial, circular orbit I lose 3 parameters?
How can the true anomaly be undefined? The body must lie at some point on
the orbit at a given time.

Right. The problem is the lack of the perihelion; there's
no point of reference for the true anomaly for a circular
orbit without explicity defining one. Often the ascending
node is chosen. For circular equatorial orbits, the
Vernal Equinox (I axis in the IJK frame) is often used.

I think for this particular problem you may want to work
in the IJK frame, and use a few rotations to construct
various working vectors.

A handy routine to implement is the rotation of a point
about an arbitrary line, centered on the origin. See,
for example,

http://www.mines.edu/~gmurray/ArbitraryAxisRotation/ArbitraryAxisRotation.html

let's take your case of a circular orbit where you are
given the semimajor axis (orbit radius, a), inclination (i),
longitude of ascending node (Om), and mean longitude at epoch
(M).

A vector parallel to the node vector follows directly
from Om:

n = cos(Om)I + sin(Om)J (I, J, and K are axis unit vectors)

A vector normal to the plane of the orbit (and hence parallel
to the angular momentum vector) is then obtained by rotating
a unit vector on the K axis by i about vector n; use the
rotation routine mentioned above to rotate the end point of
the unit vector (0,0,1) around the vector n. Call the
resulting vector h (note that it's not really the specific
angular momentum vector, but it is parallel to it).

As I mentioned previously, M is the angle between the Vernal
Equinox (I-axis)and the radius vector at epoch, measured
eastwards to the ascending node, and then in the orbital plane
to the radius vector. If you subtract Om from this, you are
left with the angle in the plane of the orbit from the line of
nodes (vector n) to the radius vector. Call this angle u, the
argument of latitude at epoch.

So a vector in the direction of the radius vector can be had by
rotating n around h by angle u. Normalize it and multiply it by
the semimajor axis and you've got the radius vector in the IJK
frame.

The velocity vector for a circular orbit is perpendicular to the
radius vector, and lies in the plane of the orbit. A vector in
the right direction will be given by h x r. Normalize this and
multiply by the orbital speed, sqrt(mu/r).

I think I got that right. It's been a while.

I prefer to use a rotation matrix
approach to do coordinate transformations. But you
still need the appropriate angles to plug in.

I am using a rotation matrix. Why did you think I was not?

Semi major axis, eccentricity, inclination, longitude of ascending node,
longitude of perihelion, mean longitude [at epoch].

So, for your circular orbit the longitude of perihelion
will be undefined (no perihelion). If it happens to be an
equitorial orbit, you lose the ascending node, too.

Just to be clear, I am going from the position and velocity at a given time,
to the orbital elements.

The mean longitude at epoch will help. For the circular
orbit it should be the same as the true longitude at epoch,
which is the angle between the Vernal Equinox and the
radius vector at epoch, measured eastwards to the
ascending node, and then in the orbital plane to the
radius vector.

....

The problem is the lack of the perihelion; there's
no point of reference for the true anomaly for a circular
orbit without explicity defining one. Often the ascending
node is chosen. For circular equatorial orbits, the
Vernal Equinox (I axis in the IJK frame) is often used.

The ijk frame you mention, is that in the 'q-plane' of the orbit (the orbit
before the rotations that would be applied by I,w and Om) or the 'real' axis
system?

let's take your case of a circular orbit where you are
given the semimajor axis (orbit radius, a), inclination (i),
longitude of ascending node (Om), and mean longitude at epoch
(M).

Sorry, going this way round is easy. It's going back that's the problem.
If I'm given a position and a velocity at a point on a elliptical orbit,
everything's fine. If the orbit turns out to be circular though, then I
can't calculate those undefined parameters!
Presumably, in this case (e=0) I could just set them to zero?

"Makhno" wrote in message ...
I prefer to use a rotation matrix
approach to do coordinate transformations. But you
still need the appropriate angles to plug in.

I am using a rotation matrix. Why did you think I was not?

All I had to go on was your one web page reference showing
a particular example. Sorry, I didn't realize that the
R(x)'s were rotation matrix elements.


Semi major axis, eccentricity, inclination, longitude of ascending node,
longitude of perihelion, mean longitude [at epoch].

So, for your circular orbit the longitude of perihelion
will be undefined (no perihelion). If it happens to be an
equitorial orbit, you lose the ascending node, too.

Just to be clear, I am going from the position and velocity at a given time,
to the orbital elements.

That's not what you said before! You said that you
were trying to find the position and velocity vectors
from a given set of orbital elements (I asked you what
parameters you were given to start with).


The mean longitude at epoch will help. For the circular
orbit it should be the same as the true longitude at epoch,
which is the angle between the Vernal Equinox and the
radius vector at epoch, measured eastwards to the
ascending node, and then in the orbital plane to the
radius vector.

...

The problem is the lack of the perihelion; there's
no point of reference for the true anomaly for a circular
orbit without explicity defining one. Often the ascending
node is chosen. For circular equatorial orbits, the
Vernal Equinox (I axis in the IJK frame) is often used.

The ijk frame you mention, is that in the 'q-plane' of the orbit (the orbit
before the rotations that would be applied by I,w and Om) or the 'real' axis
system?

The IJK frame is the inertial frame defined by the
ecliptic plane as the IJ-plane, and the Vernal Equinox
(First Point of Aries) as the direction of the I-axis.


let's take your case of a circular orbit where you are
given the semimajor axis (orbit radius, a), inclination (i),
longitude of ascending node (Om), and mean longitude at epoch
(M).

Sorry, going this way round is easy. It's going back that's the problem.
If I'm given a position and a velocity at a point on a elliptical orbit,
everything's fine. If the orbit turns out to be circular though, then I
can't calculate those undefined parameters!
Presumably, in this case (e=0) I could just set them to zero?

No, they're actually undefined! It takes a set of six
parameters to completely define an orbit, but they are
not always the same six. In certain cases you must
choose other parameters to replace the undefined ones.

If I may suggest, you might want to get your hands on a
copy of "Fundamentals of Astrodynamics" by Bate, Mueller,
and White. It's a rather inexpensive paperback (you can
get it used for less than ten dollars, new for less than
twelve), but quite excellent. They present a foolproof
procedure for going from position and velocity to
orbital elements, in just a few pages.