Calculation of coordinates from JPL Small-body database parameters?

Question

I would like to write a script that can calculate the distance between two asteroids at some moment in time, given some object in the JPL database of asteroids and comets, so available parameters are restricted to the menu they provide there. Just getting consistent spatial coordinates for objects (in any coordinate system, I don't care about the definition of the system) at a time value would be sufficient to do this.

Trivial conversions are not going to bother me like going between spherical and Cartesian coordinates, or converting date/month/year to floating point. I have more trouble making sense of what those parameters are, and what to do with them.

In terms of degrees of freedom, I ultimately expect 3 functions that take one input parameter (time), and give a spatial coordinate (x,y,z). Given this problem specification, I think need 5 free parameters. I say this in order to cover the extent along each coordinate, the time offset for that period, and the angular orientation of the orbital ellipse. Orbital period is redundant with the spatial data since I know GM of the sun, so I don't need another one for that. The JPL database has a great number of parameters with little obvious meaning. My question is about usage of these. Which ones do I need, and how should I use them?

Here are some of the most promising ones, I will need to sort them into categories. Now, I understand that JED is probably a time unit (not hard). I understand that epoch is probably a reference time and might be needed but won't count toward my 5 expected free variables.

time related

• tp, time of perihelion passage, JED
• epoch, epoch of osculation, JED (not counting)

traditional orbital parameters

• a, semi-major axis, AU
• e, eccentricity

spatial offset concepts

• i, inclination, deg
• peri, argument of perihelion, deg
• node, longitude of the ascending node, deg

I know that this question sounds like Wikipedia/textbook kind of stuff that I could easily look up. So sure, let's go to Wikipedia for their method of calculating coordinates of an orbiting body. They rely on true anomaly (or mean anomaly) for setting the temporal start point. JPL does seem to have mean anomaly, but I must be missing something for that method. That parameter deals with the location of the object at some point in time, and the parameter (variable M) doesn't give any kind of reference point to use. In fact, my idea that epoch is a reference time for other parameters is entirely speculation on my part. That parameter (and only that parameter) is the same for all objects in the database, so I thought that interpretation made sense, but I could be wrong.

I suspect that some parameters are redundant, and that either the time of perihelion passage or the Mean anomaly could be used in its own unique calculation sequence. I don't care which approach I use, I just need one that works.

What is the easiest way for me to go about getting a cookbook to start writing this script to get coordinates of objects in the database?

Thanks in advance for entertaining such a detail-oriented question.

---

Let me take a little more space (in an already crowded question) to respond to a comment.

you are asking us to write a chapter of a book. That's a bit much for a stackexchange Q&A type of question.

That textbook would have tremendously short chapters. But I recognize the hazard of this question looking like "explain orbital mechanics to me!" because that's not what I have in mind. So instead of trying to say this in words, let me just say it in equations.

$$ R = a \frac{ 1 - e ^ 2 }{ 1 + e \cos{ TA } } \\ X = R \left( \cos{ N} \cos{ (TA + w ) } - \sin{N} \sin{ ( TA+w) } \right) \cos{ i} \\ Y = R \left( \sin{ N } \cos{ (TA+w) } + \cos{N} \sin{ (TA+w) } \right) \cos{i} \\ Z = R \sin{ (TA+w) } \sin{i}$$

I just found this from physics forms, and it almost answers my question, aside from a few processing issues... if I interpreted it correctly. The source didn't have consistent parenthesis, but the above form fits my mathematical picture of the problem. Getting TA shouldn't be hard either, considering that all I need is a scriptable solution. But I still haven't seen a form for TA that could obviously be obtained with the current time.

That should be missing just one additional expression or a few other similar details. I wish I could just get some code from the SPICE toolkit that does the same thing, but its source is probably 10,000s of lines of code, doing all kinds of irrelevant things. If the above equations are right, this should really only require just 1 or 2 more small, concise, things.

Answer 1

Use this page to generate SPICE kernels for the bodies of interest, and then use SPICE routines to calculate whatever you like to your heart's content.

---

Aarrr. I don't be knowin why ye be reinventing th' wheel.

Alrighty ya scurvy bilge rat, here it be:

$$x=a\left(\cos\tau-e\right)$$ $$y=a\sqrt{1-e^2}\sin\tau$$ $$z=0$$

That be givin it t' ye in th' plane. Then ya be rotatin th' plane with that land lubber Euler's transform in $\Omega$, $i$, $\omega$ (longitude of ascending node, inclination, and argument of periapsis), by multiplyin his matrix by yer vector up there:

$$\left( \begin{array}{ccc} \cos\omega \cos\Omega-\cos i \sin\omega \sin\Omega & -\cos\Omega \sin\omega -\cos i \cos\omega \sin\Omega & \sin i \sin\Omega \\ \cos i \cos\Omega \sin\omega+\cos\omega\sin\Omega & \cos i \cos\omega\cos \Omega-\sin\omega\sin\Omega & -\cos\Omega \sin i \\ \sin i \sin\omega & \cos\omega \sin i & \cos i \\ \end{array} \right)$$

Ye also be needin th' time it be ($\mu$ being th' $GM$ of the Sun):

$$t=\sqrt{a^3\over\mu}\left(\tau-e\sin\tau\right)$$

Yer $\tau=0$ and $t=0$ being yer time a periapsis. Ye be seein that ev'ry $2\pi$ in $\tau$, it be one a yer spins about th' Sun.

If yer scurvy rock be gettin close to a planet, then ye be wastin yer time, as th' orbit be changin on ye.

Answer 2

!$$SOF
    COMMAND = 'DES=2010 VQ'
    CENTER =  '500@0'
    MAKE_EPHEM = 'YES'
    TABLE_TYPE = 'VECTORS'
    START_TIME = '2017-04-01 00:00:00'
    STOP_TIME = '2017-04-01 00:00:01'
    STEP_SIZE = '1'
    OUT_UNITS =  'AU-D'
    VECT_TABLE =  '3'
    REF_PLANE =  'ECLIPTIC'
    REF_SYSTEM =  'J2000'
    VECT_CORR =  'NONE'
    VEC_LABELS =  'NO'
    CSV_FORMAT =  'YES'
    OBJ_DATA =  'YES'
    !$$EOF

returns heliocentic coordinates:

$$SOE
2457844.500000000, A.D. 2017-Apr-01 00:00:00.0000,  7.539123579055962E-01,  7.098001692200981E-01,  1.408370538460704E-03, -1.019002972030220E-02,  1.127277705358026E-02,  9.280758228718769E-05,  5.980387083611182E-03,  1.035471924841630E+00,  3.082270281389151E-04,
2457844.500011574, A.D. 2017-Apr-01 00:00:01.0000,  7.539122399654240E-01,  7.098002996920419E-01,  1.408371612622512E-03, -1.019003207292571E-02,  1.127277484368998E-02,  9.280757740530809E-05,  5.980387104214990E-03,  1.035471928409068E+00,  3.082263804494916E-04,
$$EOE

Here XYZ coordinates is like:

JDTDB,            Calendar Date (TDB),                      X,                      Y,                      Z,                     VX,                     VY,                     VZ,                     LT,                     RG,                     RR,

The other asteroid:

COMMAND = 'DES=2008 LD'
CENTER =  '500@0'
MAKE_EPHEM = 'YES'
TABLE_TYPE = 'VECTORS'
START_TIME = '2017-04-01 00:00:00'
STOP_TIME = '2017-04-01 00:00:01'
STEP_SIZE = '1'
OUT_UNITS =  'AU-D'
VECT_TABLE =  '3'
REF_PLANE =  'ECLIPTIC'
REF_SYSTEM =  'J2000'
VECT_CORR =  'NONE'
VEC_LABELS =  'NO'
CSV_FORMAT =  'YES'
OBJ_DATA =  'YES'
!$$EOF

....
then do sqrt(dx^2+dy^2+dz^2)