Saturday, October 25, 2014

Saturday: This is apparently why I write all my own programs.

Because computer programmers are basically awful at everything.

Let's start at the beginning.  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.
And this is where things become me bashing my head against fucking programmers.  I did the initial coding in octave, as after prototyping it, I could simply convert my scalar values to vectors and use the element-wise math in there to parallel a lot of orbits together.  Plus: built in linear algebra and function solvers.  Great.

Except generating a plot in octave in 3d is apparently ruled by the whims of lunatics.  The manual says you can call the axis() function to define the axes of the system, and if you specify "manual," that's retained.  However, every call to scatter3() seemed to reset to an auto set of axes, so even though I generated the thousand frames for my movie, since each frame had unique axes, it looks like shit.  No amount of dicking with the axes elements could fix this.

Ok.  So let's try perl.  Except Gtk is fucking stupid.  Just brain-numbingly "no, you have to do things this way, because there is no other way to do it."  Everything has to have a handler, and a caller, and they all have to be manually packed into other objects, and even then, getting an object you can draw to and clear...  Whatever.  Gtk is broken by design, and I hope the people who wrote it feel bad.

Googling points me to the GD library for perl.  Great.  Allocate an image object of size blah, fill with color whatever, draw to it in arbitrary fashions as you'd expect.  Wonderful.  And it has a 3d mode?  Excellent.

Except 3d mode can do lines, and curves, and polygons, but not a fucking point.  You can't just render a point in 3d space.

Fine, we'll do step five above, because it's functionally equivalent to the R_Omega transform.  You project to a new axis, and then dump the data that's not perpendicular to that axis, as it's just the depth information that you can't see in a 2d figure.  Not that bad of an issue.  

Let's check with a quick test to make sure things work.  There's even an animated gif output product so I can make the movie of the orbits directly.  But, because captain computer decided to increment the version from 2.0.xxx to 2.1.yyy such that yyy < xxx, and then included a check in the makefile that animated gif support is only available when xxx > 33, this doesn't work.  There's a bug report on it in ubuntu, but it hasn't been propagated to debian, so maybe install it from source locally, and fix the bug yourself, and then.

No.  Fuck it.  It shouldn't be my job to fix things because every damn computer programmer elsewhere can't adequately test their shit.  I'm reasonably sure my math is fine, but without a decent data visualization code, I can't check.

So I went and ate a burger and had some cheese tots.

No comments:

Post a Comment