KStars GSoC: Progress Update

Here’s a brief update on the progress with KStars. The main accomplishment so far is writing new implementations of the coordinate conversions that KStars uses to calculate the positions of sky objects. The new code uses linear algebra instead of spherical trigonometry, so we represent all of our points as vectors of length one, and our coordinate conversions become linear maps.

In the old code, everything was done with spherical trigonometry, which means that doing any calculation requires many calls to sin, cos, tan and friends, and also causes a lot of problems because all of the coordinate systems have singularities at the north and south poles. The biggest disadvantage, however, is that it means that it’s impossible to seperate calculating what the transformation is from actually applying the transformation.

For example, consider the old implementation of nutation:

void SkyPoint::nutate(const KSNumbers *num) {
    double cosRA, sinRA, cosDec, sinDec, tanDec;
    double cosOb, sinOb;

    RA.SinCos( sinRA, cosRA );
    Dec.SinCos( sinDec, cosDec );

    num->obliquity()->SinCos( sinOb, cosOb );

    //Step 2: Nutation
    if ( fabs( Dec.Degrees() ) < 80.0 ) { //approximate method
        tanDec = sinDec/cosDec;

        double dRA  = num->dEcLong()*( cosOb + sinOb*sinRA*tanDec ) - num->dObliq()*cosRA*tanDec;
        double dDec = num->dEcLong()*( sinOb*cosRA ) + num->dObliq()*sinRA;

        RA.setD( RA.Degrees() + dRA );
        Dec.setD( Dec.Degrees() + dDec );
    } else { //exact method
        dms EcLong, EcLat;
        findEcliptic( num->obliquity(), EcLong, EcLat );

        //Add dEcLong to the Ecliptic Longitude
        dms newLong( EcLong.Degrees() + num->dEcLong() );
        setFromEcliptic( num->obliquity(), newLong, EcLat );
    }
}

There are a few things to notice here:

1. Because of the overuse of trig functions, we introduce a lot of useless variables just so that we can use some GNU extension to compute sin and cos at the same time for a slight speed boost, costing us in readability.

2. Again because of speed considerations, we need to use an approximate method instead of an exact method when we can get away with it.

3. This method doesn’t actually take the date we want to “nutate to” as a parameter. Instead, you have to pass in a KSNumbers class, which is basically a huge mess of unrelated variables that depend on time and are used for various computations.

Here is the new implementation:

namespace Convert {
...
CoordConversion Nutate(const JulianDate jd)
{
    double dEcLong, dObliq;
    AstroVars::nutationVars(jd, &dEcLong, &dObliq);
    //Add dEcLong to the Ecliptic Longitude
    CoordConversion rot = AngleAxisd(dEcLong*DEG2RAD,Vector3d::UnitY()).matrix();
    return EclToEq(jd) * rot * EqToEcl(jd);
}
...
}

which gives a matrix that can be used as follows:

JulianDate jd = ...;
EquatorialCoord point = ...;
EquatorialCoord nutated = Convert::Nutate(jd) * point;

Notice:

1. When we want to compute the nutation for a particular date, we just give that date, and not a huge bundle of numbers. The particular numbers that we need are in their own function, nutationVars.

2. It’s roughly the same amount of work to compute Convert::nutationVars as to compute one nutated point with the exact method. But then for every point, we just have to multiply a vector by a matrix – 9 multiplications and 6 additions – instead of having to redo all the work.

3. Because the matrix is computed only once, there’s no additional cost to using the exact method instead of the approximation.

So, by doing things in a more elegant way, we gain readability, cut the size of the codebase, improve speed, and allow for further optimization by batch processing.

So far, I’ve done this for all of the coordinate conversions, except for the computation of aberration, which has a new implementation but is still trig-based, since unlike the other computations, aberration is not an orthogonal map. It should be possible to optimize it more, though. Also, all of the new code has unit tests, so that we can check whether or not it is correct, unlike the old code which has neither tests nor a coherent idea of what “correct behaviour” might be.

The next step is to set up code that stores and works with arrays of objects, so that we can actually do batch processing.

Posted July 9, 2013 under kstars, code, planetkde.