Part of my summer of code project involves rewriting the KStars code so that it uses 3D cartesian coordinates to represent points on a sphere, instead of using spherical coordinates and having to do spherical trigonometry. I talked about this a little bit in my last update, and mentioned that the code to compute stellar aberration was still using spherical coordinates. It’s now been replaced by a new algorithm that uses stereographic projection.
Consider an observer on earth, looking at some star. Relative to the earth, the observer is stationary, but relative to the solar system the observer is moving at the same speed as the earth. The Earth is moving quite quickly around the sun, and this speed is large enough that we can see the effects of relativity: the angle of the beam changes, because the speed of light is constant. (Imagine how rain appears to fall diagonally when travelling in a car, but here instead of just adding the velocities, we use relativity, since the speed of light is constant).
The consequence is that points in the sky appear to move in ellipses as the year passes. Points near the ecliptic plane will travel in flattened ellipses, moving back and forth in a line, while points at the poles will travel in circles.
In the existing implementation, every time we want to compute the position of a point, we use four parameters to estimate the effect of abberation, and then do a lengthy calculation along these lines:
double dRA = -1.0 * K * ( cosRA * cosL * cosOb + sinRA * sinL )/cosDec
+ e * K * ( cosRA * cosP * cosOb + sinRA * sinP )/cosDec;
double dDec = -1.0 * K * ( cosL * cosOb * ( tanOb * cosDec - sinRA * sinDec ) + cosRA * sinDec * sinL )
+ e * K * ( cosP * cosOb * ( tanOb * cosDec - sinRA * sinDec ) + cosRA * sinDec * sinP );This requires a ton of trigonometry calls, it’s difficult to understand what all the terms are, and there’s no work being shared between calculations.
The new algorithm uses a stereographic projection to do a lot of the work. If you don’t know what a stereographic projection is, it’s very simple: to project a point from the sphere onto the plane, you simply draw the line that passes through that point and the north pole, and look at where that line intersects the plane. The only snag is that the north pole itself gets sent to infinity, and though this isn’t a big problem mathematically (we can work in projective space) it’s not good for computations, so we want to avoid doing that. Wikipedia has more information.
My new implementation is roughly along the lines of this paper, except for some minor details. Geometrically, the effect of aberration is to shift the position of a point towards the direction of motion so that \[ \tan \frac{\theta'}{2} = \sqrt{\frac{c-v}{c+v}}\tan\frac{\theta}{2}, \] where \(\theta\) and \(\theta'\) are respectively the true and apparent angles between the point and the direction of motion. (Note that MathJaX isn’t rendered on PlanetKDE or RSS, if the formulas are not displayed properly). It turns out that if the direction of motion is aligned with the south pole, then after doing a stereographic projection, this effect is just scaling by \[ \sqrt{\frac{c-v}{c+v}}. \]
Thus we can compute aberration by rotating our coordinate system to align it with the Earth’s motion, projecting, scaling, and deprojecting. However, KStars already has an implementation of a very accurate method for computing the motion of the Earth at any point in time, so this is even more accurate, because we’re doing the full relativistic calculation instead of just using a special-case method.
What’s more, we don’t need to do any trigonometry at all, just a few simple multiplications and divisions. And unlike the old method, we can share work between projections, computing the velocity and scaling factor just once.
The only worry is that we want to avoid the singularity, but this too is not a big deal, since we can project through the south pole and multiply by the reciprocal. Moreover, since we’re branching based on the location of the point, if we order our data to have spatial locality, then the branch is predictable and basically goes away. (I’m saving more details on performance changes in KStars for another post).
Before, we were using spherical coordinates, and the expression to get the change caused by aberration was really awkward. But here, we’re using a coordinate system which fits the problem, and describing the change is simple: it’s just a scaling factor!
It’s a really nice example of how choosing the right language to use to describe the problem you’re having lets you get a much better solution.