Orbital mechanics

Julie mentioned a hack-mo-tron thing about asteroids last week.  I was reasonably sure that the math wasn't super difficult, and so after waking up today, I decided to sit down and see if that was true or not.  First, you need to see the diagram of things from wikipedia:

Then you need to note that the orbit is elliptical, and finally remember Kepler's second law.

And then you do math at it.

It's harder than I thought, mostly because step two doesn't have a analytic solution.  The methodology here is as follows.

1. First, work in a coordinate system in which the origin is at one focus of the elliptical orbit, and write down the equation for the distance to the orbit from that origin as a function of the true anomaly, listed here as theta.  Precalculate the integral of that because you need it later.
2. Next, we want to work in fixed (or at least known) time steps.  Therefore we need to know theta as a function of t.  The t here is actually (t_abs - t_perihelion) / period, in order to normalize it sanely.  Next, you do a bunch of calculus, and use the fact that the full area is swept out over one period to calibrate things.  This gives a mess that's effectively a constant for a given t, called kappa here.  This has to be solved analytically, but over the range [-pi:pi], there's only one non-trivial solution, so you can do a quick Newton method solution to find it.
3. Use the anomaly and orbital radius to write the vector U, which contains the (u,v,w) position in orbit coordinates.
4. Use R_omega to orient the orbit in real space using omega, the argument of perihelion.  This is just a simple spin.  Next, use R_Omega to handle the inclination, with Omega the ascending node and i the inclination angle.  This yields X, which is the (x,y,z) in a universal coordinate system.
5. Deproject (x,y,z) to (x',y'), because there is no decent 3d visualization code.

I don't have plots because I couldn't find a good arbitrary plotting library that was easy enough to use to do so.