Precessing ellipses

[chris]

I'm trying to model the orbits of Phobos and Deimos as precessing ellipses using the parameters on this page:

http://ssd.jpl.nasa.gov/sat_elem.html#legend

I'm completely frustrated . . . I can't understand why this code (written for clarity, not performance) doesn't work:

Code: Select all

    Point3d computePosition(double jd) const
    {
        double epoch = 2433283.0 - 0.5; // 00:00 1 Jan 1950
        double a     = 23460.0;
        double e     = 0.0002;
        double w0    = 290.496;
        double M0    = 296.230;
        double i     = 1.793;
        double node0 = 339.600;
        double n     = 285.1618919;
        double Pw    = 26.892;
        double Pnode = 54.536;

        double refplane_RA  = 316.700;
        double refplane_Dec =  53.564;
        double marspole_RA  = 317.681;
        double marspole_Dec =  52.886;

        double t = jd - epoch;
        double T = t / 365.25;

        double M    = degToRad(M0 + n * t);
        double w    = degToRad(w0 + 360.0 * T / Pw);
        double node = degToRad(node0 + 360.0 * T / Pnode);

        // Solve Kepler's equation--for a low eccentricity orbit, just a few
        // iterations is enough.
        double E = M;
        for (int k = 0; k < 3; k++)
            E = M + e * sin(E);

        // Compute the position in the orbital plane (y = 0)
        double x = a * (cos(E) - e);
        double z = a * sqrt(1 - square(e)) * -sin(E);

        // Orientation of the orbital plane with respect to the Laplacian plane
        Mat3d Rorbit     = (Mat3d::yrotation(node) *
                            Mat3d::xrotation(degToRad(i)) *
                            Mat3d::yrotation(w));

        // Rotate to the Earth's equatorial plane
        double N = degToRad(refplane_RA);
        double J = degToRad(90 - refplane_Dec);
        Mat3d RLaplacian = (Mat3d::yrotation( N) *
                            Mat3d::xrotation( J) *
                            Mat3d::yrotation(-N));

        // Rotate to the Martian equatorial plane
        N = degToRad(marspole_RA);
        J = degToRad(90 - marspole_Dec);
        Mat3d RMars_eq   = (Mat3d::yrotation( N) *
                            Mat3d::xrotation(-J) *
                            Mat3d::yrotation(-N));

        Point3d p = Rorbit * Point3d(x, 0.0, z);

        p = RLaplacian * p;

        p = RMars_eq * p;

        return p;
    }


Hopefully, the code is close enough to mathematical notation that it makes sense even if you're not familiar with C++. Can anyone see what I'm missing here? The precession of the node and pericenter seem wrong to me--they're causing Deimos to move slightly faster than it does in my Horizons-derived sampled orbit.

--Chris

[Spaceman Spiff]

Hi Chris,

in such a rush, it's hard to be sure this is right...

I think the error comes from effectively including the precession of pericentre twice. You work out the updated argument of pericentre in the formula

Code: Select all

        double w    = degToRad(w0 + 360.0 * T / Pw);


but then you rotate the orbit with what I guess are Euler angles, the third of which is the argument of pericentre, w.

The Mean Anomaly at epoch M0 is measured around from the argument of pericentre at epoch w0. If you use the new the pericentre w and don't trim it from M0, then when you compute M and then rotate by w you're pushing Phobos and Deimos around at twice the rate you should.

Try working out w first, then work out M as

Code: Select all

       double M    = degToRad((M0  - w + w0) + n * t);


(I think!) or rotate by w0, not w (easier?).

By the way, do you not keep track of units in your code, in comments or appending say _deg, or _year to variable names? I'm guessing: a is in km; w0, M0, i, node0 are degrees, Pw and Pnode are years (for one revolution of pericentre and node respectively). I understand n to be degrees per day. So, it must be Deimos you've showed, as Phobos orbits Mars in 7 hours, and Deimos in 30 hours.

I hope it's easy to see if this fixes it - I've no time to really think deeply about it. Also, I won't be able to correspond until Sat earliest - I'm away on mission.

G'luck.

Ch'rs,

Spiff.

[chris]

Spiff,

That was it! Or close to it . . .

I ended up calculating the mean anomaly as:

M = M0 + t * n - T * dw

and the argument of pericenter as:

w = w0 + T * dw - T * dnode

The orbits are still a bit off though; the positions of both Phobos and Deimos are about 10.5 minutes (in time) behind their Horizons locations. I think that the elements on that web page may have a light time delay included.

--Chris

[Spaceman Spiff]

Ah, yes! The argument of periwotsit is measured around from the ascending node, so that needs to be trimmed out too. The dw is new, but I guess it's Pw divided per year? Ditto dnode...? So, I could agree with those equations.

Yet Phobos and Deimos both appear to be running 10.5 mins late? Yes, very annoying - no timely transits at this rate! Plus 10.5 mins late is 8.23° behind for Phobos, and 2.08° for Deimos - visibly out of place.

Hmmm...

Some things that bother me:

1. Generally, any quantity expressed in years should mean tropical years. This is 365.2422 days, not the 365.25 days you divide by. The precession rates could then be slightly wrong. It's a potential source of error, but I work out it's accumulated only 0.023° for Deimos, and 0.553° for Phobos. Still not enough to account for the above effect.

2. I looked at that Horizons web site. Yes, it appears light time is included in the Horizons data - you see the moons where they were when the light left them. Webpage http://ssd.jpl.nasa.gov/eph_help.html says 'corrected for light-time'. However, today (17 Mar 2004), Mars is 1.8 AU away, and light delay would be more like 15 mins. Check what date you compare positions for.

3. What celestia data do you compare with Horizons - it only gives R.A. and Dec from a place on Earth, not orbit-centred angles, doesn't it? Do you make Celestia output R.A. and Dec and compare to Horizons?

4. I'm surprised the Horizons data is so imprecise for the orbital periods: with several hundred thousand orbits of Phobos since discovery, you'd have thought we get more than three decimal places. Still, the orbit periods are to six decimal places in the solarsys.scc file - but, how do we know those match the Horizons data? Or do you set them to 0.319 and 1.262 exactly in the solarsys.scc file?

5. Try advancing decades into the future or past, and seeing if the delay accumulates (wrong factors) or just oscillates (light delay variation).

Sorry, can't quite crack this one, ol' chap. Hints only until further inspiration turns up!

Spiff.